From 89149ea7eb44b2f83204b1486ecf87aefc489077 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Thu, 13 Oct 2016 14:38:06 -0400 Subject: [PATCH] added treatment of delayed neutrons to the multi-group solver --- docs/source/io_formats/mgxs_library.rst | 12 +- docs/source/usersguide/input.rst | 12 +- .../python/pincell_multigroup/build-xml.py | 6 +- openmc/mgxs/library.py | 79 +- openmc/mgxs_library.py | 934 ++++++-- src/cmfd_input.F90 | 2 +- src/global.F90 | 5 +- src/hdf5_interface.F90 | 21 + src/input_xml.F90 | 58 +- src/mgxs_data.F90 | 88 +- src/mgxs_header.F90 | 1975 +++++++++++++---- src/physics_mg.F90 | 26 +- src/source.F90 | 6 +- src/tally.F90 | 596 ++++- src/tally_filter.F90 | 4 +- 15 files changed, 3019 insertions(+), 805 deletions(-) diff --git a/docs/source/io_formats/mgxs_library.rst b/docs/source/io_formats/mgxs_library.rst index 04474c01b..a66764a83 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 MeV. The length of this array should be the number of groups plus 1. @@ -110,7 +112,7 @@ Temperature-dependent data, provided for temperature K. 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 + - **inverse-velocity** (*double[]*) -- Average inverse velocity for each of the groups in the library. This dataset is optional. **//K/scatter data/** @@ -159,4 +161,4 @@ Data specific to neutron scattering for the temperature K The former is used when `representation` is "isotropic", and the latter when `representation` is "angle". Finally, the g_out dimension has a dimensionality of `g_out bounds`[0] to - `g_out bounds`[1]. \ No newline at end of file + `g_out bounds`[1]. diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 9e3211819..a7d38f01c 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1681,8 +1681,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 @@ -1863,12 +1861,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. | +----------------------+---------------------------------------------------+ @@ -1897,8 +1893,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 | @@ -1934,8 +1928,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/python/pincell_multigroup/build-xml.py b/examples/python/pincell_multigroup/build-xml.py index 216c56539..31a6c52f7 100644 --- a/examples/python/pincell_multigroup/build-xml.py +++ b/examples/python/pincell_multigroup/build-xml.py @@ -19,7 +19,7 @@ groups = openmc.mgxs.EnergyGroups(group_edges=[1E-11, 0.0635E-6, 10.0E-6, 1.0E-4, 1.0E-3, 0.5, 1.0, 20.0]) # Instantiate the 7-group (C5G7) cross section data -uo2_xsdata = openmc.XSdata('UO2', groups) +uo2_xsdata = openmc.XSdata('UO2', groups, delayed_groups=8) uo2_xsdata.order = 0 uo2_xsdata.set_total( [0.1779492, 0.3298048, 0.4803882, 0.5543674, 0.3118013, 0.3951678, @@ -42,6 +42,8 @@ uo2_xsdata.set_nu_fission([2.005998E-02, 2.027303E-03, 1.570599E-02, 5.257105E-01]) uo2_xsdata.set_chi([5.8791E-01, 4.1176E-01, 3.3906E-04, 1.1761E-07, 0.0000E+00, 0.0000E+00, 0.0000E+00]) +uo2_xsdata.set_beta([2.13333e-4, 1.04514e-3, 6.03969e-4, 1.33963e-3, + 2.29386e-3, 7.05174e-4, 6.00381e-4, 2.07736e-4]) h2o_xsdata = openmc.XSdata('LWTR', groups) h2o_xsdata.order = 0 @@ -59,7 +61,7 @@ h2o_xsdata.set_scatter_matrix( [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0022157, 0.6999130, 0.5373200], [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.1324400, 2.4807000]]]) -mg_cross_sections_file = openmc.MGXSLibrary(groups) +mg_cross_sections_file = openmc.MGXSLibrary(groups, delayed_groups=8) mg_cross_sections_file.add_xsdatas([uo2_xsdata, h2o_xsdata]) mg_cross_sections_file.export_to_hdf5() diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index 02c63f50b..e91f3a245 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -898,6 +898,7 @@ class Library(object): if nuclide is not 'total': name += '_' + nuclide xsdata = openmc.XSdata(name, self.energy_groups) + xsdata.delayed_groups = self.num_delayed_groups if order is None: # Set the order to the Library's order (the defualt behavior) @@ -924,41 +925,79 @@ class Library(object): mymgxs = self.get_mgxs(domain, 'nu-transport') xsdata.set_total_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide], subdomains=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') @@ -966,6 +1005,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))): @@ -976,6 +1016,7 @@ class Library(object): nuclide=[nuclide], subdomain=subdomain) using_multiplicity = True + else: using_multiplicity = False @@ -991,14 +1032,13 @@ class Library(object): nuclide=[nuclide], subdomain=subdomain) - # Since we are not using multiplicity, then - # scattering multiplication (nu-scatter) must be - # accounted for approximately by using an adjusted - # absorption cross section. + # Since we are not using multiplicity, then scattering + # multiplication (nu-scatter) must be accounted for + # approximately by using an adjusted absorption cross section. if 'total' 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 +1091,8 @@ class Library(object): xs_type = 'macro' # Initialize file - mgxs_file = openmc.MGXSLibrary(self.energy_groups) + mgxs_file = openmc.MGXSLibrary(self.energy_groups, + 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_library.py b/openmc/mgxs_library.py index cba5c20a1..c661bb52f 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -34,6 +34,8 @@ class XSdata(object): temperatures : numpy.ndarray Temperatures (in units of Kelvin) of the provided datasets. Defaults to a single temperature at 294K. + delayed_groups : int + Number of delayed groups Attributes ---------- @@ -47,6 +49,8 @@ class XSdata(object): to a single temperature at 294K. energy_groups : openmc.mgxs.EnergyGroups Energy group structure + delayed_groups : int + Num delayed groups fissionable : bool Whether or not this is a fissionable data set. scatter_format : {'legendre', 'histogram', or 'tabular'} @@ -64,20 +68,6 @@ class XSdata(object): num_polar : int Number of equal width angular bins that the polar angular domain is subdivided into. This only applies when ``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 - ``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 ``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 ``representation``. total : dict of numpy.ndarray Group-wise total cross section. absorption : dict of numpy.ndarray @@ -102,37 +92,37 @@ 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 ``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 ``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. - - Notes - ----- - The parameters containing cross section data have dimensionalities which - depend upon the value of ``representation`` as well as the number of - Legendre or other angular dimensions as described by ``order``. The - ``vector_shape``, ``matrix_shape``, and ``pn_matrix_shape`` properties are - provided to obtain the dimensionality of the data for each temperature. - - The following are cross sections which should use each of these properties: - - 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 + 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 cross section vector. + inverse_velocity : dict of numpy.ndarray + Inverse of velocity, in units of sec/cm. """ def __init__(self, name, energy_groups, temperatures=[294.], - representation='isotropic'): + representation='isotropic', delayed_groups=0): + # Initialize class attributes self.name = name self.energy_groups = energy_groups + self.delayed_groups = delayed_groups self.temperatures = temperatures self.representation = representation self._awr = None @@ -141,16 +131,21 @@ class XSdata(object): 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] @property def name(self): @@ -160,6 +155,10 @@ class XSdata(object): def energy_groups(self): return self._energy_groups + @property + def delayed_groups(self): + return self._delayed_groups + @property def representation(self): return self._representation @@ -192,10 +191,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 @@ -220,6 +215,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 @@ -228,6 +231,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: @@ -236,34 +247,6 @@ class XSdata(object): else: 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) - - @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) - - @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) - @name.setter def name(self, name): check_type('name for XSdata', name, basestring) @@ -271,6 +254,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) @@ -278,16 +262,28 @@ class XSdata(object): msg = 'Unable to assign an EnergyGroups object ' \ 'with uninitialized group edges' raise ValueError(msg) + self._energy_groups = energy_groups + @delayed_groups.setter + def delayed_groups(self, delayed_groups): + + # Check validity of delayed_groups + check_type('delayed_groups', delayed_groups, int) + check_greater_than('delayed_groups', delayed_groups, 0, equality=True) + + self._delayed_groups = delayed_groups + @representation.setter def representation(self, representation): + # Check it is of valid type. check_value('representation', representation, _REPRESENTATIONS) self._representation = representation @awr.setter def awr(self, awr): + # Check validity of type and that the awr value is > 0 check_type('awr', awr, Real) check_greater_than('awr', awr, 0.0) @@ -300,8 +296,10 @@ class XSdata(object): @temperatures.setter def temperatures(self, temperatures): + check_type('temperatures', temperatures, Iterable, expected_iter_type=Real) + # Convert to a numpy array so we can easily get the shape for checking nptemperatures = np.asarray(temperatures) @@ -311,12 +309,14 @@ class XSdata(object): @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 @order.setter def order(self, order): + # Check type and value check_type('order', order, Integral) check_greater_than('order', order, 0, equality=True) @@ -324,6 +324,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) @@ -331,15 +332,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 @@ -351,6 +348,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) @@ -362,9 +360,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 @@ -383,10 +387,19 @@ 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 + if self.representation is 'isotropic': + shapes = [(self.energy_groups.num_groups,)] + else: + shapes = [(self.num_polar, self.num_azimuthal, + self.energy_groups.num_groups)] + # Convert to a numpy array so we can easily get the shape for checking nptotal = np.asarray(total) - check_value('total shape', nptotal.shape, [self.vector_shape]) + check_value('total shape', nptotal.shape, shapes) check_type('temperature', temperature, Real) check_value('temperature', temperature, self.temperatures) @@ -410,11 +423,19 @@ 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 + if self.representation is 'isotropic': + shapes = [(self.energy_groups.num_groups,)] + else: + shapes = [(self.num_polar, self.num_azimuthal, + self.energy_groups.num_groups)] + # Convert to a numpy array so we can easily get the shape for checking npabsorption = np.asarray(absorption) - check_value('absorption shape', npabsorption.shape, - [self.vector_shape]) + check_value('absorption shape', npabsorption.shape, shapes) check_type('temperature', temperature, Real) check_value('temperature', temperature, self.temperatures) @@ -438,10 +459,19 @@ 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 + if self.representation is 'isotropic': + shapes = [(self.energy_groups.num_groups,)] + else: + shapes = [(self.num_polar, self.num_azimuthal, + self.energy_groups.num_groups)] + # Convert to a numpy array so we can easily get the shape for checking npfission = np.asarray(fission) - check_value('fission shape', npfission.shape, [self.vector_shape]) + check_value('fission shape', npfission.shape, shapes) check_type('temperature', temperature, Real) check_value('temperature', temperature, self.temperatures) @@ -468,12 +498,20 @@ 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 + if self.representation is 'isotropic': + shapes = [(self.energy_groups.num_groups,)] + else: + shapes = [(self.num_polar, self.num_azimuthal, + self.energy_groups.num_groups)] + # Convert to a numpy array so we can easily get the shape for checking npkappa_fission = np.asarray(kappa_fission) - check_value('kappa fission shape', npkappa_fission.shape, - [self.vector_shape]) + check_value('kappa fission shape', npkappa_fission.shape, shapes) check_type('temperature', temperature, Real) check_value('temperature', temperature, self.temperatures) @@ -500,27 +538,167 @@ 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 + if self.representation is 'isotropic': + shapes = [(self.energy_groups.num_groups,)] + else: + shapes = [(self.num_polar, self.num_azimuthal, + self.energy_groups.num_groups)] + # Convert to a numpy array so we can easily get the shape for checking npchi = np.asarray(chi) - # Check the shape - if npchi.shape != self.vector_shape: - msg = 'Provided chi iterable does not have the expected shape.' - raise ValueError(msg) + check_value('chi shape', npchi.shape, shapes) check_type('temperature', temperature, Real) check_value('temperature', temperature, self.temperatures) i = self.temperatures.tolist().index(temperature) self._chi[i] = npchi - 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 + if self.representation is 'isotropic': + shapes = [(self.energy_groups.num_groups,)] + else: + shapes = [(self.num_polar, self.num_azimuthal, + self.energy_groups.num_groups)] + + # Convert to a numpy array so we can easily get the shape for checking + npchi_prompt = np.asarray(chi_prompt) + check_value('chi prompt shape', npchi_prompt.shape, shapes) + check_type('temperature', temperature, Real) + check_value('temperature', temperature, self.temperatures) + + i = self.temperatures.tolist().index(temperature) + self._chi_prompt[i] = npchi_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 + if self.representation is 'isotropic': + shapes = [(self.energy_groups.num_groups,), + (self.delayed_groups, self.energy_groups.num_groups)] + else: + shapes = [(self.num_polar, self.num_azimuthal, + self.energy_groups.num_groups), + (self.delayed_groups, self.num_polar, self.num_azimuthal, + self.energy_groups.num_groups)] + + # Convert to a numpy array so we can easily get the shape for checking + npchi_delayed = np.asarray(chi_delayed) + check_value('chi delayed shape', npchi_delayed.shape, shapes) + check_type('temperature', temperature, Real) + check_value('temperature', temperature, self.temperatures) + + i = self.temperatures.tolist().index(temperature) + self._chi_delayed[i] = npchi_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 + if self.representation is 'isotropic': + shapes = [(self.delayed_groups,), + (self.delayed_groups, self.energy_groups.num_groups)] + else: + shapes = [(self.delayed_groups, self.num_polar, self.num_azimuthal, + self.energy_groups.num_groups), + (self.delayed_groups, self.num_polar, self.num_azimuthal)] + + # Convert to a numpy array so we can easily get the shape for checking + npbeta = np.asarray(beta) + check_value('beta shape', npbeta.shape, shapes) + check_type('temperature', temperature, Real) + check_value('temperature', temperature, self.temperatures) + + i = self.temperatures.tolist().index(temperature) + self._beta[i] = npbeta + + 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() + + """ + + check_type('decay_rate', decay_rate, Iterable, expected_iter_type=Real) + + # Get the accepted shapes for this xs + if self.representation is 'isotropic': + shapes = [(self.delayed_groups,), + (self.delayed_groups, self.energy_groups.num_groups)] + else: + shapes = [(self.delayed_groups, self.num_polar, self.num_azimuthal, + self.energy_groups.num_groups), + (self.delayed_groups, self.num_polar, self.num_azimuthal)] + + # Convert to a numpy array so we can easily get the shape for checking + npdecay_rate = np.asarray(decay_rate) + check_value('decay rate shape', npdecay_rate.shape, shapes) + check_type('temperature', temperature, Real) + check_value('temperature', temperature, self.temperatures) + + i = self.temperatures.tolist().index(temperature) + self._decay_rate[i] = npdecay_rate def set_scatter_matrix(self, scatter, temperature=294.): """This method sets the cross section for this XSdata object at the @@ -539,11 +717,21 @@ class XSdata(object): openmc.mgxs_library.set_scatter_matrix_mgxs() """ + + # Get the accepted shapes for this xs + if self.representation is 'isotropic': + shapes = [(self.num_orders, self.energy_groups.num_groups, + self.energy_groups.num_groups)] + else: + shapes = [(self.num_polar, self.num_azimuthal, self.num_orders, + self.energy_groups.num_groups, + self.energy_groups.num_groups)] + # Convert to a numpy array so we can easily get the shape for checking npscatter = np.asarray(scatter) check_iterable_type('scatter', npscatter, Real, max_depth=len(npscatter.shape)) - check_value('scatter shape', npscatter.shape, [self.pn_matrix_shape]) + check_value('scatter shape', npscatter.shape, shapes) check_type('temperature', temperature, Real) check_value('temperature', temperature, self.temperatures) @@ -567,12 +755,21 @@ class XSdata(object): openmc.mgxs_library.set_multiplicity_matrix_mgxs() """ + + # Get the accepted shapes for this xs + if self.representation is 'isotropic': + shapes = [(self.energy_groups.num_groups, + self.energy_groups.num_groups)] + else: + shapes = [(self.num_polar, self.num_azimuthal, + self.energy_groups.num_groups, + self.energy_groups.num_groups)] + # Convert to a numpy array so we can easily get the shape for checking npmultiplicity = np.asarray(multiplicity) check_iterable_type('multiplicity', npmultiplicity, Real, max_depth=len(npmultiplicity.shape)) - check_value('multiplicity shape', npmultiplicity.shape, - [self.matrix_shape]) + check_value('multiplicity shape', npmultiplicity.shape, shapes) check_type('temperature', temperature, Real) check_value('temperature', temperature, self.temperatures) @@ -596,82 +793,153 @@ 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 + if self.representation is 'isotropic': + shapes = [(self.energy_groups.num_groups,), + (self.energy_groups.num_groups, + self.energy_groups.num_groups)] + else: + shapes = [(self.num_polar, self.num_azimuthal, + self.energy_groups.num_groups), + (self.num_polar, self.num_azimuthal, + self.energy_groups.num_groups, + self.energy_groups.num_groups)] # Convert to a numpy array so we can easily get the shape for checking npnu_fission = np.asarray(nu_fission) - + check_value('nu_fission shape', npnu_fission.shape, shapes) check_iterable_type('nu_fission', npnu_fission, Real, max_depth=len(npnu_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', npnu_fission.shape, - [self.vector_shape]) - else: - check_value('nu_fission shape', npnu_fission.shape, - [self.matrix_shape]) - else: - check_value('nu_fission shape', npnu_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 npnu_fission.shape == self.vector_shape: - self.use_chi = True - else: - self.use_chi = False - i = self.temperatures.tolist().index(temperature) self._nu_fission[i] = npnu_fission if np.sum(npnu_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 + if self.representation is 'isotropic': + shapes = [(self.energy_groups.num_groups,)] + else: + shapes = [(self.num_polar, self.num_azimuthal, + self.energy_groups.num_groups)] + + # Convert to a numpy array so we can easily get the shape for checking + npprompt_nu_fission = np.asarray(prompt_nu_fission) + check_value('prompt_nu_fission shape', npprompt_nu_fission.shape, shapes) + check_iterable_type('prompt_nu_fission', npprompt_nu_fission, Real, + max_depth=len(npprompt_nu_fission.shape)) + check_type('temperature', temperature, Real) + check_value('temperature', temperature, self.temperatures) + + i = self.temperatures.tolist().index(temperature) + self._prompt_nu_fission[i] = npprompt_nu_fission + if np.sum(npprompt_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 + if self.representation is 'isotropic': + shapes = [(self.delayed_groups, self.energy_groups.num_groups,)] + else: + shapes = [(self.delayed_groups, self.num_polar, self.num_azimuthal, + self.energy_groups.num_groups)] + + # Convert to a numpy array so we can easily get the shape for checking + npdelayed_nu_fission = np.asarray(delayed_nu_fission) + check_value('delayed_nu_fission shape', npdelayed_nu_fission.shape, + shapes) + check_iterable_type('delayed_nu_fission', npdelayed_nu_fission, Real, + max_depth=len(npdelayed_nu_fission.shape)) + check_type('temperature', temperature, Real) + check_value('temperature', temperature, self.temperatures) + + i = self.temperatures.tolist().index(temperature) + self._delayed_nu_fission[i] = npdelayed_nu_fission + if np.sum(npdelayed_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 units of Kelvin) of the provided dataset. Defaults to 294K """ - check_type('inverse velocities', inv_vel, Iterable, + + check_type('inverse_velocity', inv_vel, Iterable, expected_iter_type=Real) + + # Get the accepted shapes for this xs + if self.representation is 'isotropic': + shapes = [(self.energy_groups.num_groups,)] + else: + shapes = [(self.num_polar, self.num_azimuthal, + self.energy_groups.num_groups)] + # Convert to a numpy array so we can easily get the shape for checking npinv_vel = np.asarray(inv_vel) - check_value('inverse velocities shape', npinv_vel.shape, - [self.vector_shape]) + check_value('inverse_velocity shape', npinv_vel.shape, shapes) check_type('temperature', temperature, Real) check_value('temperature', temperature, self.temperatures) i = self.temperatures.tolist().index(temperature) - self._inverse_velocities[i] = npinv_vel + self._inverse_velocity[i] = npinv_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 units of Kelvin) of the provided dataset. Defaults to 294K @@ -850,14 +1118,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, + ['universe', 'cell', 'material', 'mesh']) + check_type('temperature', temperature, Real) + check_value('temperature', temperature, self.temperatures) + + i = self.temperatures.tolist().index(temperature) + 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('delayed_groups', delayed_nu_fission.num_delayed_groups, + [self.delayed_groups]) + check_value('domain_type', delayed_nu_fission.domain_type, + ['universe', 'cell', 'material', 'mesh']) + check_type('temperature', temperature, Real) + check_value('temperature', temperature, self.temperatures) + + i = self.temperatures.tolist().index(temperature) + 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): @@ -936,12 +1307,6 @@ 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, @@ -957,8 +1322,191 @@ 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, + ['universe', 'cell', 'material', 'mesh']) + check_type('temperature', temperature, Real) + check_value('temperature', temperature, self.temperatures) + + i = self.temperatures.tolist().index(temperature) + 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('delayed_groups', chi_delayed.num_delayed_groups, + [self.delayed_groups]) + check_value('domain_type', chi_delayed.domain_type, + ['universe', 'cell', 'material', 'mesh']) + check_type('temperature', temperature, Real) + check_value('temperature', temperature, self.temperatures) + + i = self.temperatures.tolist().index(temperature) + 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('delayed_groups', beta.num_delayed_groups, + [self.delayed_groups]) + check_value('domain_type', beta.domain_type, + ['universe', 'cell', 'material', 'mesh']) + check_type('temperature', temperature, Real) + check_value('temperature', temperature, self.temperatures) + + i = self.temperatures.tolist().index(temperature) + 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('delayed_groups', decay_rate.num_delayed_groups, + [self.delayed_groups]) + check_value('domain_type', decay_rate.domain_type, + ['universe', 'cell', 'material', 'mesh']) + check_type('temperature', temperature, Real) + check_value('temperature', temperature, self.temperatures) + + i = self.temperatures.tolist().index(temperature) + 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', @@ -1114,22 +1662,29 @@ class XSdata(object): np.nan_to_num(self._multiplicity_matrix[i]) def _get_xsdata_group(self, file): + grp = file.create_group(self.name) + if self.awr is not None: grp.attrs['awr'] = self.awr + if self.fissionable is not None: grp.attrs['fissionable'] = self.fissionable + if self.representation is not None: grp.attrs['representation'] = np.array(self.representation, dtype='S') 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 + if self.scatter_format is not None: grp.attrs['scatter-format'] = np.array(self.scatter_format, dtype='S') + if self.order is not None: grp.attrs['order'] = self.order @@ -1141,27 +1696,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 ' @@ -1170,11 +1761,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 = [] @@ -1183,12 +1777,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 = [] @@ -1196,6 +1793,7 @@ 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", data=np.array(flat_mult)) @@ -1209,6 +1807,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): @@ -1216,6 +1815,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): @@ -1226,12 +1826,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): @@ -1241,6 +1844,7 @@ 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)) @@ -1252,9 +1856,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): @@ -1266,23 +1870,32 @@ class MGXSLibrary(object): ---------- energy_groups : openmc.mgxs.EnergyGroups Energy group structure + delayed_groups : int + Num delayed groups Attributes ---------- energy_groups : openmc.mgxs.EnergyGroups Energy group structure. + 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, delayed_groups=0): self.energy_groups = energy_groups + self.delayed_groups = delayed_groups self._xsdatas = [] @property def energy_groups(self): return self._energy_groups + @property + def delayed_groups(self): + return self._delayed_groups + @property def temperatures(self): return self._temperatures @@ -1296,6 +1909,11 @@ class MGXSLibrary(object): check_type('energy groups', energy_groups, openmc.mgxs.EnergyGroups) self._energy_groups = energy_groups + @delayed_groups.setter + def delayed_groups(self, delayed_groups): + check_type('delayed groups', delayed_groups, int) + self._delayed_groups = delayed_groups + def add_xsdata(self, xsdata): """Add an XSdata entry to the file. @@ -1305,10 +1923,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) @@ -1324,6 +1944,7 @@ class MGXSLibrary(object): XSdatas to add """ + check_iterable_type('xsdatas', xsdatas, XSdata) for xsdata in xsdatas: @@ -1360,7 +1981,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.delayed_groups file.attrs['group structure'] = self.energy_groups.group_edges for xsdata in self._xsdatas: diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 327b55b42..a60fc9a89 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 1e837b8c0..48105bbf9 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -115,7 +115,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 262bbc232..da8ebe120 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -99,6 +99,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 @@ -2601,6 +2602,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 4702c3bcf..536d94ad9 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3060,8 +3060,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 @@ -3086,8 +3086,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 @@ -3101,14 +3101,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) @@ -3599,7 +3591,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 @@ -3612,12 +3610,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') @@ -4654,15 +4646,25 @@ contains ! Open file for reading file_id = file_open(path_cross_sections, 'r', parallel=.true.) - if (check_attribute(file_id, "groups")) then - ! Get neutron group count - call read_attribute(energy_groups, file_id, "groups") + if (check_attribute(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 (check_attribute(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 + call write_message("WARNING: delayed_groups element not provided so & + &number of delayed groups set to 0") + end if + + allocate(rev_energy_bins(num_energy_groups + 1)) + allocate(energy_bins(num_energy_groups + 1)) + if (check_attribute(file_id, "group structure")) then ! Get neutron group structure call read_attribute(energy_bins, file_id, "group structure") @@ -4671,10 +4673,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 2fbfdfd42..f598958d4 100644 --- a/src/mgxs_data.F90 +++ b/src/mgxs_data.F90 @@ -21,26 +21,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!") @@ -60,23 +59,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 @@ -103,7 +85,9 @@ contains ! First find out the data representation if (check_attribute(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 @@ -118,16 +102,19 @@ 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 @@ -147,6 +134,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. @@ -161,12 +149,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)) @@ -174,27 +162,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 !=============================================================================== @@ -203,16 +196,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 MeV + real(8) :: kT ! temperature in MeV 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 @@ -243,17 +237,18 @@ contains !=============================================================================== subroutine get_temperatures(temps) - type(VectorReal), allocatable, intent(out) :: temps(:) + type(VectorReal), allocatable, intent(out) :: temps(:) integer :: i, j, k integer :: i_nuclide ! index in nuclides array - integer :: i_material + integer :: i_material ! Index in materials array real(8) :: temperature ! temperature in Kelvin allocate(temps(n_nuclides_total)) 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 @@ -268,6 +263,7 @@ contains i_material = material_dict % get_key(cells(i) % material(j)) associate (mat => materials(i_material)) NUC_NAMES_LOOP: do k = 1, size(mat % names) + ! Get index in temps array i_nuclide = nuclide_dict % get_key(to_lower(mat % names(k))) diff --git a/src/mgxs_header.F90 b/src/mgxs_header.F90 index 193f337f1..a26453a8f 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 (Gin x Dg) + 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 !=============================================================================== @@ -66,12 +72,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 +99,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 +138,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 @@ -393,29 +406,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), allocatable :: temp_beta(:, :) real(8) :: dmu, mu, norm 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(:) @@ -427,99 +441,458 @@ module mgxs_header ! 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 data for all the cross sections + allocate(xs % total(energy_groups)) + allocate(xs % absorption(energy_groups)) + allocate(xs % delayed_nu_fission(energy_groups, delayed_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(energy_groups, energy_groups, & + delayed_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 + if (this % fissionable) then - allocate(xs % nu_fission(groups)) - allocate(xs % chi(groups, groups)) - if (check_dataset(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(energy_groups, delayed_groups)) + + ! Set beta + if (check_dataset(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(gin, dg) = 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 (check_dataset(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, (/energy_groups, delayed_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 (check_dataset(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 - - ! Set the vector nu-fission from the matrix nu-fission - do gin = 1, groups - xs % nu_fission(gin) = sum(temp_2d(:, gin)) - 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)) - end do - deallocate(temp_2d) + temp_beta = ZERO 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 (check_dataset(xsdata_grp, "fission")) then - allocate(xs % fission(groups)) - call read_dataset(xs % fission, xsdata_grp, "fission") + ! If chi provided, set chi_prompt and chi_delayed + if (check_dataset(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 + if (sum(xs % chi_prompt(:, gin)) == ZERO) then + call fatal_error("Encountered chi for a group that sums to & + &zero") + else + xs % chi_prompt(:, gin) = xs % chi_prompt(:, gin) / & + sum(xs % chi_prompt(:, gin)) + end if + end do + + ! Set chi_delayed to chi_prompt + do dg = 1, delayed_groups + xs % chi_delayed(:, :, dg) = xs % chi_prompt(:, :) + end do + + ! Deallocate temporary chi array + deallocate(temp_arr) + end if + + ! If nu-fission provided, set prompt_nu_fission and + ! delayed_nu_fission. If nu-fission is a matrix, set chi_prompt and + ! chi_delayed. + if (check_dataset(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(gin, dg) = temp_beta(gin, dg) * & + xs % prompt_nu_fission(gin) + end do + + ! Correct prompt_nu_fission using delayed neutron fraction + xs % prompt_nu_fission(gin) = (1 - sum(temp_beta(gin, :))) * & + xs % prompt_nu_fission(gin) + 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(gin, dg) = temp_beta(gin, dg) * & + xs % prompt_nu_fission(gin) + end do + + ! Correct prompt_nu_fission using delayed neutron fraction + xs % prompt_nu_fission(gin) = (1 - sum(temp_beta(gin, :))) * & + xs % prompt_nu_fission(gin) + end do + + ! Now pull out information needed for chi + xs % chi_prompt(:, :) = temp_2d + + ! Deallocate temporary 2D array for nu_fission matrix + deallocate(temp_2d) + + ! Normalize chi so its CDF goes to 1 + do gin = 1, energy_groups + if (sum(xs % chi_prompt(:, gin)) == ZERO) then + call fatal_error("Encountered chi for a group that sums to & + &zero") + else + xs % chi_prompt(:, gin) = xs % chi_prompt(:, gin) / & + sum(xs % chi_prompt(:, gin)) + 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 (check_dataset(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 (check_dataset(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 + if (sum(xs % chi_prompt(:, gin)) == ZERO) then + call fatal_error("Encountered chi prompt for a group that & + &sums to zero") + else + xs % chi_prompt(:, gin) = xs % chi_prompt(:, gin) / & + sum(xs % chi_prompt(:, gin)) + end if + end do + + ! Deallocate temporary array for chi_prompt + deallocate(temp_arr) + end if + + ! If chi_delayed provided, set chi_delayed + if (check_dataset(xsdata_grp, "chi-delayed")) then + + ! Get the dimensions of the nu-fission 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(gout, gin, dg) = temp_arr(gout) + end do + + ! Normalize chi so its CDF goes to 1 + if (sum(xs % chi_delayed(:, gin, dg)) == ZERO) then + call fatal_error("Encountered chi delayed for a group & + &that sums to zero") + else + xs % chi_delayed(:, gin, dg) = & + xs % chi_delayed(:, gin, dg) / & + sum(xs % chi_delayed(:, gin, dg)) + 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(energy_groups, delayed_groups)) + temp_2d = reshape(temp_arr, (/energy_groups, delayed_groups/)) + + do dg = 1, delayed_groups + do gin = 1, energy_groups + do gout = 1, energy_groups + xs % chi_delayed(gout, gin, dg) = temp_2d(gout, dg) + end do + + ! Normalize chi so its CDF goes to 1 + if (sum(xs % chi_delayed(:, gin, dg)) == ZERO) then + call fatal_error("Encountered chi delayed for a group & + &that sums to zero") + else + xs % chi_delayed(:, gin, dg) = & + xs % chi_delayed(:, gin, dg) / & + sum(xs % chi_delayed(:, gin, dg)) + 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 + + ! If prompt-nu-fission present, set prompt_nu_fission + if (check_dataset(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 (check_dataset(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(gin, dg) = temp_beta(gin, dg) * & + 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) /= energy_groups) then + call fatal_error("The delayed-nu-fission matrix was input & + &with a 1st dimension not equal to the number of & + &energy groups.") + end if + + ! Issue error if 2nd dimension not correct + if (dims(2) /= delayed_groups) then + call fatal_error("The delayed-nu-fission matrix was input & + &with a 2nd dimension not equal to the number of & + &delayed 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 (energy_groups, & + &delayed_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, (/energy_groups, & + delayed_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 (check_dataset(xsdata_grp, "fission")) then + call read_dataset(xs % fission, xsdata_grp, "fission") + end if + + ! Get kappa_fission xs + if (check_dataset(xsdata_grp, "kappa_fission")) then + call read_dataset(xs % kappa_fission, xsdata_grp, "kappa_fission") + end if + + ! Get decay rate xs + if (check_dataset(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 (check_dataset(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 (check_dataset(xsdata_grp, "inverse_velocity")) then + call read_dataset(xs % inverse_velocity, xsdata_grp, & + "inverse_velocity") + else + xs % inverse_velocity = ZERO + end if + ! Get scattering data if (.not. check_group(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 (check_dataset(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 matrix must be provided") end if + if (check_dataset(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 matrix must be provided") @@ -528,11 +901,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. check_dataset(scatt_grp, 'scatter matrix')) & call fatal_error("'scatter matrix' must be provided") call read_dataset(temp_arr, scatt_grp, "scatter matrix") @@ -544,11 +920,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 @@ -557,21 +934,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 @@ -579,11 +963,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 * & @@ -591,8 +978,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) * & @@ -601,33 +990,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 (check_dataset(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") - ! Convert temp_arr to a jagged array ((gin) % data(gout)) for passing - ! to ScattData - allocate(temp_mult(groups)) + ! Convert temp_arr to a jagged array ((gin) % data(gout)) for + ! passing to ScattData + allocate(temp_mult(energy_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 @@ -635,9 +1032,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 @@ -658,12 +1057,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 (check_dataset(xsdata_grp, "total")) then call read_dataset(xs % total, xsdata_grp, "total") else @@ -672,16 +1070,10 @@ 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 (check_dataset(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) @@ -690,28 +1082,29 @@ module mgxs_header 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), allocatable :: temp_1d(:), temp_2d(:, :), temp_3d(:, :, :) + real(8), allocatable :: temp_4d(:, :, :, :), temp_beta(:, :, :, :) real(8) :: dmu, mu, norm - integer :: order, order_dim, gin, gout, l, imu + 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(:, :, :) @@ -725,154 +1118,509 @@ module mgxs_header ! 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)) + + ! Load the more specific data + allocate(xs % prompt_nu_fission(energy_groups, this % n_azi, & + this % n_pol)) + allocate(xs % delayed_nu_fission(energy_groups, delayed_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(energy_groups, energy_groups, & + delayed_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 + if (this % fissionable) then - allocate(xs % nu_fission(groups, this % n_azi, this % n_pol)) - allocate(xs % chi(groups, groups, this % n_azi, this % n_pol)) - if (check_dataset(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(energy_groups, delayed_groups, this % n_azi, & + this % n_pol)) + + ! Set beta + if (check_dataset(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(gin, dg, 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 (check_dataset(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, (/energy_groups, & + delayed_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 (check_dataset(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 (check_dataset(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 + if (sum(xs % chi_prompt(:, gin, iazi, ipol)) == 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) / & + sum(xs % chi_prompt(:, gin, iazi, ipol)) + 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 chi_prompt provided, set chi_prompt + if (check_dataset(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 + if (sum(xs % chi_prompt(:, gin, iazi, ipol)) == 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) / & + sum(xs % chi_prompt(:, gin, iazi, ipol)) + 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 (check_dataset(xsdata_grp, "chi-delayed")) then + + ! Get the dimensions of the nu-fission 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_prompt + allocate(temp_1d(delayed_groups * this % n_azi * this % n_pol)) + allocate(temp_3d(delayed_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, (/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 + do gout = 1, energy_groups + xs % chi_delayed(gout, gin, dg, iazi, ipol) = & + temp_3d(gout, iazi, ipol) + end do + + ! Normalize chi so its CDF goes to 1 + if (sum(xs % chi_delayed(:, gin, dg, iazi, ipol)) & + == ZERO) then + call fatal_error("Encountered chi delayed for a group& + & that sums to zero") + else + xs % chi_delayed(:, gin, dg, iazi, ipol) = & + xs % chi_delayed(:, gin, dg, iazi, ipol) / & + sum(xs % chi_delayed(:, gin, dg, iazi, ipol)) + end if + end do + end do + end do + end do + + ! Deallocate temporary arrays for chi_delayed + deallocate(temp_1d) + deallocate(temp_3d) + + else if (ndims == 2) then + + ! Allocate temporary array for chi_delayed + allocate(temp_1d(delayed_groups * energy_groups * this % n_azi & + * this % n_pol)) + allocate(temp_4d(energy_groups, delayed_groups, this % n_azi, & + this % n_pol)) + + ! Get chi-delayed + call read_dataset(temp_1d, xsdata_grp, "chi-delayed") + temp_4d = reshape(temp_1d, (/energy_groups, delayed_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(gout, gin, dg, iazi, ipol) = & + temp_4d(gout, dg, iazi, ipol) + end do + + ! Normalize chi so its CDF goes to 1 + if (sum(xs % chi_delayed(:, gin, dg, iazi, ipol)) & + == ZERO) then + call fatal_error("Encountered chi delayed for a group& + & that sums to zero") + else + xs % chi_delayed(:, gin, dg, iazi, ipol) = & + xs % chi_delayed(:, gin, dg, iazi, ipol) / & + sum(xs % chi_delayed(:, gin, dg, iazi, ipol)) + 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 (check_dataset(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 (check_dataset(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 == 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 1D & + &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(gin, dg, iazi, ipol) = & + temp_beta(gin, dg, 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) + ! 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) /= energy_groups) then + call fatal_error("The delayed-nu-fission matrix was input & + &with a 1st dimension not equal to the number of & + &energy groups.") + end if + + ! Issue error if 2nd dimension not correct + if (dims(2) /= delayed_groups) then + call fatal_error("The delayed-nu-fission matrix was input & + &with a 2nd dimension not equal to the number of & + &delayed 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 (energy_groups, & + &delayed_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, (/energy_groups, & + delayed_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 (check_dataset(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 (check_dataset(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 (check_dataset(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 (check_dataset(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 (check_dataset(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 end if + ! All the XS unique to fissionable isotopes have been set; Now set all + ! the generation XS + if (check_dataset(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 (check_dataset(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") + + allocate(xs % inverse_velocity(energy_groups, this % n_azi, & + this % n_pol)) + 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. check_group(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 (check_dataset(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 matrix must be provided") end if + if (check_dataset(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 matrix must be provided") @@ -883,17 +1631,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. check_dataset(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. @@ -902,41 +1652,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 @@ -947,12 +1703,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 * & @@ -960,6 +1719,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 @@ -976,7 +1736,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(:, :) = & @@ -985,46 +1745,54 @@ module mgxs_header end do end do end if + deallocate(input_scatt) ! Now get the multiplication matrix if (check_dataset(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 @@ -1035,8 +1803,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) @@ -1058,7 +1828,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 @@ -1066,13 +1836,14 @@ module mgxs_header end do end do - allocate(xs % total(groups, this % n_azi, this % n_pol)) if (check_dataset(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, & + + 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_arr) + + deallocate(temp_1d) else do ipol = 1, this % n_pol do iazi = 1, this % n_azi @@ -1086,7 +1857,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 @@ -1094,20 +1865,11 @@ module mgxs_header end do end do - ! Get kinetics data - if (check_dataset(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 @@ -1134,6 +1896,7 @@ module mgxs_header else this % name = trim(mat % name) end if + this % fissionable = mat % fissionable ! The following info we should initialize, but we dont need it nor @@ -1141,6 +1904,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 @@ -1156,12 +1920,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 @@ -1177,8 +1943,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 @@ -1193,11 +1961,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) @@ -1217,19 +1989,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 + 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) :: groups ! Number of E groups + 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 @@ -1237,6 +2011,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 @@ -1257,31 +2032,54 @@ 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)) + 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(energy_groups, & + delayed_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(energy_groups, energy_groups, & + delayed_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 @@ -1293,16 +2091,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 " // & @@ -1310,11 +2112,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 @@ -1331,7 +2135,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 @@ -1345,25 +2151,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. @@ -1383,12 +2202,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) * & @@ -1404,9 +2226,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) = & @@ -1422,8 +2246,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 @@ -1442,12 +2266,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(:, gin, dg)) + if (norm > ZERO) then + this % xs(t) % chi_delayed(:, gin, dg) = & + this % xs(t) % chi_delayed(:, gin, dg) / norm + end if + end do + end do end if ! Deallocate temporaries @@ -1459,13 +2294,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 @@ -1473,6 +2309,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 @@ -1490,20 +2327,24 @@ module mgxs_header ! Set the meta-data call mgxs_combine(this, temps, mat, nuclides, max_order, scatter_format, & - order_dim) + order_dim) ! 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 @@ -1516,31 +2357,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(energy_groups, & + delayed_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(energy_groups, energy_groups, & + delayed_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 @@ -1562,11 +2427,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 " // & @@ -1574,11 +2442,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 @@ -1598,7 +2468,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 @@ -1607,30 +2479,44 @@ 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. @@ -1652,14 +2538,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 * & @@ -1679,9 +2568,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) @@ -1703,8 +2594,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) / & @@ -1732,15 +2623,31 @@ 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_prompt(:, gin, iazi, ipol)) + if (norm > ZERO) then + this % xs(t) % chi_delayed(:, gin, dg, iazi, ipol) = & + this % xs(t) % chi_delayed(:, gin, dg, iazi, ipol)& + / norm + end if + end do + end do + end do + end do + end if ! Deallocate temporaries @@ -1755,44 +2662,110 @@ 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) - else - xs = ZERO - end if + xs = this % xs(t) % fission(gin) + case('kappa_fission') - if (allocated(this % xs(t) % k_fission)) then - xs = this % xs(t) % k_fission(gin) + 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 = this % xs(t) % decay_rate(1) + end if + + 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(gin, dg) + else + xs = sum(this % xs(t) % delayed_nu_fission(gin, :)) + end if + + case('nu_fission') + xs = this % xs(t) % prompt_nu_fission(gin) + & + sum(this % xs(t) % delayed_nu_fission(gin, :)) + + case('nu') + if (this % xs(t) % fission(gin) > ZERO) then + xs = (sum(this % xs(t) % delayed_nu_fission(gin, :)) + & + this % xs(t) % prompt_nu_fission(gin)) / & + this % xs(t) % fission(gin) else xs = ZERO end if - case('nu_fission') - xs = this % xs(t) % nu_fission(gin) - case('chi') + + case('nu_prompt') + if (this % xs(t) % fission(gin) > ZERO) then + xs = this % xs(t) % prompt_nu_fission(gin) / & + this % xs(t) % fission(gin) + else + xs = ZERO + end if + + case('nu_delayed') + if (this % xs(t) % fission(gin) > ZERO) then + if (present(dg)) then + xs = this % xs(t) % delayed_nu_fission(gin, dg) & + / this % xs(t) % fission(gin) + else + xs = sum(this % xs(t) % delayed_nu_fission(gin, :)) & + / this % xs(t) % fission(gin) + end if + else + xs = ZERO + end if + + 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(gout, gin, dg) + else + xs = this % xs(t) % chi_delayed(gout, gin, 1) + end if + else + if (present(dg)) then + xs = sum(this % xs(t) % chi_delayed(:, gin, dg)) + else + xs = sum(this % xs(t) % chi_delayed(:, gin, 1)) + end if + end if + case('scatter') if (present(gout)) then if (gout < this % xs(t) % scatter % gmin(gin) .or. & @@ -1805,6 +2778,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. & @@ -1820,6 +2794,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. & @@ -1840,20 +2815,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 @@ -1861,33 +2837,100 @@ 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) - else - xs = ZERO - end if + xs = this % xs(t) % fission(gin, iazi, ipol) + case('kappa_fission') - if (allocated(this % xs(t) % k_fission)) then - xs = this % xs(t) % k_fission(gin, iazi, ipol) + 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(gin, dg, iazi, ipol) + else + xs = sum(this % xs(t) % delayed_nu_fission(gin, :, iazi, ipol)) + end if + + case('nu_fission') + xs = this % xs(t) % prompt_nu_fission(gin, iazi, ipol) + & + sum(this % xs(t) % delayed_nu_fission(gin, :, iazi, ipol)) + + case('nu') + if (this % xs(t) % fission(gin, iazi, ipol) > ZERO) then + xs = (sum(this % xs(t) % delayed_nu_fission(gin, :, iazi, ipol)) + & + this % xs(t) % prompt_nu_fission(gin, iazi, ipol)) / & + this % xs(t) % 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_prompt') + if (this % xs(t) % fission(gin, iazi, ipol) > ZERO) then + xs = this % xs(t) % prompt_nu_fission(gin, iazi, ipol) / & + this % xs(t) % fission(gin, iazi, ipol) + else + xs = ZERO + end if + + case('nu_delayed') + if (this % xs(t) % fission(gin, iazi, ipol) > ZERO) then + if (present(dg)) then + xs = this % xs(t) % delayed_nu_fission(gin, dg, iazi, ipol) / & + this % xs(t) % fission(gin, iazi, ipol) + else + xs = sum(this % xs(t) % delayed_nu_fission(gin, :, iazi, ipol)) & + / this % xs(t) % fission(gin, iazi, ipol) + end if + else + xs = ZERO + end if + + 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(gout, gin, dg, iazi, ipol) + else + xs = this % xs(t) % chi_delayed(gout, gin, 1, iazi, ipol) + end if + else + if (present(dg)) then + xs = sum(this % xs(t) % chi_delayed(:, gin, dg, iazi, ipol)) + else + xs = sum(this % xs(t) % chi_delayed(:, gin, 1, 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. & @@ -1900,6 +2943,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. & @@ -1915,6 +2959,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. & @@ -1936,60 +2981,150 @@ 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_pd ! Running probability for prompt/delayed + 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) :: nu, nu_prompt + nu = this % get_xs('nu', gin) + nu_prompt = this % get_xs('nu_prompt', gin) - end function mgxsiso_sample_fission_energy + ! Sample random numbers + xi_pd = prn() + xi_gout = prn() - function mgxsang_sample_fission_energy(this, gin, uvw) result(gout) + ! Neutron is born prompt + if (xi_pd < nu_prompt / nu) then + + ! set the delayed group for the particle born from fission to 0 + dg = 0 + + gout = 1 + prob_gout = this % xs(this % index_temp) % chi_prompt(gout, gin) + + do while (prob_gout < xi_gout) + gout = gout + 1 + prob_gout = prob_gout + & + this % xs(this % index_temp) % chi_prompt(gout, gin) + end do + + ! Neutron is born delayed + else + + ! Get the delayed group + dg = 0 + prob_pd = nu_prompt / nu + + do while (xi_pd >= prob_pd) + dg = dg + 1 + prob_pd = prob_pd + & + this % get_xs('nu_delayed', gin, dg=dg) & + / nu + end do + + ! Adjust dg in case of round off error + dg = min(dg, num_delayed_groups) + + ! Get the outgoing group + gout = 1 + prob_gout = this % xs(this % index_temp) % chi_delayed(gout, gin, dg) + + do while (prob_gout < xi_gout) + gout = gout + 1 + prob_gout = prob_gout + & + this % xs(this % index_temp) % chi_delayed(gout, gin, 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_pd ! Running probability for prompt/delayed + real(8) :: prob_gout ! Running probability for gout + real(8) :: nu, nu_prompt - call find_angle(this % polar, this % azimuthal, uvw, iazi, ipol) + ! Get nu and nu_prompt + nu = this % get_xs('nu', gin, uvw=uvw) + nu_prompt = this % get_xs('nu_prompt', 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 < nu_prompt / nu) 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 + prob_pd = nu_prompt / nu + + do while (xi_pd < prob_pd) + dg = dg + 1 + prob_pd = prob_pd + & + this % get_xs('nu_delayed', gin, uvw=uvw, dg=dg) / nu + end do + + ! Adjust dg in case of round off error + dg = min(dg, 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 = 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 @@ -2038,7 +3173,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 @@ -2060,7 +3197,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 c60fc6e2e..776c7905a 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 9aeccde15..cc85bf6ec 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 24c1b9495..c5a1b8352 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 @@ -1106,6 +1110,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 @@ -1117,7 +1126,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 @@ -1128,16 +1139,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 @@ -1146,6 +1160,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) & + / nucxs % get_xs('total', 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('total', 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 @@ -1441,12 +1473,418 @@ contains 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) / & + nucxs % 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) / nucxs % 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) / nucxs % 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) / nucxs % 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 @@ -1472,10 +1910,10 @@ contains 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 @@ -1931,7 +2369,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 @@ -1950,6 +2388,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) @@ -1976,14 +2415,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. & @@ -2046,91 +2502,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 4dd9bf129..37265cc5f 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.