diff --git a/docs/source/io_formats/nuclear_data.rst b/docs/source/io_formats/nuclear_data.rst index 8d8da8eda..060e96c0c 100644 --- a/docs/source/io_formats/nuclear_data.rst +++ b/docs/source/io_formats/nuclear_data.rst @@ -26,8 +26,9 @@ Incident Neutron Data temperature-dependent data set. For example, the data set corresponding to 300 Kelvin would be located at `300K`. -:Datasets: - **K** (*double*) -- kT values (in MeV) for each Temperature - TTT (in Kelvin) +:Datasets: + - **K** (*double*) -- kT values (in MeV) for each Temperature + TTT (in Kelvin) **//reactions/reaction_/** @@ -44,11 +45,14 @@ temperature-dependent data set. For example, the data set corresponding to temperature-dependent data set. For example, the data set corresponding to 300 Kelvin would be located at `300K`. -:Attributes: - **threshold_idx** (*int*) -- Index on the energy grid that the - reaction threshold corresponds to for temperature TTT (in Kelvin) +:Datasets: + - **xs** (*double[]*) -- Cross section values tabulated against the + nuclide energy grid for temperature TTT (in Kelvin) -:Datasets: - **xs** (*double[]*) -- Cross section values tabulated against the - nuclide energy grid for temperature TTT (in Kelvin) + :Attributes: + - **threshold_idx** (*int*) -- Index on the energy + grid that the reaction threshold corresponds to for + temperature TTT (in Kelvin) **//reactions/reaction_/product_/** @@ -112,7 +116,7 @@ Thermal Neutron Scattering Data **//** :Attributes: - **atomic_weight_ratio** (*double*) -- Mass in units of neutron masses - - **zaids** (*int[]*) -- ZAID identifiers for which the thermal + - **nuclides** (*char[][]*) -- Names of nuclides for which the thermal scattering data applies to - **secondary_mode** (*char[]*) -- Indicates how the inelastic outgoing angle-energy distributions are represented ('equal', @@ -124,8 +128,9 @@ Thermal Neutron Scattering Data temperature-dependent data set. For example, the data set corresponding to 300 Kelvin would be located at `300K`. -:Datasets: - **K** (*double*) -- kT values (in MeV) for each Temperature - TTT (in Kelvin) +:Datasets: + - **K** (*double*) -- kT values (in MeV) for each Temperature + TTT (in Kelvin) **//elastic/K/** diff --git a/docs/source/methods/cross_sections.rst b/docs/source/methods/cross_sections.rst index 4ed777536..a4d0f7d24 100644 --- a/docs/source/methods/cross_sections.rst +++ b/docs/source/methods/cross_sections.rst @@ -53,12 +53,12 @@ speed up the calculation. Logarithmic Mapping +++++++++++++++++++ -To speed up energy grid searches, OpenMC uses logarithmic mapping technique -[Brown]_ to limit the range of energies that must be searched for each -nuclide. The entire energy range is divided up into equal-lethargy segments, and -the bounding energies of each segment are mapped to bounding indices on each of -the nuclide energy grids. By default, OpenMC uses 8000 equal-lethargy segments -as recommended by Brown. +To speed up energy grid searches, OpenMC uses a `logarithmic mapping technique`_ +to limit the range of energies that must be searched for each nuclide. The +entire energy range is divided up into equal-lethargy segments, and the bounding +energies of each segment are mapped to bounding indices on each of the nuclide +energy grids. By default, OpenMC uses 8000 equal-lethargy segments as +recommended by Brown. Other Methods +++++++++++++ @@ -74,9 +74,9 @@ offers support for an experimental data format called windowed multipole (WMP). This data format requires less memory than pointwise cross sections, and it allows on-the-fly Doppler broadening to arbitrary temperature. -The multipole method was introduced by [Hwang]_ and the faster windowed -multipole method by [Josey]_. In the multipole format, cross section resonances -are represented by poles, :math:`p_j`, and residues, :math:`r_j`, in the complex +The multipole method was introduced by Hwang_ and the faster windowed multipole +method by Josey_. In the multipole format, cross section resonances are +represented by poles, :math:`p_j`, and residues, :math:`r_j`, in the complex plane. The 0K cross sections in the resolved resonance region can be computed by summing up a contribution from each pole: @@ -232,21 +232,10 @@ sections. This allows flexibility for the model to use highly anisotropic scattering information in the water while the fuel can be simulated with linear or even isotropic scattering. -.. only:: html - - .. rubric:: References - -.. [Brown] Forrest B. Brown, "New Hash-based Energy Lookup Algorithm for Monte - Carlo codes," LA-UR-14-24530, Los Alamos National Laboratory (2014). - -.. [Hwang] R. N. Hwang, "A Rigorous Pole Representation of Multilevel Cross - Sections and Its Practical Application," *Nucl. Sci. Eng.*, **96**, - 192-209 (1987). - -.. [Josey] Colin Josey, Pablo Ducru, Benoit Forget, and Kord Smith, "Windowed - Multipole for Cross Section Doppler Broadening," *J. Comp. Phys*, - **307**, 715-727 (2016). http://dx.doi.org/10.1016/j.jcp.2015.08.013 - +.. _logarithmic mapping technique: + https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-14-24530.pdf +.. _Hwang: http://www.ans.org/pubs/journals/nse/a_16381 +.. _Josey: http://dx.doi.org/10.1016/j.jcp.2015.08.013 .. _MCNP: http://mcnp.lanl.gov .. _Serpent: http://montecarlo.vtt.fi .. _NJOY: http://t2.lanl.gov/codes.shtml diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 57663594d..41623b6a1 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -281,6 +281,8 @@ based on the recommended value in LA-UR-14-24530_. .. note:: This element is not used in the multi-group :ref:`energy_mode`. +.. _multipole_library: + ```` Element ------------------------------- @@ -290,8 +292,8 @@ OpenMC can use it for on-the-fly Doppler-broadening of resolved resonance range cross sections. If this element is absent from the settings.xml file, the :envvar:`OPENMC_MULTIPOLE_LIBRARY` environment variable will be used. - .. note:: The element must also be set to "true" - for windowed multipole functionality. + .. note:: The :ref:`temperature_method` must also be set to "multipole" for + windowed multipole functionality. ```` Element --------------------------- @@ -707,6 +709,47 @@ survival biasing, otherwise known as implicit capture or absorption. *Default*: false +.. _temperature_default: + +```` Element +--------------------------------- + +The ```` element specifies a default temperature in Kelvin +that is to be applied to cells in the absence of an explicit cell temperature or +a material default temperature. + + *Default*: 293.6 K + +.. _temperature_method: + +```` Element +-------------------------------- + +The ```` element has an accepted value of "nearest" or +"interpolation". A value of "nearest" indicates that for each cell, the nearest +temperature at which cross sections are given is to be applied, within a given +tolerance (see :ref:`temperature_tolerance`). A value of "interpolation" +indicates that cross sections are to be interpolated between temperatures at +which nuclear data are present. A value of "multipole" indicates that the +windowed multipole method should be used to evaluate temperature-dependent cross +sections in the resolved resonance range (a :ref:`windowed multipole library +` must also be available). + + *Default*: "nearest" + +.. _temperature_tolerance: + +```` Element +----------------------------------- + +The ```` element specifies a tolerance in Kelvin that is +to be applied when the "nearest" temperature method is used. For example, if a +cell temperature is 340 K and the tolerance is 15 K, then the closest +temperature in the range of 325 K to 355 K will be used to evaluate cross +sections. + + *Default*: 10 K + ```` Element --------------------- @@ -1083,7 +1126,9 @@ Each ```` element can have the following attributes or sub-elements: specified for the "distributed temperature" feature. This will give each unique instance of the cell its own temperature. - *Default*: The temperature of the coldest nuclide in the cell's material(s) + *Default*: If a material default temperature is supplied, it is used. In the + absence of a material default temperature, the :ref:`global default + temperature ` is used. :rotation: If the cell is filled with a universe, this element specifies the angles in @@ -1289,11 +1334,12 @@ Each ``material`` element can have the following attributes or sub-elements: *Default*: "" :temperature: + An element with no attributes which is used to set the default temperature + of the material in Kelvin. - An element with no attributes which is used to set the temperature of the - material. This element accepts a maximum 6-character string that indicates - the default temperature rounded to the nearest integer in units of Kelvin, - e.g. "294K". + *Default*: If a material default temperature is not given and a cell + temperature is not specified, the :ref:`global default temperature + ` is used. :density: An element with attributes/sub-elements called ``value`` and ``units``. The @@ -1409,19 +1455,6 @@ Each ``material`` element can have the following attributes or sub-elements: .. _IUPAC Isotopic Compositions of the Elements 2009: http://pac.iupac.org/publications/pac/pdf/2011/pdf/8302x0397.pdf -```` Element ------------------------- - -In some circumstances, the temperature may be the same for many or -all materials in a given problem. In this case, rather than specifying the -```` element on every material, a ```` -element can be used to set the default material temperature for any material -without an explicitly provided temperature. This element has no attributes and -accepts a maximum 6-character string that indicates the default temperature -rounded to the nearest integer in units of Kelvin, e.g. "294K". - - *Default*: None - ------------------------------------ Tallies Specification -- tallies.xml ------------------------------------ diff --git a/examples/python/basic/build-xml.py b/examples/python/basic/build-xml.py index 802debf8c..022737fdb 100644 --- a/examples/python/basic/build-xml.py +++ b/examples/python/basic/build-xml.py @@ -33,7 +33,6 @@ fuel.add_nuclide(u235, 1.) # Instantiate a Materials collection and export to XML materials_file = openmc.Materials([moderator, fuel]) -materials_file.default_temperature = '294K' materials_file.export_to_xml() diff --git a/examples/python/boxes/build-xml.py b/examples/python/boxes/build-xml.py index 175d3f511..308019e7d 100644 --- a/examples/python/boxes/build-xml.py +++ b/examples/python/boxes/build-xml.py @@ -38,7 +38,6 @@ moderator.add_s_alpha_beta('c_H_in_H2O') # Instantiate a Materials collection and export to XML materials_file = openmc.Materials([fuel1, fuel2, moderator]) -materials_file.default_temperature = '294K' materials_file.export_to_xml() diff --git a/examples/python/lattice/hexagonal/build-xml.py b/examples/python/lattice/hexagonal/build-xml.py index b131e9dea..cca072ba7 100644 --- a/examples/python/lattice/hexagonal/build-xml.py +++ b/examples/python/lattice/hexagonal/build-xml.py @@ -37,7 +37,6 @@ iron.add_nuclide(fe56, 1.) # Instantiate a Materials collection and export to XML materials_file = openmc.Materials([moderator, fuel, iron]) -materials_file.default_temperature = '294K' materials_file.export_to_xml() diff --git a/examples/python/lattice/nested/build-xml.py b/examples/python/lattice/nested/build-xml.py index 59962c2c3..3189641e7 100644 --- a/examples/python/lattice/nested/build-xml.py +++ b/examples/python/lattice/nested/build-xml.py @@ -32,7 +32,6 @@ moderator.add_s_alpha_beta('c_H_in_H2O') # Instantiate a Materials collection and export to XML materials_file = openmc.Materials((moderator, fuel)) -materials_file.default_temperature = '294K' materials_file.export_to_xml() diff --git a/examples/python/lattice/simple/build-xml.py b/examples/python/lattice/simple/build-xml.py index ad1c465d4..2f1f8e76e 100644 --- a/examples/python/lattice/simple/build-xml.py +++ b/examples/python/lattice/simple/build-xml.py @@ -32,7 +32,6 @@ moderator.add_s_alpha_beta('c_H_in_H2O') # Instantiate a Materials collection and export to XML materials_file = openmc.Materials([moderator, fuel]) -materials_file.default_temperature = '294K' materials_file.export_to_xml() diff --git a/examples/python/pincell/build-xml.py b/examples/python/pincell/build-xml.py index a11ddb7cd..fc91ae693 100644 --- a/examples/python/pincell/build-xml.py +++ b/examples/python/pincell/build-xml.py @@ -102,7 +102,6 @@ borated_water.add_s_alpha_beta('c_H_in_H2O') # Instantiate a Materials collection and export to XML materials_file = openmc.Materials([uo2, helium, zircaloy, borated_water]) -materials_file.default_temperature = '294K' materials_file.export_to_xml() diff --git a/examples/python/pincell_multigroup/build-xml.py b/examples/python/pincell_multigroup/build-xml.py index 6dbfa336b..2b08e8274 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.300K', groups) +uo2_xsdata = openmc.XSdata('UO2', groups) uo2_xsdata.order = 0 uo2_xsdata.total = [0.1779492, 0.3298048, 0.4803882, 0.5543674, 0.3118013, 0.3951678, 0.5644058] @@ -41,7 +41,7 @@ uo2_xsdata.nu_fission = [2.005998E-02, 2.027303E-03, 1.570599E-02, uo2_xsdata.chi = [5.8791E-01, 4.1176E-01, 3.3906E-04, 1.1761E-07, 0.0000E+00, 0.0000E+00, 0.0000E+00] -h2o_xsdata = openmc.XSdata('LWTR.300K', groups) +h2o_xsdata = openmc.XSdata('LWTR', groups) h2o_xsdata.order = 0 h2o_xsdata.total = [0.15920605, 0.412969593, 0.59030986, 0.58435, 0.718, 1.2544497, 2.650379] @@ -66,8 +66,8 @@ mg_cross_sections_file.export_to_xml() ############################################################################### # Instantiate some Macroscopic Data -uo2_data = openmc.Macroscopic('UO2', '300K') -h2o_data = openmc.Macroscopic('LWTR', '300K') +uo2_data = openmc.Macroscopic('UO2') +h2o_data = openmc.Macroscopic('LWTR') # Instantiate some Materials and register the appropriate Macroscopic objects uo2 = openmc.Material(material_id=1, name='UO2 fuel') @@ -80,7 +80,6 @@ water.add_macroscopic(h2o_data) # Instantiate a Materials collection and export to XML materials_file = openmc.Materials([uo2, water]) -materials_file.default_xs = '300K' materials_file.export_to_xml() diff --git a/examples/python/reflective/build-xml.py b/examples/python/reflective/build-xml.py index 098fa8620..af86e446a 100644 --- a/examples/python/reflective/build-xml.py +++ b/examples/python/reflective/build-xml.py @@ -25,7 +25,6 @@ fuel.add_nuclide(u235, 1.) # Instantiate a Materials collection and export to XML materials_file = openmc.Materials([fuel]) -materials_file.default_temperature = '294K' materials_file.export_to_xml() diff --git a/examples/xml/basic/materials.xml b/examples/xml/basic/materials.xml index b7bc2e4e5..606c676df 100644 --- a/examples/xml/basic/materials.xml +++ b/examples/xml/basic/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/examples/xml/boxes/materials.xml b/examples/xml/boxes/materials.xml index 417fb83ba..1d0ab4a1c 100644 --- a/examples/xml/boxes/materials.xml +++ b/examples/xml/boxes/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/examples/xml/lattice/nested/materials.xml b/examples/xml/lattice/nested/materials.xml index ddd932fcb..222272195 100644 --- a/examples/xml/lattice/nested/materials.xml +++ b/examples/xml/lattice/nested/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/examples/xml/lattice/simple/materials.xml b/examples/xml/lattice/simple/materials.xml index ddd932fcb..222272195 100644 --- a/examples/xml/lattice/simple/materials.xml +++ b/examples/xml/lattice/simple/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/examples/xml/pincell/materials.xml b/examples/xml/pincell/materials.xml index 020c0e676..9f9afa384 100644 --- a/examples/xml/pincell/materials.xml +++ b/examples/xml/pincell/materials.xml @@ -1,9 +1,6 @@ - - 294K - - 300K - diff --git a/examples/xml/pincell_multigroup/mg_cross_sections.xml b/examples/xml/pincell_multigroup/mg_cross_sections.xml index 3c671a192..af1f0072b 100644 --- a/examples/xml/pincell_multigroup/mg_cross_sections.xml +++ b/examples/xml/pincell_multigroup/mg_cross_sections.xml @@ -11,8 +11,8 @@ --> - UO2.300K - UO2.300K + UO2 + UO2 2.53E-8 0 true @@ -67,8 +67,8 @@ - MOX1.300K - MOX1.300K + MOX1 + MOX1 2.53E-8 0 true @@ -124,8 +124,8 @@ - MOX2.300K - MOX2.300K + MOX2 + MOX2 2.53E-8 0 true @@ -180,8 +180,8 @@ - MOX3.300K - MOX3.300K + MOX3 + MOX3 2.53E-8 0 true @@ -236,8 +236,8 @@ - FC.300K - FC.300K + FC + FC 2.53E-8 0 true @@ -286,8 +286,8 @@ - GT.300K - GT.300K + GT + GT 2.53E-8 0 false @@ -318,8 +318,8 @@ - LWTR.300K - LWTR.300K + LWTR + LWTR 2.53E-8 0 false @@ -351,8 +351,8 @@ - CR.300K - CR.300K + CR + CR 2.53E-8 0 false diff --git a/examples/xml/reflective/materials.xml b/examples/xml/reflective/materials.xml index 37cfcdc6e..2472a7471 100644 --- a/examples/xml/reflective/materials.xml +++ b/examples/xml/reflective/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index c629dbc2a..7568a24ba 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -100,7 +100,7 @@ class IncidentNeutron(EqualityMixin): Parameters ---------- name : str - Name of the table + Name of the nuclide using the GND naming convention atomic_number : int Number of protons in the nucleus mass_number : int @@ -133,7 +133,7 @@ class IncidentNeutron(EqualityMixin): metastable : int Metastable state of the nucleus. A value of zero indicates ground state. name : str - ZAID identifier of the table, e.g. 92235. + Name of the nuclide using the GND naming convention reactions : collections.OrderedDict Contains the cross sections, secondary angle and energy distributions, and other associated data for each reaction. The keys are the MT values @@ -534,11 +534,6 @@ class IncidentNeutron(EqualityMixin): # Assign temperature to the running list kTs = [ace.temperature] - # If mass number hasn't been specified, make an educated guess - zaid, xs = ace.name.split('.') - name, element, Z, mass_number, metastable = \ - _get_metadata(int(zaid), metastable_scheme) - data = cls(name, Z, mass_number, metastable, ace.atomic_weight_ratio, kTs) diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index 21edab686..58c527605 100644 --- a/openmc/data/thermal.py +++ b/openmc/data/thermal.py @@ -8,7 +8,7 @@ import h5py import openmc.checkvalue as cv from openmc.mixin import EqualityMixin -from .data import K_BOLTZMANN +from .data import K_BOLTZMANN, ATOMIC_SYMBOL from .ace import Table, get_table from .angle_energy import AngleEnergy from .function import Tabulated1D @@ -156,7 +156,7 @@ class ThermalScattering(EqualityMixin): Parameters ---------- name : str - ZAID identifier of the table, e.g. lwtr.10t. + Name of the material using GND convention, e.g. c_H_in_H2O atomic_weight_ratio : float Atomic mass ratio of the target nuclide. kTs : Iterable of float @@ -174,7 +174,7 @@ class ThermalScattering(EqualityMixin): Inelastic scattering cross section derived in the incoherent approximation name : str - Name of the table, e.g. lwtr.20t. + Name of the material using GND convention, e.g. c_H_in_H2O temperatures : Iterable of str List of string representations the temperatures of the target nuclide in the data set. The temperatures are strings of the temperature, @@ -182,8 +182,8 @@ class ThermalScattering(EqualityMixin): kTs : Iterable of float List of temperatures of the target nuclide in the data set. The temperatures have units of MeV. - zaids : Iterable of int - ZAID identifiers that the thermal scattering data applies to + nuclides : Iterable of str + Nuclide names that the thermal scattering data applies to """ @@ -200,7 +200,7 @@ class ThermalScattering(EqualityMixin): self.inelastic_mu_out = {} self.inelastic_dist = {} self.secondary_mode = None - self.zaids = [] + self.nuclides = [] def __repr__(self): if hasattr(self, 'name'): @@ -226,7 +226,7 @@ class ThermalScattering(EqualityMixin): # Write basic data g = f.create_group(self.name) g.attrs['atomic_weight_ratio'] = self.atomic_weight_ratio - g.attrs['zaids'] = self.zaids + g.attrs['nuclides'] = np.string_(self.nuclides) g.attrs['secondary_mode'] = np.string_(self.secondary_mode) ktg = g.create_group('kTs') for i, temperature in enumerate(self.temperatures): @@ -446,7 +446,7 @@ class ThermalScattering(EqualityMixin): temperatures = [str(int(round(kT / K_BOLTZMANN))) + "K" for kT in kTs] table = cls(name, atomic_weight_ratio, kTs) - table.zaids = group.attrs['zaids'] + table.nuclides = [nuc.decode() for nuc in group.attrs['nuclides']] table.secondary_mode = group.attrs['secondary_mode'].decode() # Read thermal elastic scattering @@ -628,8 +628,10 @@ class ThermalScattering(EqualityMixin): table.elastic_mu_out[temperatures[0]].shape = \ (n_energy, n_mu) - # Get relevant ZAIDs - pairs = np.fromiter(map(lambda p: p[0], ace.pairs), int) - table.zaids = pairs[np.nonzero(pairs)] + # Get relevant nuclides + for zaid, awr in ace.pairs: + if zaid > 0: + Z, A = divmod(zaid, 1000) + table.nuclides.append(ATOMIC_SYMBOL[Z] + str(A)) return table diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py index a3a2187b7..37ad6c1be 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -99,9 +99,6 @@ class XSdata(object): Unique identifier for the xsdata object alias : str Separate unique identifier for the xsdata object - zaid : int - 1000*(atomic number) + mass number. As an example, the zaid of U235 - would be 92235. awr : float Atomic weight ratio of an isotope. That is, the ratio of the mass of the isotope to the mass of a single neutron. @@ -227,7 +224,6 @@ class XSdata(object): self._energy_groups = energy_groups self._representation = representation self._alias = None - self._zaid = None self._awr = None self._kT = None self._fissionable = False @@ -262,10 +258,6 @@ class XSdata(object): def alias(self): return self._alias - @property - def zaid(self): - return self._zaid - @property def awr(self): return self._awr @@ -396,13 +388,6 @@ class XSdata(object): else: self._alias = self._name - @zaid.setter - def zaid(self, zaid): - # Check type and value - check_type('zaid', zaid, Integral) - check_greater_than('zaid', zaid, 0) - self._zaid = zaid - @awr.setter def awr(self, awr): # Check validity of type and that the awr value is > 0 @@ -1013,18 +998,10 @@ class XSdata(object): subelement = ET.SubElement(element, 'kT') subelement.text = str(self._kT) - if self._zaid is not None: - subelement = ET.SubElement(element, 'zaid') - subelement.text = str(self._zaid) - if self._awr is not None: subelement = ET.SubElement(element, 'awr') subelement.text = str(self._awr) - if self._kT is not None: - subelement = ET.SubElement(element, 'kT') - subelement.text = str(self._kT) - if self._fissionable is not None: subelement = ET.SubElement(element, 'fissionable') subelement.text = str(self._fissionable) diff --git a/openmc/nuclide.py b/openmc/nuclide.py index 186292a87..11c59e687 100644 --- a/openmc/nuclide.py +++ b/openmc/nuclide.py @@ -20,9 +20,6 @@ class Nuclide(object): ---------- name : str Name of the nuclide, e.g. U235 - zaid : int - 1000*(atomic number) + mass number. As an example, the zaid of U235 - would be 92235. scattering : 'data' or 'iso-in-lab' or None The type of angular scattering distribution to use @@ -31,7 +28,6 @@ class Nuclide(object): def __init__(self, name=''): # Initialize class attributes self._name = '' - self._zaid = None self._scattering = None # Set the Material class attributes @@ -62,8 +58,6 @@ class Nuclide(object): def __repr__(self): string = 'Nuclide - {0}\n'.format(self._name) - if self.zaid is not None: - string += '{0: <16}{1}{2}\n'.format('\tZAID', '=\t', self.zaid) if self.scattering is not None: string += '{0: <16}{1}{2}\n'.format('\tscattering', '=\t', self.scattering) @@ -73,10 +67,6 @@ class Nuclide(object): def name(self): return self._name - @property - def zaid(self): - return self._zaid - @property def scattering(self): return self._scattering @@ -96,14 +86,8 @@ class Nuclide(object): '"{}" is being renamed as "{}".'.format(name, self._name) warnings.warn(msg) - @zaid.setter - def zaid(self, zaid): - check_type('zaid', zaid, Integral) - self._zaid = zaid - @scattering.setter def scattering(self, scattering): - if not scattering in ['data', 'iso-in-lab']: msg = 'Unable to set scattering for Nuclide to {0} ' \ 'which is not "data" or "iso-in-lab"'.format(scattering) diff --git a/openmc/settings.py b/openmc/settings.py index c1d4adc24..86aac9471 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -1,4 +1,4 @@ -from collections import Iterable, MutableSequence +from collections import Iterable, MutableSequence, Mapping from numbers import Real, Integral import warnings from xml.etree import ElementTree as ET @@ -78,8 +78,6 @@ class Settings(object): cross section library. If it is not set, the :envvar:`OPENMC_MULTIPOLE_LIBRARY` environment variable will be used. A multipole library is optional. - energy_grid : {'nuclide', 'logarithm', 'material-union'} - Set the method used to search energy grids. energy_mode : {'continuous-energy', 'multi-group'} Set whether the calculation should be continuous-energy or multi-group. max_order : int @@ -103,6 +101,14 @@ class Settings(object): Coordinates of the lower-left point of the Shannon entropy mesh entropy_upper_right : tuple or list Coordinates of the upper-right point of the Shannon entropy mesh + temperature : dict + Defines a default temperature and method for treating intermediate + temperatures at which nuclear data doesn't exist. Accepted keys are + 'default', 'method', and 'tolerance'. The value for 'default' should be + a float representing the default temperature in Kelvin. The value for + 'method' should be 'nearest' or 'interpolation'. If the method is + 'nearest', 'tolerance' indicates a range of temperature within which + cross sections may be used. trigger_active : bool Indicate whether tally triggers are used trigger_max_batches : int @@ -130,9 +136,6 @@ class Settings(object): Coordinates of the lower-left point of the UFS mesh ufs_upper_right : tuple or list Coordinates of the upper-right point of the UFS mesh - use_windowed_multipole : bool - Whether or not windowed multipole can be used to evaluate resolved - resonance cross sections. resonance_scattering : ResonanceScattering or iterable of ResonanceScattering The elastic scattering model to use for resonant isotopes volume_calculations : VolumeCalculation or iterable of VolumeCalculation @@ -160,7 +163,6 @@ class Settings(object): self._confidence_intervals = None self._cross_sections = None self._multipole_library = None - self._energy_grid = None self._ptables = None self._run_cmfd = None self._seed = None @@ -197,6 +199,8 @@ class Settings(object): self._trace = None self._track = None + self._temperature = {} + # Cutoff subelement self._weight = None self._weight_avg = None @@ -216,7 +220,6 @@ class Settings(object): self._settings_file = ET.Element("settings") self._run_mode_subelement = None - self._multipole_active = None self._resonance_scattering = cv.CheckedList( ResonanceScattering, 'resonance scattering models') @@ -271,10 +274,6 @@ class Settings(object): def multipole_library(self): return self._multipole_library - @property - def energy_grid(self): - return self._energy_grid - @property def ptables(self): return self._ptables @@ -363,6 +362,10 @@ class Settings(object): def verbosity(self): return self._verbosity + @property + def temperature(self): + return self._temperature + @property def trace(self): return self._trace @@ -415,10 +418,6 @@ class Settings(object): def dd_count_interactions(self): return self._dd_count_interactions - @property - def use_windowed_multipole(self): - return self._multipole_active - @property def resonance_scattering(self): return self._resonance_scattering @@ -593,12 +592,6 @@ class Settings(object): cv.check_type('cross sections', multipole_library, basestring) self._multipole_library = multipole_library - @energy_grid.setter - def energy_grid(self, energy_grid): - cv.check_value('energy grid', energy_grid, - ['nuclide', 'logarithm', 'material-union']) - self._energy_grid = energy_grid - @ptables.setter def ptables(self, ptables): cv.check_type('probability tables', ptables, bool) @@ -674,6 +667,21 @@ class Settings(object): cv.check_type('no reduction option', no_reduce, bool) self._no_reduce = no_reduce + @temperature.setter + def temperature(self, temperature): + cv.check_type('temperature settings', temperature, Mapping) + for key, value in temperature.items(): + cv.check_value('temperature key', key, + ['default', 'method', 'tolerance']) + if key == 'default': + cv.check_type('default temperature', value, Real) + elif key == 'method': + cv.check_value('temperature method', value, + ['nearest', 'interpolation', 'multipole']) + elif key == 'tolerance': + cv.check_type('temperature tolerance', value, Real) + self._temperature = temperature + @threads.setter def threads(self, threads): cv.check_type('number of threads', threads, Integral) @@ -801,11 +809,6 @@ class Settings(object): self._dd_count_interactions = interactions - @use_windowed_multipole.setter - def use_windowed_multipole(self, active): - cv.check_type('use_windowed_multipole', active, bool) - self._multipole_active = active - @resonance_scattering.setter def resonance_scattering(self, res): if not isinstance(res, MutableSequence): @@ -963,11 +966,6 @@ class Settings(object): element = ET.SubElement(self._settings_file, "multipole_library") element.text = str(self._multipole_library) - def _create_energy_grid_subelement(self): - if self._energy_grid is not None: - element = ET.SubElement(self._settings_file, "energy_grid") - element.text = str(self._energy_grid) - def _create_ptables_subelement(self): if self._ptables is not None: element = ET.SubElement(self._settings_file, "ptables") @@ -1050,6 +1048,13 @@ class Settings(object): element = ET.SubElement(self._settings_file, "no_reduce") element.text = str(self._no_reduce).lower() + def _create_temperature_subelements(self): + if self.temperature: + for key, value in self.temperature.items(): + element = ET.SubElement(self._settings_file, + "temperature_{}".format(key)) + element.text = str(value) + def _create_threads_subelement(self): if self._threads is not None: element = ET.SubElement(self._settings_file, "threads") @@ -1107,20 +1112,10 @@ class Settings(object): subelement = ET.SubElement(element, "count_interactions") subelement.text = str(self._dd_count_interactions).lower() - def _create_use_multipole_subelement(self): - if self._multipole_active is not None: - element = ET.SubElement(self._settings_file, - "use_windowed_multipole") - element.text = str(self._multipole_active) - def _create_resonance_scattering_subelement(self): if len(self.resonance_scattering) > 0: elem = ET.SubElement(self._settings_file, 'resonance_scattering') for r in self.resonance_scattering: - if r.nuclide.name != r.nuclide_0K.name: - raise ValueError("The nuclide and nuclide_0K attributes of " - "a ResonantScattering object must have " - "identical names.") elem.append(r.to_xml_element()) def export_to_xml(self): @@ -1142,7 +1137,6 @@ class Settings(object): self._create_confidence_intervals() self._create_cross_sections_subelement() self._create_multipole_library_subelement() - self._create_energy_grid_subelement() self._create_energy_mode_subelement() self._create_max_order_subelement() self._create_ptables_subelement() @@ -1155,11 +1149,11 @@ class Settings(object): self._create_no_reduce_subelement() self._create_threads_subelement() self._create_verbosity_subelement() + self._create_temperature_subelements() self._create_trace_subelement() self._create_track_subelement() self._create_ufs_subelement() self._create_dd_subelement() - self._create_use_multipole_subelement() self._create_resonance_scattering_subelement() self._create_volume_calcs_subelement() @@ -1175,14 +1169,26 @@ class Settings(object): class ResonanceScattering(object): """Specification of the elastic scattering model for resonant isotopes + Parameters + ---------- + nuclide : openmc.Nuclide + The nuclide affected by this resonance scattering treatment. + method : {'ARES', 'CXS', 'DBRC', 'WCM'} + The method used to sample outgoing scattering energies. Valid options + are 'ARES', 'CXS' (constant cross section), 'DBRC' (Doppler broadening + rejection correction), and 'WCM' (weight correction method). + E_min : float + The minimum energy above which the specified method is applied. By + default, CXS will be used below E_min. + E_max : float + The maximum energy below which the specified method is applied. By + default, the asymptotic target-at-rest model is applied above E_max. + Attributes ---------- nuclide : openmc.Nuclide The nuclide affected by this resonance scattering treatment. - nuclide_0K : openmc.Nuclide - This should be the same isotope as the nuclide attribute above, but it - should have an xs attribute that identifies 0 Kelvin data. - method : str + method : {'ARES', 'CXS', 'DBRC', 'WCM'} The method used to sample outgoing scattering energies. Valid options are 'ARES', 'CXS' (constant cross section), 'DBRC' (Doppler broadening rejection correction), and 'WCM' (weight correction method). @@ -1195,21 +1201,20 @@ class ResonanceScattering(object): """ - def __init__(self): - self._nuclide = None - self._nuclide_0K = None - self._method = None + def __init__(self, nuclide, method='CXS', E_min=None, E_max=None): self._E_min = None self._E_max = None + self.nuclide = nuclide + self.method = method + if E_min is not None: + self.E_min = E_min + if E_max is not None: + self.E_max = E_max @property def nuclide(self): return self._nuclide - @property - def nuclide_0K(self): - return self._nuclide_0K - @property def method(self): return self._method @@ -1227,11 +1232,6 @@ class ResonanceScattering(object): cv.check_type('nuclide', nuc, Nuclide) self._nuclide = nuc - @nuclide_0K.setter - def nuclide_0K(self, nuc): - cv.check_type('nuclide_0K', nuc, Nuclide) - self._nuclide_0K = nuc - @method.setter def method(self, m): cv.check_value('method', m, ('ARES', 'CXS', 'DBRC', 'WCM')) diff --git a/openmc/summary.py b/openmc/summary.py index 3e25f2c26..38c0e335c 100644 --- a/openmc/summary.py +++ b/openmc/summary.py @@ -83,11 +83,10 @@ class Summary(object): n_nuclides = self._f['nuclides/n_nuclides_total'].value names = self._f['nuclides/names'].value awrs = self._f['nuclides/awrs'].value - zaids = self._f['nuclides/zaids'].value for n in range(n_nuclides): name = names[n].decode() name = name[:name.find('.')] - self.nuclides[name] = (zaids[n], awrs[n]) + self.nuclides[name] = awrs[n] def _read_geometry(self): # Read in and initialize the Materials and Geometry @@ -113,7 +112,6 @@ class Summary(object): material_id = int(key.lstrip('material ')) index = self._f['materials'][key]['index'].value name = self._f['materials'][key]['name'].value.decode() - temperature = self._f['materials'][key]['temperature'].value.decode() density = self._f['materials'][key]['atom_density'].value nuc_densities = self._f['materials'][key]['nuclide_densities'][...] nuclides = self._f['materials'][key]['nuclides'].value diff --git a/src/constants.F90 b/src/constants.F90 index 6d5f2fc10..ff2bd3710 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -268,6 +268,12 @@ module constants JENDL_33 = 7, & JENDL_40 = 8 + ! Temperature treatment method + integer, parameter :: & + TEMPERATURE_NEAREST = 1, & + TEMPERATURE_INTERPOLATION = 2, & + TEMPERATURE_MULTIPOLE = 3 + ! ============================================================================ ! TALLY-RELATED CONSTANTS @@ -401,12 +407,6 @@ module constants integer, parameter :: ERROR_INT = -huge(0) real(8), parameter :: ERROR_REAL = -huge(0.0_8) * 0.917826354_8 - ! Energy grid methods - integer, parameter :: & - GRID_NUCLIDE = 1, & ! unique energy grid for each nuclide - GRID_MAT_UNION = 2, & ! material union grids with pointers - GRID_LOGARITHM = 3 ! lethargy mapping - ! Running modes integer, parameter :: & MODE_FIXEDSOURCE = 1, & ! Fixed source mode diff --git a/src/cross_section.F90 b/src/cross_section.F90 index b67226a1b..43e854496 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -2,7 +2,7 @@ module cross_section use algorithm, only: binary_search use constants - use energy_grid, only: grid_method, log_spacing + use energy_grid, only: log_spacing use error, only: fatal_error use global use list_header, only: ListElemInt @@ -37,7 +37,6 @@ contains ! union grid real(8) :: atom_density ! atom density of a nuclide logical :: check_sab ! should we check for S(a,b) table? - type(Material), pointer :: mat ! current material ! Set all material macroscopic cross sections to zero material_xs % total = ZERO @@ -49,89 +48,83 @@ contains ! Exit subroutine if material is void if (p % material == MATERIAL_VOID) return - mat => materials(p % material) - - ! Find energy index on energy grid - if (grid_method == GRID_MAT_UNION) then - i_grid = find_energy_index(mat, p % E) - else if (grid_method == GRID_LOGARITHM) then + associate (mat => materials(p % material)) + ! Find energy index on energy grid i_grid = int(log(p % E/energy_min_neutron)/log_spacing) - end if - ! Determine if this material has S(a,b) tables - check_sab = (mat % n_sab > 0) + ! Determine if this material has S(a,b) tables + check_sab = (mat % n_sab > 0) - ! Initialize position in i_sab_nuclides - j = 1 + ! Initialize position in i_sab_nuclides + j = 1 - ! Add contribution from each nuclide in material - do i = 1, mat % n_nuclides - ! ======================================================================== - ! CHECK FOR S(A,B) TABLE + ! Add contribution from each nuclide in material + do i = 1, mat % n_nuclides + ! ======================================================================== + ! CHECK FOR S(A,B) TABLE - i_sab = 0 + i_sab = 0 - ! Check if this nuclide matches one of the S(a,b) tables specified -- this - ! relies on i_sab_nuclides being in sorted order - if (check_sab) then - if (i == mat % i_sab_nuclides(j)) then - ! Get index in sab_tables - i_sab = mat % i_sab_tables(j) + ! Check if this nuclide matches one of the S(a,b) tables specified -- this + ! relies on i_sab_nuclides being in sorted order + if (check_sab) then + if (i == mat % i_sab_nuclides(j)) then + ! Get index in sab_tables + i_sab = mat % i_sab_tables(j) - ! If particle energy is greater than the highest energy for the S(a,b) - ! table, don't use the S(a,b) table - if (p % E > sab_tables(i_sab) % threshold_inelastic) i_sab = 0 + ! If particle energy is greater than the highest energy for the S(a,b) + ! table, don't use the S(a,b) table + if (p % E > sab_tables(i_sab) % data(1) % threshold_inelastic) i_sab = 0 - ! Increment position in i_sab_nuclides - j = j + 1 + ! Increment position in i_sab_nuclides + j = j + 1 - ! Don't check for S(a,b) tables if there are no more left - if (j > mat % n_sab) check_sab = .false. + ! Don't check for S(a,b) tables if there are no more left + if (j > mat % n_sab) check_sab = .false. + end if end if - end if - ! ======================================================================== - ! CALCULATE MICROSCOPIC CROSS SECTION + ! ======================================================================== + ! CALCULATE MICROSCOPIC CROSS SECTION - ! Determine microscopic cross sections for this nuclide - i_nuclide = mat % nuclide(i) + ! Determine microscopic cross sections for this nuclide + i_nuclide = mat % nuclide(i) - ! Calculate microscopic cross section for this nuclide - if (p % E /= micro_xs(i_nuclide) % last_E & - .or. p % sqrtkT /= micro_xs(i_nuclide) % last_sqrtkT) then - call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, & - i_grid, p % sqrtkT) - else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then - call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, & - i_grid, p % sqrtkT) - end if + ! Calculate microscopic cross section for this nuclide + if (p % E /= micro_xs(i_nuclide) % last_E & + .or. p % sqrtkT /= micro_xs(i_nuclide) % last_sqrtkT) then + call calculate_nuclide_xs(i_nuclide, i_sab, p % E, i_grid, p % sqrtkT) + else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then + call calculate_nuclide_xs(i_nuclide, i_sab, p % E, i_grid, p % sqrtkT) + end if - ! ======================================================================== - ! ADD TO MACROSCOPIC CROSS SECTION + ! ======================================================================== + ! ADD TO MACROSCOPIC CROSS SECTION - ! Copy atom density of nuclide in material - atom_density = mat % atom_density(i) + ! Copy atom density of nuclide in material + atom_density = mat % atom_density(i) - ! Add contributions to material macroscopic total cross section - material_xs % total = material_xs % total + & - atom_density * micro_xs(i_nuclide) % total + ! Add contributions to material macroscopic total cross section + material_xs % total = material_xs % total + & + atom_density * micro_xs(i_nuclide) % total - ! Add contributions to material macroscopic scattering cross section - material_xs % elastic = material_xs % elastic + & - atom_density * micro_xs(i_nuclide) % elastic + ! Add contributions to material macroscopic scattering cross section + material_xs % elastic = material_xs % elastic + & + atom_density * micro_xs(i_nuclide) % elastic - ! Add contributions to material macroscopic absorption cross section - material_xs % absorption = material_xs % absorption + & - atom_density * micro_xs(i_nuclide) % absorption + ! Add contributions to material macroscopic absorption cross section + material_xs % absorption = material_xs % absorption + & + atom_density * micro_xs(i_nuclide) % absorption - ! Add contributions to material macroscopic fission cross section - material_xs % fission = material_xs % fission + & - atom_density * micro_xs(i_nuclide) % fission + ! Add contributions to material macroscopic fission cross section + material_xs % fission = material_xs % fission + & + atom_density * micro_xs(i_nuclide) % fission - ! Add contributions to material macroscopic nu-fission cross section - material_xs % nu_fission = material_xs % nu_fission + & - atom_density * micro_xs(i_nuclide) % nu_fission - end do + ! Add contributions to material macroscopic nu-fission cross section + material_xs % nu_fission = material_xs % nu_fission + & + atom_density * micro_xs(i_nuclide) % nu_fission + end do + end associate end subroutine calculate_xs @@ -140,169 +133,160 @@ contains ! given index in the nuclides array at the energy of the given particle !=============================================================================== - subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_mat, i_nuc_mat, & - i_log_union, sqrtkT) + subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_log_union, sqrtkT) integer, intent(in) :: i_nuclide ! index into nuclides array integer, intent(in) :: i_sab ! index into sab_tables array real(8), intent(in) :: E ! energy - integer, intent(in) :: i_mat ! index into materials array - integer, intent(in) :: i_nuc_mat ! index into nuclides array for a material integer, intent(in) :: i_log_union ! index into logarithmic mapping array or ! material union energy grid real(8), intent(in) :: sqrtkT ! Square root of kT, material dependent logical :: use_mp ! true if XS can be calculated with windowed multipole + integer :: i_temp ! index for temperature integer :: i_grid ! index on nuclide energy grid integer :: i_low ! lower logarithmic mapping index integer :: i_high ! upper logarithmic mapping index real(8) :: f ! interp factor on nuclide energy grid + real(8) :: kT ! temperature in MeV real(8) :: sigT, sigA, sigF ! Intermediate multipole variables - type(Nuclide), pointer :: nuc - type(Material), pointer :: mat - ! Set pointer to nuclide and material - nuc => nuclides(i_nuclide) - mat => materials(i_mat) - - ! Check to see if there is multipole data present at this energy - use_mp = .false. - if (nuc % mp_present) then - if (E >= nuc % multipole % start_E/1.0e6_8 .and. & - E <= nuc % multipole % end_E/1.0e6_8) then - use_mp = .true. - end if - end if - - ! Evaluate multipole or interpolate - if (use_mp) then - ! Call multipole kernel - call multipole_eval(nuc % multipole, E, sqrtkT, sigT, sigA, sigF) - - micro_xs(i_nuclide) % total = sigT - micro_xs(i_nuclide) % absorption = sigA - micro_xs(i_nuclide) % elastic = sigT - sigA - - if (nuc % fissionable) then - micro_xs(i_nuclide) % fission = sigF - micro_xs(i_nuclide) % nu_fission = sigF * nuc % nu(E, EMISSION_TOTAL) + associate (nuc => nuclides(i_nuclide)) + ! Check to see if there is multipole data present at this energy + use_mp = .false. + if (nuc % mp_present) then + if (E >= nuc % multipole % start_E/1.0e6_8 .and. & + E <= nuc % multipole % end_E/1.0e6_8) then + use_mp = .true. + else + ! If using multipole data but outside the RRR, just pick the first + ! temperature + i_temp = 1 + end if else - micro_xs(i_nuclide) % fission = ZERO - micro_xs(i_nuclide) % nu_fission = ZERO + ! If not using multipole data, do a linear search on temperature + kT = sqrtkT**2 + do i_temp = 1, size(nuclides(i_nuclide) % kTs) + if (abs(nuclides(i_nuclide) % kTs(i_temp) - kT) < & + K_BOLTZMANN*temperature_tolerance) exit + end do end if - ! Ensure these values are set - ! Note, the only time either is used is in one of 4 places: - ! 1. physics.F90 - scatter - For inelastic scatter. - ! 2. physics.F90 - sample_fission - For partial fissions. - ! 3. tally.F90 - score_general - For tallying on MTxxx reactions. - ! 4. cross_section.F90 - calculate_urr_xs - For unresolved purposes. - ! It is worth noting that none of these occur in the resolved - ! resonance range, so the value here does not matter. - micro_xs(i_nuclide) % index_grid = 0 - micro_xs(i_nuclide) % interp_factor = ZERO - else - ! Determine index on nuclide energy grid - select case (grid_method) - case (GRID_MAT_UNION) + ! Evaluate multipole or interpolate + if (use_mp) then + ! Call multipole kernel + call multipole_eval(nuc % multipole, E, sqrtkT, sigT, sigA, sigF) - i_grid = mat % nuclide_grid_index(i_nuc_mat, i_log_union) + micro_xs(i_nuclide) % total = sigT + micro_xs(i_nuclide) % absorption = sigA + micro_xs(i_nuclide) % elastic = sigT - sigA - case (GRID_LOGARITHM) - ! Determine the energy grid index using a logarithmic mapping to reduce - ! the energy range over which a binary search needs to be performed - - if (E < nuc % energy(1)) then - i_grid = 1 - elseif (E > nuc % energy(nuc % n_grid)) then - i_grid = nuc % n_grid - 1 + if (nuc % fissionable) then + micro_xs(i_nuclide) % fission = sigF + micro_xs(i_nuclide) % nu_fission = sigF * nuc % nu(E, EMISSION_TOTAL) else - ! Determine bounding indices based on which equal log-spaced interval - ! the energy is in - i_low = nuc % grid_index(i_log_union) - i_high = nuc % grid_index(i_log_union + 1) + 1 - - ! Perform binary search over reduced range - i_grid = binary_search(nuc % energy(i_low:i_high), & - i_high - i_low + 1, E) + i_low - 1 + micro_xs(i_nuclide) % fission = ZERO + micro_xs(i_nuclide) % nu_fission = ZERO end if - case (GRID_NUCLIDE) - ! Perform binary search on the nuclide energy grid in order to determine - ! which points to interpolate between + ! Ensure these values are set + ! Note, the only time either is used is in one of 4 places: + ! 1. physics.F90 - scatter - For inelastic scatter. + ! 2. physics.F90 - sample_fission - For partial fissions. + ! 3. tally.F90 - score_general - For tallying on MTxxx reactions. + ! 4. cross_section.F90 - calculate_urr_xs - For unresolved purposes. + ! It is worth noting that none of these occur in the resolved + ! resonance range, so the value here does not matter. + micro_xs(i_nuclide) % index_temp = i_temp + micro_xs(i_nuclide) % index_grid = 0 + micro_xs(i_nuclide) % interp_factor = ZERO + else + associate (grid => nuc % grid(i_temp), xs => nuc % sum_xs(i_temp)) + ! Determine the energy grid index using a logarithmic mapping to reduce + ! the energy range over which a binary search needs to be performed - if (E <= nuc % energy(1)) then - i_grid = 1 - elseif (E > nuc % energy(nuc % n_grid)) then - i_grid = nuc % n_grid - 1 - else - i_grid = binary_search(nuc % energy, nuc % n_grid, E) + if (E < grid % energy(1)) then + i_grid = 1 + elseif (E > grid % energy(size(grid % energy))) then + i_grid = size(grid % energy) - 1 + else + ! Determine bounding indices based on which equal log-spaced interval + ! the energy is in + i_low = grid % grid_index(i_log_union) + i_high = grid % grid_index(i_log_union + 1) + 1 + + ! Perform binary search over reduced range + i_grid = binary_search(grid % energy(i_low:i_high), & + i_high - i_low + 1, E) + i_low - 1 + end if + + ! check for rare case where two energy points are the same + if (grid % energy(i_grid) == grid % energy(i_grid + 1)) & + i_grid = i_grid + 1 + + ! calculate interpolation factor + f = (E - grid % energy(i_grid)) / & + (grid % energy(i_grid + 1) - grid % energy(i_grid)) + + micro_xs(i_nuclide) % index_temp = i_temp + micro_xs(i_nuclide) % index_grid = i_grid + micro_xs(i_nuclide) % interp_factor = f + + ! Initialize nuclide cross-sections to zero + micro_xs(i_nuclide) % fission = ZERO + micro_xs(i_nuclide) % nu_fission = ZERO + + ! Calculate microscopic nuclide total cross section + micro_xs(i_nuclide) % total = (ONE - f) * xs % total(i_grid) & + + f * xs % total(i_grid + 1) + + ! Calculate microscopic nuclide elastic cross section + micro_xs(i_nuclide) % elastic = (ONE - f) * xs % elastic(i_grid) & + + f * xs % elastic(i_grid + 1) + + ! Calculate microscopic nuclide absorption cross section + micro_xs(i_nuclide) % absorption = (ONE - f) * xs % absorption( & + i_grid) + f * xs % absorption(i_grid + 1) + + if (nuc % fissionable) then + ! Calculate microscopic nuclide total cross section + micro_xs(i_nuclide) % fission = (ONE - f) * xs % fission(i_grid) & + + f * xs % fission(i_grid + 1) + + ! Calculate microscopic nuclide nu-fission cross section + micro_xs(i_nuclide) % nu_fission = (ONE - f) * xs % nu_fission( & + i_grid) + f * xs % nu_fission(i_grid + 1) + end if + end associate + end if + + ! Initialize sab treatment to false + micro_xs(i_nuclide) % index_sab = NONE + micro_xs(i_nuclide) % elastic_sab = ZERO + + ! Initialize URR probability table treatment to false + micro_xs(i_nuclide) % use_ptable = .false. + + ! If there is S(a,b) data for this nuclide, we need to do a few + ! things. Since the total cross section was based on non-S(a,b) data, we + ! need to correct it by subtracting the non-S(a,b) elastic cross section and + ! then add back in the calculated S(a,b) elastic+inelastic cross section. + + if (i_sab > 0) call calculate_sab_xs(i_nuclide, i_sab, E, sqrtkT) + + ! if the particle is in the unresolved resonance range and there are + ! probability tables, we need to determine cross sections from the table + + if (urr_ptables_on .and. nuc % urr_present .and. .not. use_mp) then + if (E > nuc % urr_data(i_temp) % energy(1) .and. E < nuc % & + urr_data(i_temp) % energy(nuc % urr_data(i_temp) % n_energy)) then + call calculate_urr_xs(i_nuclide, i_temp, E) end if - - end select - - ! check for rare case where two energy points are the same - if (nuc % energy(i_grid) == nuc % energy(i_grid+1)) i_grid = i_grid + 1 - - ! calculate interpolation factor - f = (E - nuc%energy(i_grid))/(nuc%energy(i_grid+1) - nuc%energy(i_grid)) - - micro_xs(i_nuclide) % index_grid = i_grid - micro_xs(i_nuclide) % interp_factor = f - - ! Initialize nuclide cross-sections to zero - micro_xs(i_nuclide) % fission = ZERO - micro_xs(i_nuclide) % nu_fission = ZERO - - ! Calculate microscopic nuclide total cross section - micro_xs(i_nuclide) % total = (ONE - f) * nuc % total(i_grid) & - + f * nuc % total(i_grid+1) - - ! Calculate microscopic nuclide elastic cross section - micro_xs(i_nuclide) % elastic = (ONE - f) * nuc % elastic(i_grid) & - + f * nuc % elastic(i_grid+1) - - ! Calculate microscopic nuclide absorption cross section - micro_xs(i_nuclide) % absorption = (ONE - f) * nuc % absorption( & - i_grid) + f * nuc % absorption(i_grid+1) - - if (nuc % fissionable) then - ! Calculate microscopic nuclide total cross section - micro_xs(i_nuclide) % fission = (ONE - f) * nuc % fission(i_grid) & - + f * nuc % fission(i_grid+1) - - ! Calculate microscopic nuclide nu-fission cross section - micro_xs(i_nuclide) % nu_fission = (ONE - f) * nuc % nu_fission( & - i_grid) + f * nuc % nu_fission(i_grid+1) end if - end if - ! Initialize sab treatment to false - micro_xs(i_nuclide) % index_sab = NONE - micro_xs(i_nuclide) % elastic_sab = ZERO - - ! Initialize URR probability table treatment to false - micro_xs(i_nuclide) % use_ptable = .false. - - ! If there is S(a,b) data for this nuclide, we need to do a few - ! things. Since the total cross section was based on non-S(a,b) data, we - ! need to correct it by subtracting the non-S(a,b) elastic cross section and - ! then add back in the calculated S(a,b) elastic+inelastic cross section. - - if (i_sab > 0) call calculate_sab_xs(i_nuclide, i_sab, E) - - ! if the particle is in the unresolved resonance range and there are - ! probability tables, we need to determine cross sections from the table - - if (urr_ptables_on .and. nuc % urr_present) then - if (E > nuc % urr_data % energy(1) .and. & - E < nuc % urr_data % energy(nuc % urr_data % n_energy)) then - call calculate_urr_xs(i_nuclide, E) - end if - end if - - micro_xs(i_nuclide) % last_E = E - micro_xs(i_nuclide) % last_index_sab = i_sab - micro_xs(i_nuclide) % last_sqrtkT = sqrtkT + micro_xs(i_nuclide) % last_E = E + micro_xs(i_nuclide) % last_index_sab = i_sab + micro_xs(i_nuclide) % last_sqrtkT = sqrtkT + end associate end subroutine calculate_nuclide_xs @@ -312,75 +296,85 @@ contains ! whatever data were taken from the normal Nuclide table. !=============================================================================== - subroutine calculate_sab_xs(i_nuclide, i_sab, E) + subroutine calculate_sab_xs(i_nuclide, i_sab, E, sqrtkT) integer, intent(in) :: i_nuclide ! index into nuclides array integer, intent(in) :: i_sab ! index into sab_tables array real(8), intent(in) :: E ! energy + real(8), intent(in) :: sqrtkT ! temperature integer :: i_grid ! index on S(a,b) energy grid + integer :: i_temp ! temperature index real(8) :: f ! interp factor on S(a,b) energy grid real(8) :: inelastic ! S(a,b) inelastic cross section real(8) :: elastic ! S(a,b) elastic cross section - type(SAlphaBeta), pointer :: sab + real(8) :: kT ! Set flag that S(a,b) treatment should be used for scattering micro_xs(i_nuclide) % index_sab = i_sab + ! Determine temperature for S(a,b) table + kT = sqrtkT**2 + do i_temp = 1, size(sab_tables(i_sab) % kTs) + if (abs(sab_tables(i_sab) % kTs(i_temp) - kT) < & + K_BOLTZMANN*temperature_tolerance) exit + end do + ! Get pointer to S(a,b) table - sab => sab_tables(i_sab) + associate (sab => sab_tables(i_sab) % data(i_temp)) - ! Get index and interpolation factor for inelastic grid - if (E < sab % inelastic_e_in(1)) then - i_grid = 1 - f = ZERO - else - i_grid = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, E) - f = (E - sab%inelastic_e_in(i_grid)) / & - (sab%inelastic_e_in(i_grid+1) - sab%inelastic_e_in(i_grid)) - end if + ! Get index and interpolation factor for inelastic grid + if (E < sab % inelastic_e_in(1)) then + i_grid = 1 + f = ZERO + else + i_grid = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, E) + f = (E - sab%inelastic_e_in(i_grid)) / & + (sab%inelastic_e_in(i_grid+1) - sab%inelastic_e_in(i_grid)) + end if - ! Calculate S(a,b) inelastic scattering cross section - inelastic = (ONE - f) * sab % inelastic_sigma(i_grid) + & - f * sab % inelastic_sigma(i_grid + 1) + ! Calculate S(a,b) inelastic scattering cross section + inelastic = (ONE - f) * sab % inelastic_sigma(i_grid) + & + f * sab % inelastic_sigma(i_grid + 1) - ! Check for elastic data - if (E < sab % threshold_elastic) then - ! Determine whether elastic scattering is given in the coherent or - ! incoherent approximation. For coherent, the cross section is - ! represented as P/E whereas for incoherent, it is simply P + ! Check for elastic data + if (E < sab % threshold_elastic) then + ! Determine whether elastic scattering is given in the coherent or + ! incoherent approximation. For coherent, the cross section is + ! represented as P/E whereas for incoherent, it is simply P - if (sab % elastic_mode == SAB_ELASTIC_EXACT) then - if (E < sab % elastic_e_in(1)) then - ! If energy is below that of the lowest Bragg peak, the elastic - ! cross section will be zero - elastic = ZERO + if (sab % elastic_mode == SAB_ELASTIC_EXACT) then + if (E < sab % elastic_e_in(1)) then + ! If energy is below that of the lowest Bragg peak, the elastic + ! cross section will be zero + elastic = ZERO + else + i_grid = binary_search(sab % elastic_e_in, & + sab % n_elastic_e_in, E) + elastic = sab % elastic_P(i_grid) / E + end if else - i_grid = binary_search(sab % elastic_e_in, & - sab % n_elastic_e_in, E) - elastic = sab % elastic_P(i_grid) / E + ! Determine index on elastic energy grid + if (E < sab % elastic_e_in(1)) then + i_grid = 1 + else + i_grid = binary_search(sab % elastic_e_in, & + sab % n_elastic_e_in, E) + end if + + ! Get interpolation factor for elastic grid + f = (E - sab%elastic_e_in(i_grid))/(sab%elastic_e_in(i_grid+1) - & + sab%elastic_e_in(i_grid)) + + ! Calculate S(a,b) elastic scattering cross section + elastic = (ONE - f) * sab % elastic_P(i_grid) + & + f * sab % elastic_P(i_grid + 1) end if else - ! Determine index on elastic energy grid - if (E < sab % elastic_e_in(1)) then - i_grid = 1 - else - i_grid = binary_search(sab % elastic_e_in, & - sab % n_elastic_e_in, E) - end if - - ! Get interpolation factor for elastic grid - f = (E - sab%elastic_e_in(i_grid))/(sab%elastic_e_in(i_grid+1) - & - sab%elastic_e_in(i_grid)) - - ! Calculate S(a,b) elastic scattering cross section - elastic = (ONE - f) * sab % elastic_P(i_grid) + & - f * sab % elastic_P(i_grid + 1) + ! No elastic data + elastic = ZERO end if - else - ! No elastic data - elastic = ZERO - end if + end associate ! Correct total and elastic cross sections micro_xs(i_nuclide) % total = micro_xs(i_nuclide) % total - & @@ -390,6 +384,9 @@ contains ! Store S(a,b) elastic cross section for sampling later micro_xs(i_nuclide) % elastic_sab = elastic + ! Save temperature index + micro_xs(i_nuclide) % index_temp_sab = i_temp + end subroutine calculate_sab_xs !=============================================================================== @@ -397,9 +394,9 @@ contains ! from probability tables !=============================================================================== - subroutine calculate_urr_xs(i_nuclide, E) - + subroutine calculate_urr_xs(i_nuclide, i_temp, E) integer, intent(in) :: i_nuclide ! index into nuclides array + integer, intent(in) :: i_temp ! temperature index real(8), intent(in) :: E ! energy integer :: i_energy ! index for energy @@ -414,7 +411,7 @@ contains micro_xs(i_nuclide) % use_ptable = .true. - associate (nuc => nuclides(i_nuclide), urr => nuclides(i_nuclide) % urr_data) + associate (nuc => nuclides(i_nuclide), urr => nuclides(i_nuclide) % urr_data(i_temp)) ! determine energy table i_energy = 1 do @@ -433,7 +430,7 @@ contains ! random number for the same nuclide at different temperatures, therefore ! preserving correlation of temperature in probability tables. call prn_set_stream(STREAM_URR_PTABLE) - r = future_prn(int(nuc_zaid_dict % get_key(nuc % zaid), 8)) + r = future_prn(int(i_nuclide, 8)) call prn_set_stream(STREAM_TRACKING) i_low = 1 @@ -497,10 +494,10 @@ contains f = micro_xs(i_nuclide) % interp_factor ! Determine inelastic scattering cross section - associate (rxn => nuc % reactions(nuc % urr_inelastic)) - if (i_energy >= rxn % threshold) then - inelastic = (ONE - f) * rxn % sigma(i_energy - rxn%threshold + 1) + & - f * rxn % sigma(i_energy - rxn%threshold + 2) + associate (xs => nuc % reactions(nuc % urr_inelastic) % xs(i_temp)) + if (i_energy >= xs % threshold) then + inelastic = (ONE - f) * xs % value(i_energy - xs % threshold + 1) + & + f * xs % value(i_energy - xs % threshold + 2) end if end associate end if diff --git a/src/energy_grid.F90 b/src/energy_grid.F90 index 66419f83c..47408943e 100644 --- a/src/energy_grid.F90 +++ b/src/energy_grid.F90 @@ -13,64 +13,18 @@ module energy_grid contains -!=============================================================================== -! UNIONIZED_GRID creates a unionized energy grid, for the entire problem or for -! each material, composed of the grids from each nuclide in the entire problem, -! or each material, respectively. Right now, the grid for each nuclide is added -! into a linked list one at a time with an effective insertion sort. Could be -! done with a hash for all energy points and then a quicksort at the end (what -! hash function to use?) -!=============================================================================== - - subroutine unionized_grid() - - integer :: i ! index in nuclides array - integer :: j ! index in materials array - type(ListReal) :: list - type(Nuclide), pointer :: nuc - type(Material), pointer :: mat - - call write_message("Creating unionized energy grid...", 5) - - ! add grid points for each nuclide in the material - do j = 1, n_materials - mat => materials(j) - do i = 1, mat % n_nuclides - nuc => nuclides(mat % nuclide(i)) - call add_grid_points(list, nuc % energy) - end do - - ! set size of unionized material energy grid - mat % n_grid = list % size() - - ! create allocated array from linked list - allocate(mat % e_grid(mat % n_grid)) - do i = 1, mat % n_grid - mat % e_grid(i) = list % get_item(i) - end do - - ! delete linked list and dictionary - call list % clear() - end do - - ! Set pointers to unionized energy grid for each nuclide - call grid_pointers() - - end subroutine unionized_grid - !=============================================================================== ! LOGARITHMIC_GRID determines a logarithmic mapping for energies to bounding ! indices on a nuclide energy grid !=============================================================================== subroutine logarithmic_grid() - integer :: i, j, k ! Loop indices + integer :: t ! temperature index integer :: M ! Number of equally log-spaced bins real(8) :: E_max ! Maximum energy in MeV real(8) :: E_min ! Minimum energy in MeV real(8), allocatable :: umesh(:) ! Equally log-spaced energy grid - type(Nuclide), pointer :: nuc ! Set minimum/maximum energies E_max = energy_max_neutron @@ -85,123 +39,29 @@ contains umesh(:) = [(i*log_spacing, i=0, M)] do i = 1, n_nuclides_total - ! Allocate logarithmic mapping for nuclide - nuc => nuclides(i) - allocate(nuc % grid_index(0:M)) + associate (nuc => nuclides(i)) + do t = 1, size(nuc % grid) + ! Allocate logarithmic mapping for nuclide + allocate(nuc % grid(t) % grid_index(0:M)) - ! Determine corresponding indices in nuclide grid to energies on - ! equal-logarithmic grid - j = 1 - do k = 0, M - do while (log(nuc%energy(j + 1)/E_min) <= umesh(k)) - ! Ensure that for isotopes where maxval(nuc % energy) << E_max - ! that there are no out-of-bounds issues. - if (j + 1 == nuc % n_grid) then - exit - end if - j = j + 1 + ! Determine corresponding indices in nuclide grid to energies on + ! equal-logarithmic grid + j = 1 + do k = 0, M + do while (log(nuc % grid(t) % energy(j + 1)/E_min) <= umesh(k)) + ! Ensure that for isotopes where maxval(nuc % energy) << E_max + ! that there are no out-of-bounds issues. + if (j + 1 == size(nuc % grid(t) % energy)) exit + j = j + 1 + end do + nuc % grid(t) % grid_index(k) = j + end do end do - nuc % grid_index(k) = j - end do + end associate end do deallocate(umesh) end subroutine logarithmic_grid -!=============================================================================== -! ADD_GRID_POINTS adds energy points from the 'energy' array into a linked list -! of points already stored from previous arrays. -!=============================================================================== - - subroutine add_grid_points(list, energy) - - type(ListReal) :: list - real(8), intent(in) :: energy(:) - - integer :: i ! index in energy array - integer :: n ! size of energy array - integer :: current ! current index - real(8) :: E ! actual energy value - - i = 1 - n = size(energy) - - ! Set current index to beginning of the list - current = 1 - - do while (i <= n) - E = energy(i) - - ! If we've reached the end of the grid energy list, add the remaining - ! energy points to the end - if (current > list % size()) then - ! Finish remaining energies - do while (i <= n) - call list % append(energy(i)) - i = i + 1 - end do - exit - end if - - if (E < list % get_item(current)) then - - ! Insert new energy in this position - call list % insert(current, E) - - ! Advance index in linked list and in new energy grid - i = i + 1 - current = current + 1 - - elseif (E == list % get_item(current)) then - ! Found the exact same energy, no need to store duplicates so just - ! skip and move to next index - i = i + 1 - current = current + 1 - else - current = current + 1 - end if - - end do - - end subroutine add_grid_points - -!=============================================================================== -! GRID_POINTERS creates an array of pointers (ints) for each nuclide to link -! each point on the nuclide energy grid to one on a unionized energy grid -!=============================================================================== - - subroutine grid_pointers() - - integer :: i ! loop index for nuclides - integer :: j ! loop index for nuclide energy grid - integer :: k ! loop index for materials - integer :: index_e ! index on union energy grid - real(8) :: union_energy ! energy on union grid - real(8) :: energy ! energy on nuclide grid - type(Nuclide), pointer :: nuc - type(Material), pointer :: mat - - do k = 1, n_materials - mat => materials(k) - allocate(mat % nuclide_grid_index(mat % n_nuclides, mat % n_grid)) - do i = 1, mat % n_nuclides - nuc => nuclides(mat % nuclide(i)) - - index_e = 1 - energy = nuc % energy(index_e) - - do j = 1, mat % n_grid - union_energy = mat % e_grid(j) - if (union_energy >= energy .and. index_e < nuc % n_grid) then - index_e = index_e + 1 - energy = nuc % energy(index_e) - end if - mat % nuclide_grid_index(i,j) = index_e - 1 - end do - end do - end do - - end subroutine grid_pointers - end module energy_grid diff --git a/src/geometry.F90 b/src/geometry.F90 index 6d4ca7767..d41ca1447 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -248,14 +248,15 @@ contains ! ====================================================================== ! AT LOWEST UNIVERSE, TERMINATE SEARCH - ! Set the particle material + ! Save previous material and temperature p % last_material = p % material - if (size(c % material) == 1) then - ! Only one material for this cell; assign that one to the particle. - p % material = c % material(1) - else - ! Distributed instances of this cell have different materials. - ! Determine which instance this is and assign the matching material. + p % last_sqrtkT = p % sqrtkT + + ! Get distributed offset + if (size(c % material) > 1 .or. size(c % sqrtkT) > 1) then + ! Distributed instances of this cell have different + ! materials/temperatures. Determine which instance this is for + ! assigning the matching material/temperature. distribcell_index = c % distribcell_index offset = 0 do k = 1, p % n_coord @@ -276,37 +277,20 @@ contains end if end if end do - p % material = c % material(offset + 1) end if - ! Set the particle temperature - if (size(c % sqrtkT) == 1) then - ! Only one temperature for this cell; assign that one to the particle. - p % sqrtkT = c % sqrtkT(1) + ! Save the material + if (size(c % material) > 1) then + p % material = c % material(offset + 1) else - ! Distributed instances of this cell have different temperatures. - ! Determine which instance this is and assign the matching temp. - distribcell_index = c % distribcell_index - offset = 0 - do k = 1, p % n_coord - if (cells(p % coord(k) % cell) % type == CELL_FILL) then - offset = offset + cells(p % coord(k) % cell) % & - offset(distribcell_index) - elseif (cells(p % coord(k) % cell) % type == CELL_LATTICE) then - if (lattices(p % coord(k + 1) % lattice) % obj & - % are_valid_indices([& - p % coord(k + 1) % lattice_x, & - p % coord(k + 1) % lattice_y, & - p % coord(k + 1) % lattice_z])) then - offset = offset + lattices(p % coord(k + 1) % lattice) % obj % & - offset(distribcell_index, & - p % coord(k + 1) % lattice_x, & - p % coord(k + 1) % lattice_y, & - p % coord(k + 1) % lattice_z) - end if - end if - end do + p % material = c % material(1) + end if + + ! Save the temperature + if (size(c % sqrtkT) > 1) then p % sqrtkT = c % sqrtkT(offset + 1) + else + p % sqrtkT = c % sqrtkT(1) end if elseif (c % type == CELL_FILL) then CELL_TYPE diff --git a/src/global.F90 b/src/global.F90 index 0fc1bc0c5..bea5f61a8 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -99,13 +99,10 @@ module global ! What to assume for expanding natural elements integer :: default_expand = ENDF_BVII1 - ! Whether or not windowed multipole cross sections should be used. - logical :: multipole_active = .false. - - ! Total amount of nuclide ZAID and dictionary of nuclide ZAID and index -- - ! this is used when sampling unresolved resonance probability tables - integer(8) :: n_nuc_zaid_total - type(DictIntInt) :: nuc_zaid_dict + ! Default temperature and method for choosing temperatures + integer :: temperature_method = TEMPERATURE_NEAREST + real(8) :: temperature_tolerance = 10.0_8 + real(8) :: temperature_default = 293.6_8 ! ============================================================================ ! MULTI-GROUP CROSS SECTION RELATED VARIABLES @@ -430,7 +427,6 @@ module global ! Various output options logical :: output_summary = .true. - logical :: output_xs = .false. logical :: output_tallies = .true. ! ============================================================================ diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index 36195d119..415ec55c9 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -73,6 +73,7 @@ module hdf5_interface module procedure read_attribute_integer_1D module procedure read_attribute_integer_2D module procedure read_attribute_string + module procedure read_attribute_string_1D end interface read_attribute interface write_attribute @@ -2425,6 +2426,67 @@ contains call h5tclose_f(memtype, hdf5_err) end subroutine read_attribute_string + subroutine read_attribute_string_1D(buffer, obj_id, name) + character(*), target, allocatable, intent(inout) :: buffer(:) + integer(HID_T), intent(in) :: obj_id + character(*), intent(in) :: name + + integer :: hdf5_err + integer(HID_T) :: space_id + integer(HID_T) :: attr_id + integer(HSIZE_T) :: dims(1) + integer(HSIZE_T) :: maxdims(1) + + call h5aopen_f(obj_id, trim(name), attr_id, hdf5_err) + + if (allocated(buffer)) then + dims(:) = shape(buffer) + else + call h5aget_space_f(attr_id, space_id, hdf5_err) + call h5sget_simple_extent_dims_f(space_id, dims, maxdims, hdf5_err) + allocate(buffer(dims(1))) + call h5sclose_f(space_id, hdf5_err) + end if + + call read_attribute_string_1D_explicit(attr_id, dims, buffer) + call h5aclose_f(attr_id, hdf5_err) + end subroutine read_attribute_string_1D + + subroutine read_attribute_string_1D_explicit(attr_id, dims, buffer) + integer(HID_T), intent(in) :: attr_id + integer(HSIZE_T), intent(in) :: dims(1) + character(*), target, intent(inout) :: buffer(dims(1)) + + integer :: hdf5_err + integer(HID_T) :: filetype + integer(HID_T) :: memtype + integer(SIZE_T) :: size + integer(SIZE_T) :: n + type(c_ptr) :: f_ptr + + ! Make sure buffer is large enough + call h5aget_type_f(attr_id, filetype, hdf5_err) + call h5tget_size_f(filetype, size, hdf5_err) + if (size > len(buffer(1)) + 1) then + print *, size, len(buffer(1)) + call fatal_error("Character buffer is not long enough to & + &read HDF5 string array.") + end if + + ! Get datatype in memory based on Fortran character + n = len(buffer(1)) + call h5tcopy_f(H5T_FORTRAN_S1, memtype, hdf5_err) + call h5tset_size_f(memtype, n, hdf5_err) + + ! Get pointer to start of string + f_ptr = c_loc(buffer(1)(1:1)) + + call h5aread_f(attr_id, memtype, f_ptr, hdf5_err) + + call h5tclose_f(filetype, hdf5_err) + call h5tclose_f(memtype, hdf5_err) + end subroutine read_attribute_string_1D_explicit + subroutine get_shape(obj_id, dims) integer(HID_T), intent(in) :: obj_id integer(HSIZE_T), intent(out) :: dims(:) diff --git a/src/initialize.F90 b/src/initialize.F90 index 99adf967f..f05b4f2b9 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -4,7 +4,7 @@ module initialize use constants use dict_header, only: DictIntInt, ElemKeyValueII use set_header, only: SetInt - use energy_grid, only: logarithmic_grid, grid_method, unionized_grid + use energy_grid, only: logarithmic_grid, grid_method ! , unionized_grid use error, only: fatal_error, warning use geometry, only: neighbor_lists, count_instance, calc_offsets, & maximum_levels @@ -17,7 +17,7 @@ module initialize use material_header, only: Material use mgxs_data, only: read_mgxs, create_macro_xs use output, only: title, header, print_version, write_message, & - print_usage, write_xs_summary, print_plot + print_usage, print_plot use random_lcg, only: initialize_prng use state_point, only: load_state_point use string, only: to_str, starts_with, ends_with, str_to_int @@ -111,20 +111,8 @@ contains if (run_mode /= MODE_PLOTTING) then ! Construct information needed for nuclear data if (run_CE) then - ! Set undefined cell temperatures to match the material data. - call lookup_material_temperatures() - - ! Construct unionized or log energy grid for cross-sections - select case (grid_method) - case (GRID_NUCLIDE) - continue - case (GRID_MAT_UNION) - call time_unionize%start() - call unionized_grid() - call time_unionize%stop() - case (GRID_LOGARITHM) - call logarithmic_grid() - end select + ! Construct log energy grid for cross-sections + call logarithmic_grid() else ! Create material macroscopic data for MGXS call time_read_xs%start() @@ -158,9 +146,6 @@ contains else ! Write summary information if (output_summary) call write_summary() - - ! Write cross section information - if (output_xs) call write_xs_summary() end if end if @@ -1005,57 +990,4 @@ contains end subroutine allocate_offsets -!=============================================================================== -! LOOKUP_MATERIAL_TEMPERATURES If any cells have undefined temperatures, try to -! find their temperatures from material data. -!=============================================================================== - - subroutine lookup_material_temperatures() - integer :: i, j, k - real(8) :: min_temp - logical :: warning_given - - warning_given = .false. - do i = 1, n_cells - ! Ignore non-normal cells and cells with defined temperature. - if (cells(i) % type /= CELL_NORMAL) cycle - if (cells(i) % sqrtkT(1) /= ERROR_REAL) cycle - - ! Set the number of temperatures equal to the number of materials. - deallocate(cells(i) % sqrtkT) - allocate(cells(i) % sqrtkT(size(cells(i) % material))) - - ! Check each of the cell materials for temperature data. - do j = 1, size(cells(i) % material) - ! Arbitrarily set void regions to 0K. - if (cells(i) % material(j) == MATERIAL_VOID) then - cells(i) % sqrtkT(j) = ZERO - cycle - end if - - associate (mat => materials(cells(i) % material(j))) - ! Find the temperature of the coldest nuclide. - min_temp = nuclides(mat % nuclide(1)) % kT - do k = 2, mat % n_nuclides - ! Warn the user if the nuclides don't have identical temperatues. - if (nuclides(mat % nuclide(k)) % kT /= min_temp & - .and. .not. warning_given .and. multipole_active) then - call warning("OpenMC cannot & - &identify the temperature of at least one cell. For the & - &purposes of multipole cross section evaluations, all cells & - &with unknown temperature will be set to the coldest & - &temperature found in the nuclear data for that cell's & - &material") - warning_given = .true. - end if - min_temp = min(min_temp, nuclides(mat % nuclide(k)) % kT) - end do - - ! Set the temperature for this cell instance. - cells(i) % sqrtkT(j) = sqrt(min_temp) - end associate - end do - end do - end subroutine lookup_material_temperatures - end module initialize diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 1cb0dc0fd..6ef8751a0 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2,6 +2,7 @@ module input_xml use hdf5 + use algorithm, only: find use cmfd_input, only: configure_cmfd use constants use dict_header, only: DictIntInt, ElemKeyValueCI @@ -23,7 +24,8 @@ module input_xml use set_header, only: SetChar use stl_vector, only: VectorInt, VectorReal, VectorChar use string, only: to_lower, to_str, str_to_int, str_to_real, & - starts_with, ends_with, tokenize, split_string + starts_with, ends_with, tokenize, split_string, & + zero_padded use tally_header, only: TallyObject use tally_filter use tally_initialize, only: add_tallies @@ -359,26 +361,6 @@ contains ! Copy random number seed if specified if (check_for_node(doc, "seed")) call get_node_value(doc, "seed", seed) - ! Energy grid methods - if (check_for_node(doc, "energy_grid")) then - call get_node_value(doc, "energy_grid", temp_str) - else - temp_str = 'logarithm' - end if - select case (trim(temp_str)) - case ('nuclide') - grid_method = GRID_NUCLIDE - case ('material-union', 'union') - grid_method = GRID_MAT_UNION - if (trim(temp_str) == 'union') & - call warning('Energy grids will be unionized by material. Global& - & energy grid unionization is no longer an allowed option.') - case ('logarithm', 'logarithmic', 'log') - grid_method = GRID_LOGARITHM - case default - call fatal_error("Unknown energy grid method: " // trim(temp_str)) - end select - ! Number of bins for logarithmic grid if (check_for_node(doc, "log_grid_bins")) then call get_node_value(doc, "log_grid_bins", n_log_bins) @@ -978,14 +960,6 @@ contains trim(temp_str) == '0') output_summary = .false. end if - ! Check for cross sections option - if (check_for_node(node_output, "cross_sections")) then - call get_node_value(node_output, "cross_sections", temp_str) - temp_str = to_lower(temp_str) - if (trim(temp_str) == 'true' .or. & - trim(temp_str) == '1') output_xs = .true. - end if - ! Check for ASCII tallies output option if (check_for_node(node_output, "tallies")) then call get_node_value(node_output, "tallies", temp_str) @@ -1089,20 +1063,6 @@ contains end select end if - ! Check to see if windowed multipole functionality is requested - if (check_for_node(doc, "use_windowed_multipole")) then - call get_node_value(doc, "use_windowed_multipole", temp_str) - select case (to_lower(temp_str)) - case ('true', '1') - multipole_active = .true. - case ('false', '0') - multipole_active = .false. - case default - call fatal_error("Unrecognized value for in & - &settings.xml") - end select - end if - call get_node_list(doc, "volume_calc", node_vol_list) n = get_list_size(node_vol_list) allocate(volume_calcs(n)) @@ -1111,6 +1071,27 @@ contains call volume_calcs(i) % from_xml(node_vol) end do + ! Get temperature settings + if (check_for_node(doc, "temperature_default")) then + call get_node_value(doc, "temperature_default", temperature_default) + end if + if (check_for_node(doc, "temperature_method")) then + call get_node_value(doc, "temperature_method", temp_str) + select case (to_lower(temp_str)) + case ('nearest') + temperature_method = TEMPERATURE_NEAREST + case ('interpolation') + temperature_method = TEMPERATURE_INTERPOLATION + case ('multipole') + temperature_method = TEMPERATURE_MULTIPOLE + case default + call fatal_error("Unknown temperature method: " // trim(temp_str)) + end select + end if + if (check_for_node(doc, "temperature_tolerance")) then + call get_node_value(doc, "temperature_tolerance", temperature_tolerance) + end if + ! Close settings XML file call close_xmldoc(doc) @@ -2050,6 +2031,9 @@ contains integer :: i, j type(DictCharInt) :: library_dict type(Library), allocatable :: libraries(:) + type(VectorReal), allocatable :: nuc_temps(:) ! List of T to read for each nuclide + type(VectorReal), allocatable :: sab_temps(:) ! List of T to read for each S(a,b) + real(8), allocatable :: material_temps(:) if (run_CE) then call read_ce_cross_sections_xml(libraries) @@ -2076,12 +2060,18 @@ contains end if ! Parse data from materials.xml - call read_materials_xml(libraries, library_dict) + call read_materials_xml(libraries, library_dict, material_temps) + + ! Assign temperatures to cells that don't have temperatures already assigned + call assign_temperatures(material_temps) + + ! Determine desired temperatures for each nuclide and S(a,b) table + call get_temperatures(nuc_temps, sab_temps) ! Read continuous-energy cross sections if (run_CE .and. run_mode /= MODE_PLOTTING) then call time_read_xs%start() - call read_ce_cross_sections(libraries, library_dict) + call read_ce_cross_sections(libraries, library_dict, nuc_temps, sab_temps) call time_read_xs%stop() end if @@ -2092,9 +2082,10 @@ contains call library_dict % clear() end subroutine read_materials - subroutine read_materials_xml(libraries, library_dict) + subroutine read_materials_xml(libraries, library_dict, material_temps) type(Library), intent(in) :: libraries(:) type(DictCharInt), intent(inout) :: library_dict + real(8), allocatable, intent(out) :: material_temps(:) integer :: i ! loop index for materials integer :: j ! loop index for nuclides @@ -2110,7 +2101,6 @@ contains logical :: file_exists ! does materials.xml exist? logical :: sum_density ! density is taken to be sum of nuclide densities character(20) :: name ! name of isotope, e.g. 92235.03c - character(6) :: default_temperature ! Default temperature, e.g., '300K' character(MAX_WORD_LEN) :: units ! units on density character(MAX_LINE_LEN) :: filename ! absolute path to materials.xml character(MAX_LINE_LEN) :: temp_str ! temporary string when reading @@ -2144,25 +2134,13 @@ contains ! Parse materials.xml file call open_xmldoc(doc, filename) - ! Copy default temperature - if (check_for_node(doc, "default_temperature")) then - call get_node_value(doc, "default_temperature", default_temperature) - else if (.not. run_CE) then - ! FIXME This is only necessary while MG mode does not have a - ! temperature dependent library implementation. - ! Set a default for MG mode to allow MG libraries to not include - ! temperatures - default_temperature = '294K' - else - default_temperature = '' - end if - ! Get pointer to list of XML call get_node_list(doc, "material", node_mat_list) ! Allocate cells array n_materials = get_list_size(node_mat_list) allocate(materials(n_materials)) + allocate(material_temps(n_materials)) ! Initialize count for number of nuclides/S(a,b) tables index_nuclide = 0 @@ -2192,14 +2170,11 @@ contains call get_node_value(node_mat, "name", mat % name) end if - ! Copy material temperature + ! Get material default temperature if (check_for_node(node_mat, "temperature")) then - call get_node_value(node_mat, "temperature", mat % temperature) - else if (default_temperature /= '') then - mat % temperature = default_temperature + call get_node_value(node_mat, "temperature", material_temps(i)) else - call fatal_error("Must specify either a material temperature or a & - &default temperature") + material_temps(i) = ERROR_REAL end if ! ======================================================================= @@ -5719,10 +5694,10 @@ contains ASSIGN_SAB: do k = 1, size(mat % i_sab_tables) ! In order to know which nuclide the S(a,b) table applies to, we need ! to search through the list of nuclides for one which has a matching - ! zaid + ! name associate (sab => sab_tables(mat % i_sab_tables(k))) FIND_NUCLIDE: do j = 1, size(mat % nuclide) - if (any(sab % zaid == nuclides(mat % nuclide(j)) % zaid)) then + if (any(sab % nuclides == nuclides(mat % nuclide(j)) % name)) then mat % i_sab_nuclides(k) = j exit FIND_NUCLIDE end if @@ -5774,16 +5749,16 @@ contains end do end subroutine assign_sab_tables - subroutine read_ce_cross_sections(libraries, library_dict) + subroutine read_ce_cross_sections(libraries, library_dict, nuc_temps, sab_temps) type(Library), intent(in) :: libraries(:) type(DictCharInt), intent(inout) :: library_dict + type(VectorReal), intent(in) :: nuc_temps(:) + type(VectorReal), intent(in) :: sab_temps(:) integer :: i, j integer :: i_library integer :: i_nuclide integer :: i_sab - integer :: index_nuc_zaid ! index in nuclide ZAID - integer :: zaid ! ZAID of nuclide integer(HID_T) :: file_id integer(HID_T) :: group_id logical :: mp_found ! if windowed multipole libraries were found @@ -5796,8 +5771,6 @@ contains allocate(micro_xs(n_nuclides_total)) !$omp end parallel - index_nuc_zaid = 0 - ! Read cross sections do i = 1, size(materials) do j = 1, size(materials(i) % names) @@ -5813,8 +5786,8 @@ contains ! Read nuclide data from HDF5 file_id = file_open(libraries(i_library) % path, 'r') group_id = open_group(file_id, name) - call nuclides(i_nuclide) % from_hdf5(group_id, & - materials(i) % temperature) + call nuclides(i_nuclide) % from_hdf5(group_id, nuc_temps(i_nuclide), & + temperature_method, temperature_tolerance) call close_group(group_id) call file_close(file_id) @@ -5824,23 +5797,19 @@ contains ! Determine if minimum/maximum energy for this nuclide is greater/less ! than the previous - energy_min_neutron = max(energy_min_neutron, nuclides(i_nuclide) % energy(1)) - energy_max_neutron = min(energy_max_neutron, nuclides(i_nuclide) % energy(& - size(nuclides(i_nuclide) % energy))) + if (size(nuclides(i_nuclide) % grid) >= 1) then + energy_min_neutron = max(energy_min_neutron, & + nuclides(i_nuclide) % grid(1) % energy(1)) + energy_max_neutron = min(energy_max_neutron, nuclides(i_nuclide) % & + grid(1) % energy(size(nuclides(i_nuclide) % grid(1) % energy))) + end if ! Add name and alias to dictionary call already_read % add(name) - ! Construct dictionary mapping nuclide zaids to [1,N] -- used for - ! unresolved resonance probability tables - zaid = nuclides(i_nuclide) % zaid - if (.not. nuc_zaid_dict % has_key(zaid)) then - index_nuc_zaid = index_nuc_zaid + 1 - call nuc_zaid_dict % add_key(zaid, index_nuc_zaid) - end if - ! Read multipole file into the appropriate entry on the nuclides array - if (multipole_active) call read_multipole_data(i_nuclide) + if (temperature_method == TEMPERATURE_MULTIPOLE) & + call read_multipole_data(i_nuclide) end if ! Check if material is fissionable @@ -5868,8 +5837,8 @@ contains ! Read S(a,b) data from HDF5 file_id = file_open(libraries(i_library) % path, 'r') group_id = open_group(file_id, name) - call sab_tables(i_sab) % from_hdf5(group_id, & - materials(i) % temperature) + call sab_tables(i_sab) % from_hdf5(group_id, sab_temps(i_sab), & + temperature_tolerance) call close_group(group_id) call file_close(file_id) @@ -5879,14 +5848,13 @@ contains end do end do - n_nuc_zaid_total = index_nuc_zaid - ! Associate S(a,b) tables with specific nuclides call assign_sab_tables() ! Show which nuclide results in lowest energy for neutron transport do i = 1, size(nuclides) - if (nuclides(i) % energy(nuclides(i) % n_grid) == energy_max_neutron) then + if (nuclides(i) % grid(1) % energy(size(nuclides(i) % grid(1) % energy)) & + == energy_max_neutron) then call write_message("Maximum neutron transport energy: " // & trim(to_str(energy_max_neutron)) // " MeV for " // & trim(adjustl(nuclides(i) % name)), 6) @@ -5895,7 +5863,7 @@ contains end do ! If the user wants multipole, make sure we found a multipole library. - if (multipole_active) then + if (temperature_method == TEMPERATURE_MULTIPOLE) then mp_found = .false. do i = 1, size(nuclides) if (nuclides(i) % mp_present) then @@ -5911,6 +5879,107 @@ contains end subroutine read_ce_cross_sections +!=============================================================================== +! ASSIGN_TEMPERATURES If any cells have undefined temperatures, try to find +! their temperatures from material or global default temperatures +!=============================================================================== + + subroutine assign_temperatures(material_temps) + real(8), intent(in) :: material_temps(:) + + integer :: i, j + integer :: i_material + + do i = 1, n_cells + ! Ignore non-normal cells and cells with defined temperature. + if (cells(i) % material(1) == NONE) cycle + if (cells(i) % sqrtkT(1) /= ERROR_REAL) cycle + + ! Set the number of temperatures equal to the number of materials. + deallocate(cells(i) % sqrtkT) + allocate(cells(i) % sqrtkT(size(cells(i) % material))) + + ! Check each of the cell materials for temperature data. + do j = 1, size(cells(i) % material) + ! Arbitrarily set void regions to 0K. + if (cells(i) % material(j) == MATERIAL_VOID) then + cells(i) % sqrtkT(j) = ZERO + cycle + end if + + ! Use material default or global default temperature + i_material = material_dict % get_key(cells(i) % material(j)) + if (material_temps(i_material) /= ERROR_REAL) then + cells(i) % sqrtkT(j) = sqrt(K_BOLTZMANN * & + material_temps(i_material)) + else + cells(i) % sqrtkT(j) = sqrt(K_BOLTZMANN * temperature_default) + end if + end do + end do + end subroutine assign_temperatures + +!=============================================================================== +! GET_TEMPERATURES returns a list of temperatures that each nuclide/S(a,b) table +! appears at in the model. Later, this list is used to determine the actual +! temperatures to read (which may be different if interpolation is used) +!=============================================================================== + + subroutine get_temperatures(nuc_temps, sab_temps) + type(VectorReal), allocatable, intent(out) :: nuc_temps(:) + type(VectorReal), allocatable, intent(out) :: sab_temps(:) + + integer :: i, j, k + integer :: i_nuclide ! index in nuclides array + integer :: i_sab ! index in S(a,b) array + integer :: i_material + real(8) :: temperature ! temperature in Kelvin + + allocate(nuc_temps(n_nuclides_total)) + allocate(sab_temps(n_sab_tables)) + + 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 + + ! Get temperature of cell (rounding to nearest integer) + if (size(cells(i) % sqrtkT) > 1) then + temperature = cells(i) % sqrtkT(j)**2 / K_BOLTZMANN + else + temperature = cells(i) % sqrtkT(1)**2 / K_BOLTZMANN + end if + + 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 nuc_temps array + i_nuclide = nuclide_dict % get_key(to_lower(mat % names(k))) + + ! Add temperature if it hasn't already been added + if (find(nuc_temps(i_nuclide), temperature) == -1) then + call nuc_temps(i_nuclide) % push_back(temperature) + end if + end do NUC_NAMES_LOOP + + if (mat % n_sab > 0) then + SAB_NAMES_LOOP: do k = 1, size(mat % sab_names) + ! Get index in nuc_temps array + i_sab = sab_dict % get_key(to_lower(mat % sab_names(k))) + + ! Add temperature if it hasn't already been added + if (find(sab_temps(i_sab), temperature) == -1) then + call sab_temps(i_sab) % push_back(temperature) + end if + end do SAB_NAMES_LOOP + end if + end associate + end do + end do + + end subroutine get_temperatures + !=============================================================================== ! READ_0K_ELASTIC_SCATTERING !=============================================================================== @@ -5927,6 +5996,9 @@ contains real(8) :: xs_cdf_sum character(MAX_WORD_LEN) :: name type(Nuclide) :: resonant_nuc + type(VectorReal) :: temperature + + call temperature % push_back(ZERO) do i = 1, size(nuclides_0K) if (nuc % name == nuclides_0K(i) % nuclide) then @@ -5946,13 +6018,14 @@ contains ! Read nuclide data from HDF5 file_id = file_open(libraries(i_library) % path, 'r') group_id = open_group(file_id, name) - call resonant_nuc % from_hdf5(group_id, '0K') + call resonant_nuc % from_hdf5(group_id, temperature, & + TEMPERATURE_NEAREST, 1000.0_8) call close_group(group_id) call file_close(file_id) ! Copy 0K energy grid and elastic scattering cross section - call move_alloc(TO=nuc % energy_0K, FROM=resonant_nuc % energy) - call move_alloc(TO=nuc % elastic_0K, FROM=resonant_nuc % elastic) + call move_alloc(TO=nuc % energy_0K, FROM=resonant_nuc % grid(1) % energy) + call move_alloc(TO=nuc % elastic_0K, FROM=resonant_nuc % sum_xs(1) % elastic) nuc % n_grid_0K = size(nuc % energy_0K) ! Build CDF for 0K elastic scattering @@ -5987,18 +6060,22 @@ contains integer, intent(in) :: i_table ! index in nuclides/sab_tables - integer :: i logical :: file_exists ! Does multipole library exist? character(7) :: readable ! Is multipole library readable? - character(6) :: zaid_string ! String of the ZAID - character(MAX_FILE_LEN+9) :: filename ! Path to multipole xs library + character(MAX_FILE_LEN) :: filename ! Path to multipole xs library ! For the time being, and I know this is a bit hacky, we just assume - ! that the file will be zaid.h5. + ! that the file will be ZZZAAAmM.h5. associate (nuc => nuclides(i_table)) - write(zaid_string, '(I6.6)') nuc % zaid - filename = trim(path_multipole) // zaid_string // ".h5" + if (nuc % metastable > 0) then + filename = trim(path_multipole) // trim(zero_padded(nuc % Z, 3)) // & + trim(zero_padded(nuc % A, 3)) // 'm' // & + trim(to_str(nuc % metastable)) // ".h5" + else + filename = trim(path_multipole) // trim(zero_padded(nuc % Z, 3)) // & + trim(zero_padded(nuc % A, 3)) // ".h5" + end if ! Check if Multipole library exists and is readable inquire(FILE=filename, EXIST=file_exists, READ=readable) @@ -6019,16 +6096,6 @@ contains call multipole_read(filename, nuc % multipole, i_table) nuc % mp_present = .true. - ! Recreate nu-fission cross section - if (nuc % fissionable) then - do i = 1, size(nuc % energy) - nuc % nu_fission(i) = nuc % nu(nuc % energy(i), EMISSION_TOTAL) * & - nuc % fission(i) - end do - else - nuc % nu_fission(:) = ZERO - end if - end associate end subroutine read_multipole_data diff --git a/src/material_header.F90 b/src/material_header.F90 index 27a0e61fe..772e3a415 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -13,9 +13,6 @@ module material_header integer, allocatable :: nuclide(:) ! index in nuclides array real(8) :: density ! total atom density in atom/b-cm real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm - character(6) :: temperature ! Temperature of the material - ! as presented in the HDF5 library; - ! e.g., "300K" ! Energy grid information integer :: n_grid ! # of union material grid points diff --git a/src/mgxs_header.F90 b/src/mgxs_header.F90 index 06dd1e213..5947d2562 100644 --- a/src/mgxs_header.F90 +++ b/src/mgxs_header.F90 @@ -19,7 +19,6 @@ module mgxs_header type, abstract :: Mgxs character(len=104) :: name ! name of dataset, e.g. 92235.03c - integer :: zaid ! Z and A identifier, e.g. 92235 real(8) :: awr ! Atomic Weight Ratio real(8) :: kT ! temperature in MeV (k*T) @@ -211,11 +210,6 @@ module mgxs_header else this % kT = ZERO end if - if (check_for_node(node_xsdata, "zaid")) then - call get_node_value(node_xsdata, "zaid", this % zaid) - else - this % zaid = 0 - end if if (check_for_node(node_xsdata, "awr")) then call get_node_value(node_xsdata, "awr", this % awr) else @@ -970,11 +964,6 @@ module mgxs_header ! Basic nuclide information write(unit_,*) 'MGXS Entry: ' // trim(this % name) - if (this % zaid > 0) then - write(unit_,*) ' ZAID = ' // trim(to_str(this % zaid)) - else if (this % zaid < 0) then - write(unit_,*) ' Material id = ' // trim(to_str(-this % zaid)) - end if if (this % awr > ZERO) then write(unit_,*) ' AWR = ' // trim(to_str(this % awr)) end if @@ -1319,7 +1308,6 @@ module mgxs_header else this % name = mat % name end if - this % zaid = -mat % id this % fissionable = mat % fissionable this % scatt_type = scatt_type diff --git a/src/multipole.F90 b/src/multipole.F90 index 770121b9d..a099e5047 100644 --- a/src/multipole.F90 +++ b/src/multipole.F90 @@ -28,13 +28,8 @@ contains integer(HID_T) :: group_id ! Intermediate loading components - character(len=10) :: version - integer :: NMT - integer :: i, j - integer, allocatable :: MT(:) - logical :: accumulated_fission - character(len=24) :: MT_n ! Takes the form '/nuclide/reactions/MT???' integer :: is_fissionable + character(len=10) :: version associate (nuc => nuclides(i_table)) @@ -80,111 +75,8 @@ contains call read_dataset(multipole % curvefit, group_id, "curvefit") - ! Delete ACE pointwise data - call read_dataset(nuc % n_grid, group_id, "n_grid") - - deallocate(nuc % energy) - deallocate(nuc % total) - deallocate(nuc % elastic) - deallocate(nuc % fission) - deallocate(nuc % nu_fission) - deallocate(nuc % absorption) - - allocate(nuc % energy(nuc % n_grid)) - allocate(nuc % total(nuc % n_grid)) - allocate(nuc % elastic(nuc % n_grid)) - allocate(nuc % fission(nuc % n_grid)) - allocate(nuc % nu_fission(nuc % n_grid)) - allocate(nuc % absorption(nuc % n_grid)) - - nuc % total(:) = ZERO - nuc % absorption(:) = ZERO - nuc % fission(:) = ZERO - - ! Read in new energy axis (converting eV to MeV) - call read_dataset(nuc % energy, group_id, "energy_points") - nuc % energy = nuc % energy / 1.0e6_8 - - ! Get count and list of MT tables - call read_dataset(NMT, group_id, "MT_count") - allocate(MT(NMT)) - - call read_dataset(MT, group_id, "MT_list") - call close_group(group_id) - accumulated_fission = .false. - - ! Loop over each MT entry and load it into a reaction. - do i = 1, NMT - write(MT_n, '(A, I3.3)') '/nuclide/reactions/MT', MT(i) - - group_id = open_group(file_id, MT_n) - - ! Each MT needs to be treated slightly differently. - select case (MT(i)) - case(ELASTIC) - call read_dataset(nuc % elastic, group_id, "MT_sigma") - nuc % total(:) = nuc % total + nuc % elastic - case(N_FISSION) - call read_dataset(nuc % fission, group_id, "MT_sigma") - nuc % total(:) = nuc % total + nuc % fission - nuc % absorption(:) = nuc % absorption + nuc % fission - accumulated_fission = .true. - case default - ! Search through all of our secondary reactions - do j = 1, size(nuc % reactions) - if (nuc % reactions(j) % MT == MT(i)) then - ! Match found - - ! Individual Fission components exist, so remove the combined - ! fission cross section. - if ( (MT(i) == N_F .or. MT(i) == N_NF .or. MT(i) == N_2NF & - .or. MT(i) == N_3NF) .and. accumulated_fission) then - nuc % total(:) = nuc % total - nuc % fission - nuc % absorption(:) = nuc % absorption - nuc % fission - nuc % fission(:) = ZERO - accumulated_fission = .false. - end if - - deallocate(nuc % reactions(j) % sigma) - allocate(nuc % reactions(j) % sigma(nuc % n_grid)) - - call read_dataset(nuc % reactions(j) % sigma, & - group_id, "MT_sigma") - call read_dataset(nuc % reactions(j) % Q_value, & - group_id, "Q_value") - call read_dataset(nuc % reactions(j) % threshold, & - group_id, "threshold") - nuc % reactions(j) % threshold = 1 ! TODO: reconsider implications. - nuc % reactions(j) % Q_value = nuc % reactions(j) % Q_value & - / 1.0e6_8 - - ! Accumulate total - if (MT(i) /= N_LEVEL .and. MT(i) <= N_DA) then - nuc % total(:) = nuc % total + nuc % reactions(j) % sigma - end if - - ! Accumulate absorption - if (MT(i) >= N_GAMMA .and. MT(i) <= N_DA) then - nuc % absorption(:) = nuc % absorption & - + nuc % reactions(j) % sigma - end if - - ! Accumulate fission (if needed) - if ( (MT(i) == N_F .or. MT(i) == N_NF .or. MT(i) == N_2NF & - .or. MT(i) == N_3NF) ) then - nuc % fission(:) = nuc % fission + nuc % reactions(j) % sigma - nuc % absorption(:) = nuc % absorption & - + nuc % reactions(j) % sigma - end if - end if - end do - end select - - call close_group(group_id) - end do - ! Close file call file_close(file_id) diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 3fffb76f3..f9e473f45 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -7,6 +7,7 @@ module nuclide_header h5lget_name_by_idx_f, H5_INDEX_NAME_F, H5_ITER_INC_F use h5lt, only: h5ltpath_valid_f + use algorithm, only: sort, find use constants use dict_header, only: DictIntInt use endf, only: reaction_name, is_fission, is_disappearance @@ -20,7 +21,7 @@ module nuclide_header use product_header, only: AngleEnergyContainer use reaction_header, only: Reaction use secondary_uncorrelated, only: UncorrelatedAngleEnergy - use stl_vector, only: VectorInt + use stl_vector, only: VectorInt, VectorReal use string use urr_header, only: UrrData use xml_interface @@ -32,29 +33,37 @@ module nuclide_header ! for continuous-energy neutron transport. !=============================================================================== - type :: Nuclide - ! Nuclide meta-data - character(20) :: name ! name of nuclide, e.g. U235.71c - integer :: zaid ! Z and A identifier, e.g. 92235 - integer :: metastable ! metastable state - real(8) :: awr ! Atomic Weight Ratio - real(8) :: kT ! temperature in MeV (k*T) - - ! Fission information - logical :: fissionable = .false. ! nuclide is fissionable? - - ! Energy grid information - integer :: n_grid ! # of nuclide grid points + type EnergyGrid integer, allocatable :: grid_index(:) ! log grid mapping indices real(8), allocatable :: energy(:) ! energy values corresponding to xs + end type EnergyGrid - ! Microscopic cross sections + type SumXS real(8), allocatable :: total(:) ! total cross section real(8), allocatable :: elastic(:) ! elastic scattering real(8), allocatable :: fission(:) ! fission real(8), allocatable :: nu_fission(:) ! neutron production real(8), allocatable :: absorption(:) ! absorption (MT > 100) real(8), allocatable :: heating(:) ! heating + end type SumXS + + type :: Nuclide + ! Nuclide meta-data + character(20) :: name ! name of nuclide, e.g. U235.71c + integer :: Z ! atomic number + integer :: A ! mass number + integer :: metastable ! metastable state + real(8) :: awr ! Atomic Weight Ratio + real(8), allocatable :: kTs(:) ! temperature in MeV (k*T) + + ! Fission information + logical :: fissionable = .false. ! nuclide is fissionable? + + ! Energy grid for each temperature + type(EnergyGrid), allocatable :: grid(:) + + ! Microscopic cross sections + type(SumXS), allocatable :: sum_xs(:) ! Resonance scattering info logical :: resonant = .false. ! resonant scatterer? @@ -77,7 +86,7 @@ module nuclide_header ! Unresolved resonance data logical :: urr_present = .false. integer :: urr_inelastic - type(UrrData), pointer :: urr_data => null() + type(UrrData), allocatable :: urr_data(:) ! Multipole data logical :: mp_present = .false. @@ -94,7 +103,6 @@ module nuclide_header contains procedure :: clear => nuclide_clear - procedure :: print => nuclide_print procedure :: from_hdf5 => nuclide_from_hdf5 procedure :: nu => nuclide_nu procedure, private :: create_derived => nuclide_create_derived @@ -131,6 +139,7 @@ module nuclide_header ! Information for S(a,b) use integer :: index_sab ! index in sab_tables (zero means no table) integer :: last_index_sab = 0 ! index in sab_tables last used by this nuclide + integer :: index_temp_sab ! temperature index for sab_tables real(8) :: elastic_sab ! microscopic elastic scattering on S(a,b) table ! Information for URR probability table use @@ -173,23 +182,25 @@ module nuclide_header subroutine nuclide_clear(this) class(Nuclide), intent(inout) :: this ! The Nuclide object to clear - if (associated(this % urr_data)) deallocate(this % urr_data) if (associated(this % multipole)) deallocate(this % multipole) end subroutine nuclide_clear - subroutine nuclide_from_hdf5(this, group_id, temperature) - class(Nuclide), intent(inout) :: this - integer(HID_T), intent(in) :: group_id - character(len=*), intent(in) :: temperature + subroutine nuclide_from_hdf5(this, group_id, temperature, method, tolerance) + class(Nuclide), intent(inout) :: this + integer(HID_T), intent(in) :: group_id + type(VectorReal), intent(in) :: temperature ! list of desired temperatures + integer, intent(in) :: method + real(8), intent(in) :: tolerance integer :: i - integer :: Z - integer :: A + integer :: n integer :: storage_type integer :: max_corder integer :: n_links integer :: hdf5_err + integer :: i_closest + integer :: n_temperature integer(HID_T) :: urr_group, nu_group integer(HID_T) :: energy_group, energy_dset integer(HID_T) :: kT_group, kT_dset @@ -201,13 +212,14 @@ module nuclide_header integer(SIZE_T) :: name_len, name_file_len integer(HSIZE_T) :: j integer(HSIZE_T) :: dims(1) - character(MAX_WORD_LEN) :: temp - character(MAX_FILE_LEN), allocatable :: temperatures(:) - integer, allocatable :: temperatures_integer(:) - character(6) :: my_temperature - integer :: temperature_integer - type(VectorInt) :: MTs + character(MAX_WORD_LEN) :: temp_str + character(MAX_FILE_LEN), allocatable :: dset_names(:) + real(8), allocatable :: temps_available(:) ! temperatures available + real(8) :: temp_desired + real(8) :: temp_actual logical :: exists + type(VectorInt) :: MTs + type(VectorInt) :: temps_to_read ! Get name of nuclide from group name_len = len(this % name) @@ -216,62 +228,90 @@ module nuclide_header ! Get rid of leading '/' this % name = trim(this % name(2:)) - call read_attribute(Z, group_id, 'Z') - call read_attribute(A, group_id, 'A') + call read_attribute(this % Z, group_id, 'Z') + call read_attribute(this % A, group_id, 'A') call read_attribute(this % metastable, group_id, 'metastable') - this % zaid = 1000 * Z + A + 400 * this % metastable call read_attribute(this % awr, group_id, 'atomic_weight_ratio') kT_group = open_group(group_id, 'kTs') - ! Before accessing the temperature data, see if the user-provied temperature - ! exists. We can find this out by looking at the datasets within kT_group - temperature_integer = & - str_to_int(temperature(1: len_trim(adjustl(temperature)) - 1)) - call get_datasets(kT_group, temperatures) - allocate(temperatures_integer(size(temperatures))) - do i = 1, size(temperatures) - temperatures_integer(i) = & - str_to_int(temperatures(i)(1: len_trim(adjustl(temperatures(i))) - 1)) + ! Determine temperatures available + call get_datasets(kT_group, dset_names) + allocate(temps_available(size(dset_names))) + do i = 1, size(dset_names) + ! Read temperature value + call read_dataset(temps_available(i), kT_group, trim(dset_names(i))) + temps_available(i) = temps_available(i) / K_BOLTZMANN + end do + + select case (method) + case (TEMPERATURE_NEAREST) + ! Determine actual temperatures to read + TEMP_LOOP: do i = 1, temperature % size() + temp_desired = temperature % data(i) + i_closest = minloc(abs(temps_available - temp_desired), dim=1) + temp_actual = temps_available(i_closest) + if (abs(temp_actual - temp_desired) < tolerance) then + if (find(temps_to_read, nint(temp_actual)) == -1) then + call temps_to_read % push_back(nint(temp_actual)) + + ! Write warning for resonance scattering data if 0K is not available + if (abs(temp_actual - temp_desired) > 0 .and. temp_desired == 0) then + call warning(trim(this % name) // " does not contain 0K data & + &needed for resonance scattering options selected. Using & + &data at " // trim(to_str(nint(temp_actual))) // " K instead.") + end if + end if + else + call fatal_error("Nuclear data library does not contain cross sections & + &for " // trim(this % name) // " at or near " // & + trim(to_str(nint(temp_desired))) // " K.") + end if + end do TEMP_LOOP + + case (TEMPERATURE_INTERPOLATION) + ! TODO: Get bounding temperatures + call fatal_error("Temperature interpolation not yet implemented") + + case (TEMPERATURE_MULTIPOLE) + ! Add first available temperature + call temps_to_read % push_back(nint(temps_available(1))) + + end select + + ! Sort temperatures to read + call sort(temps_to_read) + + n_temperature = temps_to_read % size() + allocate(this % kTs(n_temperature)) + allocate(this % grid(n_temperature)) + + do i = 1, n_temperature + ! Get temperature as a string + temp_str = trim(to_str(temps_to_read % data(i))) // "K" + + ! Read exact temperature value + call read_dataset(this % kTs(i), kT_group, trim(temp_str)) + + ! Read energy grid + energy_group = open_group(group_id, 'energy') + energy_dset = open_dataset(energy_group, temp_str) + call get_shape(energy_dset, dims) + allocate(this % grid(i) % energy(int(dims(1), 4))) + call read_dataset(this % grid(i) % energy, energy_dset) + call close_dataset(energy_dset) + call close_group(energy_group) end do - my_temperature = & - temperatures(minloc(abs(temperature_integer - temperatures_integer), & - dim=1)) - ! Now print a warning if there is no matching temperature and then use the - ! closest temperature - if (temperature /= my_temperature) then - if (temperature == '0K') then - call warning(trim(this % name) // " does not contain 0K data & - &needed for the resonance scattering options selected") - else - call warning(trim(this % name) // " does not contain data at a & - &temperature of " // trim(temperature) // "; using the & - &nearest available temperature of " // trim(my_temperature)) - end if - end if - kT_dset = open_dataset(kT_group, my_temperature) - call read_dataset(this % kT, kT_dset) - call close_dataset(kT_dset) call close_group(kT_group) - ! Read energy grid - energy_group = open_group(group_id, 'energy') - energy_dset = open_dataset(energy_group, my_temperature) - call get_shape(energy_dset, dims) - this % n_grid = int(dims(1), 4) - allocate(this % energy(this % n_grid)) - call read_dataset(this % energy, energy_dset) - call close_dataset(energy_dset) - call close_group(energy_group) - ! Get MT values based on group names rxs_group = open_group(group_id, 'reactions') call h5gget_info_f(rxs_group, storage_type, n_links, max_corder, hdf5_err) do j = 0, n_links - 1 call h5lget_name_by_idx_f(rxs_group, ".", H5_INDEX_NAME_F, H5_ITER_INC_F, & - j, temp, hdf5_err, name_len) - if (starts_with(temp, "reaction_")) then - call MTs % push_back(int(str_to_int(temp(10:12)))) + j, temp_str, hdf5_err, name_len) + if (starts_with(temp_str, "reaction_")) then + call MTs % push_back(int(str_to_int(temp_str(10:12)))) end if end do @@ -281,7 +321,7 @@ module nuclide_header rx_group = open_group(rxs_group, 'reaction_' // trim(& zero_padded(MTs % data(i), 3))) - call this % reactions(i) % from_hdf5(rx_group, my_temperature) + call this % reactions(i) % from_hdf5(rx_group, temps_to_read) call close_group(rx_group) end do call close_group(rxs_group) @@ -290,32 +330,42 @@ module nuclide_header call h5ltpath_valid_f(group_id, 'urr', .true., exists, hdf5_err) if (exists) then this % urr_present = .true. - allocate(this % urr_data) - urr_group = open_group(group_id, 'urr/' // trim(my_temperature)) - call this % urr_data % from_hdf5(urr_group) + allocate(this % urr_data(n_temperature)) + + do i = 1, n_temperature + ! Get temperature as a string + temp_str = trim(to_str(temps_to_read % data(i))) // "K" + + ! Read probability tables for i-th temperature + urr_group = open_group(group_id, 'urr/' // trim(temp_str)) + call this % urr_data(i) % from_hdf5(urr_group) + call close_group(urr_group) + + ! Check for negative values + if (any(this % urr_data(i) % prob < ZERO)) then + call warning("Negative value(s) found on probability table & + &for nuclide " // this % name // " at " // trim(temp_str)) + end if + end do ! if the inelastic competition flag indicates that the inelastic cross ! section should be determined from a normal reaction cross section, we ! need to get the index of the reaction - if (this % urr_data % inelastic_flag > 0) then - do i = 1, size(this % reactions) - if (this % reactions(i) % MT == this % urr_data % inelastic_flag) then - this % urr_inelastic = i + if (n_temperature > 0) then + if (this % urr_data(1) % inelastic_flag > 0) then + do i = 1, size(this % reactions) + if (this % reactions(i) % MT == this % urr_data(1) % inelastic_flag) then + this % urr_inelastic = i + end if + end do + + ! Abort if no corresponding inelastic reaction was found + if (this % urr_inelastic == NONE) then + call fatal_error("Could not find inelastic reaction specified on & + &unresolved resonance probability table.") end if - end do - - ! Abort if no corresponding inelastic reaction was found - if (this % urr_inelastic == NONE) then - call fatal_error("Could not find inelastic reaction specified on & - &unresolved resonance probability table.") end if end if - - ! Check for negative values - if (any(this % urr_data % prob < ZERO)) then - call warning("Negative value(s) found on probability table & - &for nuclide " // this % name) - end if end if ! Check for nu-total @@ -325,8 +375,8 @@ module nuclide_header ! Read total nu data total_nu = open_dataset(nu_group, 'yield') - call read_attribute(temp, total_nu, 'type') - select case (temp) + call read_attribute(temp_str, total_nu, 'type') + select case (temp_str) case ('Tabulated1D') allocate(Tabulated1D :: this % total_nu) case ('Polynomial') @@ -346,8 +396,8 @@ module nuclide_header ! Check to see if this is polynomial or tabulated data fer_dset = open_dataset(fer_group, 'q_prompt') - call read_attribute(temp, fer_dset, 'type') - if (temp == 'Polynomial') then + call read_attribute(temp_str, fer_dset, 'type') + if (temp_str == 'Polynomial') then ! Read the prompt Q-value allocate(Polynomial :: this % fission_q_prompt) call this % fission_q_prompt % from_hdf5(fer_dset) @@ -358,7 +408,7 @@ module nuclide_header fer_dset = open_dataset(fer_group, 'q_recoverable') call this % fission_q_recov % from_hdf5(fer_dset) call close_dataset(fer_dset) - else if (temp == 'Tabulated1D') then + else if (temp_str == 'Tabulated1D') then ! Read the prompt Q-value allocate(Tabulated1D :: this % fission_q_prompt) call this % fission_q_prompt % from_hdf5(fer_dset) @@ -383,108 +433,125 @@ module nuclide_header subroutine nuclide_create_derived(this) class(Nuclide), intent(inout) :: this - integer :: i - integer :: j - integer :: k + integer :: i, j, k + integer :: t integer :: m integer :: n + integer :: n_grid integer :: i_fission - type(ListInt) :: MTs + integer :: n_temperature + type(VectorInt) :: MTs - ! Allocate and initialize derived cross sections - allocate(this % total(this % n_grid)) - allocate(this % elastic(this % n_grid)) - allocate(this % fission(this % n_grid)) - allocate(this % nu_fission(this % n_grid)) - allocate(this % absorption(this % n_grid)) - this % total(:) = ZERO - this % elastic(:) = ZERO - this % fission(:) = ZERO - this % nu_fission(:) = ZERO - this % absorption(:) = ZERO + n_temperature = size(this % kTs) + allocate(this % sum_xs(n_temperature)) + + do i = 1, n_temperature + ! Allocate and initialize derived cross sections + n_grid = size(this % grid(i) % energy) + allocate(this % sum_xs(i) % total(n_grid)) + allocate(this % sum_xs(i) % elastic(n_grid)) + allocate(this % sum_xs(i) % fission(n_grid)) + allocate(this % sum_xs(i) % nu_fission(n_grid)) + allocate(this % sum_xs(i) % absorption(n_grid)) + this % sum_xs(i) % total(:) = ZERO + this % sum_xs(i) % elastic(:) = ZERO + this % sum_xs(i) % fission(:) = ZERO + this % sum_xs(i) % nu_fission(:) = ZERO + this % sum_xs(i) % absorption(:) = ZERO + end do i_fission = 0 do i = 1, size(this % reactions) - call MTs % append(this % reactions(i) % MT) + call MTs % push_back(this % reactions(i) % MT) call this % reaction_index % add_key(this % reactions(i) % MT, i) associate (rx => this % reactions(i)) - j = rx % threshold - n = size(rx % sigma) - ! Skip total inelastic level scattering, gas production cross sections ! (MT=200+), etc. if (rx % MT == N_LEVEL .or. rx % MT == N_NONELASTIC) cycle if (rx % MT > N_5N2P .and. rx % MT < N_P0) cycle ! Skip level cross sections if total is available - if (rx % MT >= N_P0 .and. rx % MT <= N_PC .and. MTs % contains(N_P)) cycle - if (rx % MT >= N_D0 .and. rx % MT <= N_DC .and. MTs % contains(N_D)) cycle - if (rx % MT >= N_T0 .and. rx % MT <= N_TC .and. MTs % contains(N_T)) cycle - if (rx % MT >= N_3HE0 .and. rx % MT <= N_3HEC .and. MTs % contains(N_3HE)) cycle - if (rx % MT >= N_A0 .and. rx % MT <= N_AC .and. MTs % contains(N_A)) cycle - if (rx % MT >= N_2N0 .and. rx % MT <= N_2NC .and. MTs % contains(N_2N)) cycle + if (rx % MT >= N_P0 .and. rx % MT <= N_PC .and. find(MTs, N_P) /= -1) cycle + if (rx % MT >= N_D0 .and. rx % MT <= N_DC .and. find(MTs, N_D) /= -1) cycle + if (rx % MT >= N_T0 .and. rx % MT <= N_TC .and. find(MTs, N_T) /= -1) cycle + if (rx % MT >= N_3HE0 .and. rx % MT <= N_3HEC .and. find(MTs, N_3HE) /= -1) cycle + if (rx % MT >= N_A0 .and. rx % MT <= N_AC .and. find(MTs, N_A) /= -1) cycle + if (rx % MT >= N_2N0 .and. rx % MT <= N_2NC .and. find(MTs, N_2N) /= -1) cycle - ! Copy elastic - if (rx % MT == ELASTIC) this % elastic(:) = rx % sigma + do t = 1, n_temperature + j = rx % xs(t) % threshold + n = size(rx % xs(t) % value) - ! Add contribution to total cross section - this % total(j:j+n-1) = this % total(j:j+n-1) + rx % sigma + ! Copy elastic + if (rx % MT == ELASTIC) this % sum_xs(t) % elastic(:) = rx % xs(t) % value - ! Add contribution to absorption cross section - if (is_disappearance(rx % MT)) then - this % absorption(j:j+n-1) = this % absorption(j:j+n-1) + rx % sigma - end if + ! Add contribution to total cross section + this % sum_xs(t) % total(j:j+n-1) = this % sum_xs(t) % total(j:j+n-1) + & + rx % xs(t) % value - ! Information about fission reactions - if (rx % MT == N_FISSION) then - allocate(this % index_fission(1)) - elseif (rx % MT == N_F) then - allocate(this % index_fission(PARTIAL_FISSION_MAX)) - this % has_partial_fission = .true. - end if + ! Add contribution to absorption cross section + if (is_disappearance(rx % MT)) then + this % sum_xs(t) % absorption(j:j+n-1) = this % sum_xs(t) % & + absorption(j:j+n-1) + rx % xs(t) % value + end if - ! Add contribution to fission cross section - if (is_fission(rx % MT)) then - this % fissionable = .true. - this % fission(j:j+n-1) = this % fission(j:j+n-1) + rx % sigma - - ! Also need to add fission cross sections to absorption - this % absorption(j:j+n-1) = this % absorption(j:j+n-1) + rx % sigma - - ! If total fission reaction is present, there's no need to store the - ! reaction cross-section since it was copied to this % fission - if (rx % MT == N_FISSION) deallocate(rx % sigma) - - ! Keep track of this reaction for easy searching later - i_fission = i_fission + 1 - this % index_fission(i_fission) = i - this % n_fission = this % n_fission + 1 - - ! <<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<< - ! Before the secondary distribution refactor, when the angle/energy - ! distribution was uncorrelated, no angle was actually sampled. With - ! the refactor, an angle is always sampled for an uncorrelated - ! distribution even when no angle distribution exists in the ACE file - ! (isotropic is assumed). To preserve the RNG stream, we explicitly - ! mark fission reactions so that we avoid the angle sampling. - do k = 1, size(rx % products) - if (rx % products(k) % particle == NEUTRON) then - do m = 1, size(rx % products(k) % distribution) - associate (aedist => rx % products(k) % distribution(m) % obj) - select type (aedist) - type is (UncorrelatedAngleEnergy) - aedist % fission = .true. - end select - end associate - end do + ! Information about fission reactions + if (t == 1) then + if (rx % MT == N_FISSION) then + allocate(this % index_fission(1)) + elseif (rx % MT == N_F) then + allocate(this % index_fission(PARTIAL_FISSION_MAX)) + this % has_partial_fission = .true. end if - end do - ! <<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<< - end if - end associate - end do + end if + + ! Add contribution to fission cross section + if (is_fission(rx % MT)) then + this % fissionable = .true. + this % sum_xs(t) % fission(j:j+n-1) = this % sum_xs(t) % & + fission(j:j+n-1) + rx % xs(t) % value + + ! Also need to add fission cross sections to absorption + this % sum_xs(t) % absorption(j:j+n-1) = this % sum_xs(t) % & + absorption(j:j+n-1) + rx % xs(t) % value + + ! If total fission reaction is present, there's no need to store the + ! reaction cross-section since it was copied to this % fission + if (rx % MT == N_FISSION) deallocate(rx % xs(t) % value) + + ! Keep track of this reaction for easy searching later + if (t == 1) then + i_fission = i_fission + 1 + this % index_fission(i_fission) = i + this % n_fission = this % n_fission + 1 + + ! <<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<< + ! Before the secondary distribution refactor, when the angle/energy + ! distribution was uncorrelated, no angle was actually sampled. With + ! the refactor, an angle is always sampled for an uncorrelated + ! distribution even when no angle distribution exists in the ACE file + ! (isotropic is assumed). To preserve the RNG stream, we explicitly + ! mark fission reactions so that we avoid the angle sampling. + do k = 1, size(rx % products) + if (rx % products(k) % particle == NEUTRON) then + do m = 1, size(rx % products(k) % distribution) + associate (aedist => rx % products(k) % distribution(m) % obj) + select type (aedist) + type is (UncorrelatedAngleEnergy) + aedist % fission = .true. + end select + end associate + end do + end if + end do + ! <<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<< + end if + end if ! fission + end do ! temperature + end associate ! rx + end do ! reactions ! Determine number of delayed neutron precursors if (this % fissionable) then @@ -497,17 +564,16 @@ module nuclide_header end if ! Calculate nu-fission cross section - if (this % fissionable) then - do i = 1, size(this % energy) - this % nu_fission(i) = this % nu(this % energy(i), EMISSION_TOTAL) * & - this % fission(i) - end do - else - this % nu_fission(:) = ZERO - end if - - ! Clear MTs set - call MTs % clear() + do t = 1, n_temperature + if (this % fissionable) then + do i = 1, size(this % sum_xs(t) % fission) + this % sum_xs(t) % nu_fission(i) = this % nu(this % grid(t) % energy(i), & + EMISSION_TOTAL) * this % sum_xs(t) % fission(i) + end do + else + this % sum_xs(t) % nu_fission(:) = ZERO + end if + end do end subroutine nuclide_create_derived !=============================================================================== @@ -575,86 +641,4 @@ module nuclide_header end function nuclide_nu - -!=============================================================================== -! NUCLIDE*_PRINT displays information about a continuous-energy neutron -! cross_section table and its reactions and secondary angle/energy distributions -!=============================================================================== - - subroutine nuclide_print(this, unit) - class(Nuclide), intent(in) :: this - integer, intent(in), optional :: unit - - integer :: i ! loop index over nuclides - integer :: unit_ ! unit to write to - integer :: size_xs ! memory used for cross-sections (bytes) - integer :: size_urr ! memory used for probability tables (bytes) - - ! set default unit for writing information - if (present(unit)) then - unit_ = unit - else - unit_ = OUTPUT_UNIT - end if - - ! Initialize totals - size_urr = 0 - size_xs = 0 - - ! Basic nuclide information - write(unit_,*) 'Nuclide ' // trim(this % name) - write(unit_,*) ' zaid = ' // trim(to_str(this % zaid)) - write(unit_,*) ' awr = ' // trim(to_str(this % awr)) - write(unit_,*) ' kT = ' // trim(to_str(this % kT)) - write(unit_,*) ' # of grid points = ' // trim(to_str(this % n_grid)) - write(unit_,*) ' Fissionable = ', this % fissionable - write(unit_,*) ' # of fission reactions = ' // trim(to_str(this % n_fission)) - write(unit_,*) ' # of reactions = ' // trim(to_str(size(this % reactions))) - - ! Information on each reaction - write(unit_,*) ' Reaction Q-value COM IE' - do i = 1, size(this % reactions) - associate (rxn => this % reactions(i)) - write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,I6)') & - reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, & - rxn % threshold - - ! Accumulate data size - size_xs = size_xs + (this % n_grid - rxn%threshold + 1) * 8 - end associate - end do - - ! Add memory required for summary reactions (total, absorption, fission, - ! nu-fission) - size_xs = 8 * this % n_grid * 4 - - ! Write information about URR probability tables - size_urr = 0 - if (this % urr_present) then - associate(urr => this % urr_data) - write(unit_,*) ' Unresolved resonance probability table:' - write(unit_,*) ' # of energies = ' // trim(to_str(urr % n_energy)) - write(unit_,*) ' # of probabilities = ' // trim(to_str(urr % n_prob)) - write(unit_,*) ' Interpolation = ' // trim(to_str(urr % interp)) - write(unit_,*) ' Inelastic flag = ' // trim(to_str(urr % inelastic_flag)) - write(unit_,*) ' Absorption flag = ' // trim(to_str(urr % absorption_flag)) - write(unit_,*) ' Multiply by smooth? ', urr % multiply_smooth - write(unit_,*) ' Min energy = ', trim(to_str(urr % energy(1))) - write(unit_,*) ' Max energy = ', trim(to_str(urr % energy(urr % n_energy))) - - ! Calculate memory used by probability tables and add to total - size_urr = urr % n_energy * (urr % n_prob * 6 + 1) * 8 - end associate - end if - - ! Write memory used - write(unit_,*) ' Memory Requirements' - write(unit_,*) ' Cross sections = ' // trim(to_str(size_xs)) // ' bytes' - write(unit_,*) ' Probability Tables = ' // & - trim(to_str(size_urr)) // ' bytes' - - ! Blank line at end of nuclide - write(unit_,*) - end subroutine nuclide_print - end module nuclide_header diff --git a/src/output.F90 b/src/output.F90 index 9e23f11d8..f8b0c5ab0 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -317,57 +317,6 @@ contains end subroutine print_particle -!=============================================================================== -! WRITE_XS_SUMMARY writes information about each nuclide and S(a,b) table to a -! file called cross_sections.out. This file shows the list of reactions as well -! as information about their secondary angle/energy distributions, how much -! memory is consumed, thresholds, etc. -!=============================================================================== - - subroutine write_xs_summary() - - integer :: i ! loop index - integer :: unit_xs ! cross_sections.out file unit - character(MAX_FILE_LEN) :: path ! path of summary file - - ! Create filename for log file - path = trim(path_output) // "cross_sections.out" - - ! Open log file for writing - open(NEWUNIT=unit_xs, FILE=path, STATUS='replace', ACTION='write') - - if (run_CE) then - ! Write header - call header("CROSS SECTION TABLES", unit=unit_xs) - - NUCLIDE_LOOP: do i = 1, n_nuclides_total - ! Print information about nuclide - call nuclides(i) % print(unit=unit_xs) - end do NUCLIDE_LOOP - - SAB_TABLES_LOOP: do i = 1, n_sab_tables - ! Print information about S(a,b) table - call sab_tables(i) % print(unit=unit_xs) - end do SAB_TABLES_LOOP - else - ! Write header - call header("MGXS LIBRARY TABLES", unit=unit_xs) - NuclideMG_LOOP: do i = 1, n_nuclides_total - ! Print information about nuclide - call nuclides_mg(i) % obj % print(unit=unit_xs) - end do NuclideMG_LOOP - call header("MATERIAL MGXS TABLES", unit=unit_xs) - MATERIAL_LOOP: do i = 1, n_materials - ! Print information about Materials - call macro_xs(i) % obj % print(unit=unit_xs) - end do MATERIAL_LOOP - end if - - ! Close cross section summary file - close(unit_xs) - - end subroutine write_xs_summary - !=============================================================================== ! PRINT_COLUMNS displays a header listing what physical values will displayed ! below them diff --git a/src/particle_header.F90 b/src/particle_header.F90 index ee854aea3..53a503ce2 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -88,6 +88,7 @@ module particle_header ! Temperature of the current cell real(8) :: sqrtkT ! sqrt(k_Boltzmann * temperature) in MeV + real(8) :: last_sqrtKT ! last temperature ! Statistical data integer :: n_collision ! # of collisions @@ -129,6 +130,7 @@ contains this % cell_born = NONE this % material = NONE this % last_material = NONE + this % last_sqrtkT = NONE this % wgt = ONE this % last_wgt = ONE this % absorb_wgt = ZERO diff --git a/src/physics.F90 b/src/physics.F90 index 00bd2c20c..f1f673309 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -59,7 +59,7 @@ contains ! Advance URR seed stream 'N' times after energy changes if (p % E /= p % last_E) then call prn_set_stream(STREAM_URR_PTABLE) - call advance_prn_seed(n_nuc_zaid_total) + call advance_prn_seed(size(nuclides, kind=8)) call prn_set_stream(STREAM_TRACKING) endif @@ -200,6 +200,7 @@ contains integer :: i integer :: i_grid + integer :: i_temp real(8) :: f real(8) :: prob real(8) :: cutoff @@ -219,6 +220,7 @@ contains end if ! Get grid index and interpolatoin factor and sample fission cdf + i_temp = micro_xs(i_nuclide) % index_temp i_grid = micro_xs(i_nuclide) % index_grid f = micro_xs(i_nuclide) % interp_factor cutoff = prn() * micro_xs(i_nuclide) % fission @@ -229,13 +231,13 @@ contains FISSION_REACTION_LOOP: do i = 1, nuc % n_fission i_reaction = nuc % index_fission(i) - associate (rxn => nuc % reactions(i_reaction)) + associate (xs => nuc % reactions(i_reaction) % xs(i_temp)) ! if energy is below threshold for this reaction, skip it - if (i_grid < rxn % threshold) cycle + if (i_grid < xs % threshold) cycle ! add to cumulative probability - prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) & - + f*(rxn%sigma(i_grid - rxn%threshold + 2))) + prob = prob + ((ONE - f) * xs % value(i_grid - xs % threshold + 1) & + + f*(xs % value(i_grid - xs % threshold + 2))) end associate ! Create fission bank sites if fission occurs @@ -294,6 +296,7 @@ contains integer, intent(in) :: i_nuc_mat integer :: i + integer :: i_temp integer :: i_grid real(8) :: f real(8) :: prob @@ -301,6 +304,7 @@ contains real(8) :: uvw_new(3) ! outgoing uvw for iso-in-lab scattering real(8) :: uvw_old(3) ! incoming uvw for iso-in-lab scattering real(8) :: phi ! azimuthal angle for iso-in-lab scattering + real(8) :: kT ! temperature in MeV type(Nuclide), pointer :: nuc ! copy incoming direction @@ -308,6 +312,7 @@ contains ! Get pointer to nuclide and grid index/interpolation factor nuc => nuclides(i_nuclide) + i_temp = micro_xs(i_nuclide) % index_temp i_grid = micro_xs(i_nuclide) % index_grid f = micro_xs(i_nuclide) % interp_factor @@ -328,8 +333,15 @@ contains p % E, p % coord(1) % uvw, p % mu) else + ! Determine temperature +!!$ if (temperature_method == TEMPERATURE_MULTIPOLE) then +!!$ kT = p % sqrtkT**2 +!!$ else + kT = nuc % kTs(micro_xs(i_nuclide) % index_temp) +!!$ end if + ! Perform collision physics for elastic scattering - call elastic_scatter(i_nuclide, nuc % reactions(1), & + call elastic_scatter(i_nuclide, nuc % reactions(1), kT, & p % E, p % coord(1) % uvw, p % mu, p % wgt) end if @@ -352,22 +364,24 @@ contains &// trim(nuc % name)) end if - associate (rxn => nuc % reactions(i)) + associate (rx => nuc % reactions(i)) ! Skip fission reactions - if (rxn % MT == N_FISSION .or. rxn % MT == N_F .or. rxn % MT == N_NF & - .or. rxn % MT == N_2NF .or. rxn % MT == N_3NF) cycle + if (rx % MT == N_FISSION .or. rx % MT == N_F .or. rx % MT == N_NF & + .or. rx % MT == N_2NF .or. rx % MT == N_3NF) cycle ! some materials have gas production cross sections with MT > 200 that ! are duplicates. Also MT=4 is total level inelastic scattering which ! should be skipped - if (rxn % MT >= 200 .or. rxn % MT == N_LEVEL) cycle + if (rx % MT >= 200 .or. rx % MT == N_LEVEL) cycle - ! if energy is below threshold for this reaction, skip it - if (i_grid < rxn % threshold) cycle + associate (xs => rx % xs(i_temp)) + ! if energy is below threshold for this reaction, skip it + if (i_grid < xs % threshold) cycle - ! add to cumulative probability - prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) & - + f*(rxn%sigma(i_grid - rxn%threshold + 2))) + ! add to cumulative probability + prob = prob + ((ONE - f)*xs % value(i_grid - xs % threshold + 1) & + + f*(xs % value(i_grid - xs % threshold + 2))) + end associate end associate end do @@ -401,9 +415,10 @@ contains ! target. !=============================================================================== - subroutine elastic_scatter(i_nuclide, rxn, E, uvw, mu_lab, wgt) + subroutine elastic_scatter(i_nuclide, rxn, kT, E, uvw, mu_lab, wgt) integer, intent(in) :: i_nuclide type(Reaction), intent(in) :: rxn + real(8), intent(in) :: kT ! temperature in MeV real(8), intent(inout) :: E real(8), intent(inout) :: uvw(3) real(8), intent(out) :: mu_lab @@ -430,7 +445,7 @@ contains ! Sample velocity of target nucleus if (.not. micro_xs(i_nuclide) % use_ptable) then call sample_target_velocity(nuc, v_t, E, uvw, v_n, wgt, & - & micro_xs(i_nuclide) % elastic) + micro_xs(i_nuclide) % elastic, kT) else v_t = ZERO end if @@ -494,6 +509,7 @@ contains integer :: i ! incoming energy bin integer :: j ! outgoing energy bin integer :: k ! outgoing cosine bin + integer :: i_temp ! temperature index integer :: n_energy_out ! number of outgoing energy bins real(8) :: f ! interpolation factor real(8) :: r ! used for skewed sampling & continuous @@ -502,7 +518,6 @@ contains real(8) :: mu_ijk ! outgoing cosine k for E_in(i) and E_out(j) real(8) :: mu_i1jk ! outgoing cosine k for E_in(i+1) and E_out(j) real(8) :: prob ! probability for sampling Bragg edge - type(SAlphaBeta), pointer :: sab ! Following are needed only for SAB_SECONDARY_CONT scattering integer :: l ! sampled incoming E bin (is i or i + 1) real(8) :: E_i_1, E_i_J ! endpoints on outgoing grid i @@ -514,213 +529,216 @@ contains real(8) :: frac ! interpolation factor on outgoing energy real(8) :: r1 ! RNG for outgoing energy + i_temp = micro_xs(i_nuclide) % index_temp_sab + ! Get pointer to S(a,b) table - sab => sab_tables(i_sab) + associate (sab => sab_tables(i_sab) % data(i_temp)) - ! Determine whether inelastic or elastic scattering will occur - if (prn() < micro_xs(i_nuclide) % elastic_sab / & - micro_xs(i_nuclide) % elastic) then - ! elastic scattering + ! Determine whether inelastic or elastic scattering will occur + if (prn() < micro_xs(i_nuclide) % elastic_sab / & + micro_xs(i_nuclide) % elastic) then + ! elastic scattering - ! Get index and interpolation factor for elastic grid - if (E < sab % elastic_e_in(1)) then - i = 1 - f = ZERO - else - i = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, E) - f = (E - sab%elastic_e_in(i)) / & - (sab%elastic_e_in(i+1) - sab%elastic_e_in(i)) - end if - - ! Select treatment based on elastic mode - if (sab % elastic_mode == SAB_ELASTIC_DISCRETE) then - ! With this treatment, we interpolate between two discrete cosines - ! corresponding to neighboring incoming energies. This is used for - ! data derived in the incoherent approximation - - ! Sample outgoing cosine bin - k = 1 + int(prn() * sab % n_elastic_mu) - - ! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1) - mu_ijk = sab % elastic_mu(k,i) - mu_i1jk = sab % elastic_mu(k,i+1) - - ! Cosine of angle between incoming and outgoing neutron - mu = (1 - f)*mu_ijk + f*mu_i1jk - - elseif (sab % elastic_mode == SAB_ELASTIC_EXACT) then - ! This treatment is used for data derived in the coherent - ! approximation, i.e. for crystalline structures that have Bragg - ! edges. - - ! Sample a Bragg edge between 1 and i - prob = prn() * sab % elastic_P(i+1) - if (prob < sab % elastic_P(1)) then - k = 1 + ! Get index and interpolation factor for elastic grid + if (E < sab % elastic_e_in(1)) then + i = 1 + f = ZERO else - k = binary_search(sab % elastic_P(1:i+1), i+1, prob) + i = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, E) + f = (E - sab%elastic_e_in(i)) / & + (sab%elastic_e_in(i+1) - sab%elastic_e_in(i)) end if - ! Characteristic scattering cosine for this Bragg edge - mu = ONE - TWO*sab % elastic_e_in(k) / E + ! Select treatment based on elastic mode + if (sab % elastic_mode == SAB_ELASTIC_DISCRETE) then + ! With this treatment, we interpolate between two discrete cosines + ! corresponding to neighboring incoming energies. This is used for + ! data derived in the incoherent approximation - end if + ! Sample outgoing cosine bin + k = 1 + int(prn() * sab % n_elastic_mu) - ! Outgoing energy is same as incoming energy -- no need to do anything + ! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1) + mu_ijk = sab % elastic_mu(k,i) + mu_i1jk = sab % elastic_mu(k,i+1) - else - ! Perform inelastic calculations + ! Cosine of angle between incoming and outgoing neutron + mu = (1 - f)*mu_ijk + f*mu_i1jk - ! Get index and interpolation factor for inelastic grid - if (E < sab % inelastic_e_in(1)) then - i = 1 - f = ZERO - else - i = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, E) - f = (E - sab%inelastic_e_in(i)) / & - (sab%inelastic_e_in(i+1) - sab%inelastic_e_in(i)) - end if + elseif (sab % elastic_mode == SAB_ELASTIC_EXACT) then + ! This treatment is used for data derived in the coherent + ! approximation, i.e. for crystalline structures that have Bragg + ! edges. - ! Now that we have an incoming energy bin, we need to determine the - ! outgoing energy bin. This will depend on the "secondary energy - ! mode". If the mode is 0, then the outgoing energy bin is chosen from a - ! set of equally-likely bins. If the mode is 1, then the first - ! two and last two bins are skewed to have lower probabilities than the - ! other bins (0.1 for the first and last bins and 0.4 for the second and - ! second to last bins, relative to a normal bin probability of 1). - ! Finally, if the mode is 2, then a continuous distribution (with - ! accompanying PDF and CDF is utilized) - - if ((sab % secondary_mode == SAB_SECONDARY_EQUAL) .or. & - (sab % secondary_mode == SAB_SECONDARY_SKEWED)) then - if (sab % secondary_mode == SAB_SECONDARY_EQUAL) then - ! All bins equally likely - - j = 1 + int(prn() * sab % n_inelastic_e_out) - elseif (sab % secondary_mode == SAB_SECONDARY_SKEWED) then - ! Distribution skewed away from edge points - - ! Determine number of outgoing energy and angle bins - n_energy_out = sab % n_inelastic_e_out - - r = prn() * (n_energy_out - 3) - if (r > ONE) then - ! equally likely N-4 middle bins - j = int(r) + 2 - elseif (r > 0.6_8) then - ! second to last bin has relative probability of 0.4 - j = n_energy_out - 1 - elseif (r > HALF) then - ! last bin has relative probability of 0.1 - j = n_energy_out - elseif (r > 0.1_8) then - ! second bin has relative probability of 0.4 - j = 2 + ! Sample a Bragg edge between 1 and i + prob = prn() * sab % elastic_P(i+1) + if (prob < sab % elastic_P(1)) then + k = 1 else - ! first bin has relative probability of 0.1 - j = 1 + k = binary_search(sab % elastic_P(1:i+1), i+1, prob) end if + + ! Characteristic scattering cosine for this Bragg edge + mu = ONE - TWO*sab % elastic_e_in(k) / E + end if - ! Determine outgoing energy corresponding to E_in(i) and E_in(i+1) - E_ij = sab % inelastic_e_out(j,i) - E_i1j = sab % inelastic_e_out(j,i+1) - - ! Outgoing energy - E = (1 - f)*E_ij + f*E_i1j - - ! Sample outgoing cosine bin - k = 1 + int(prn() * sab % n_inelastic_mu) - - ! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1) - mu_ijk = sab % inelastic_mu(k,j,i) - mu_i1jk = sab % inelastic_mu(k,j,i+1) - - ! Cosine of angle between incoming and outgoing neutron - mu = (1 - f)*mu_ijk + f*mu_i1jk - - else if (sab % secondary_mode == SAB_SECONDARY_CONT) then - ! Continuous secondary energy - this is to be similar to - ! Law 61 interpolation on outgoing energy - - ! Sample between ith and (i+1)th bin - r = prn() - if (f > r) then - l = i + 1 - else - l = i - end if - - ! Determine endpoints on grid i - n_energy_out = sab % inelastic_data(i) % n_e_out - E_i_1 = sab % inelastic_data(i) % e_out(1) - E_i_J = sab % inelastic_data(i) % e_out(n_energy_out) - - ! Determine endpoints on grid i + 1 - n_energy_out = sab % inelastic_data(i + 1) % n_e_out - E_i1_1 = sab % inelastic_data(i + 1) % e_out(1) - E_i1_J = sab % inelastic_data(i + 1) % e_out(n_energy_out) - - E_1 = E_i_1 + f * (E_i1_1 - E_i_1) - E_J = E_i_J + f * (E_i1_J - E_i_J) - - ! Determine outgoing energy bin - ! (First reset n_energy_out to the right value) - n_energy_out = sab % inelastic_data(l) % n_e_out - r1 = prn() - c_j = sab % inelastic_data(l) % e_out_cdf(1) - do j = 1, n_energy_out - 1 - c_j1 = sab % inelastic_data(l) % e_out_cdf(j + 1) - if (r1 < c_j1) exit - c_j = c_j1 - end do - - ! check to make sure k is <= n_energy_out - 1 - j = min(j, n_energy_out - 1) - - ! Get the data to interpolate between - E_l_j = sab % inelastic_data(l) % e_out(j) - p_l_j = sab % inelastic_data(l) % e_out_pdf(j) - - ! Next part assumes linear-linear interpolation in standard - E_l_j1 = sab % inelastic_data(l) % e_out(j + 1) - p_l_j1 = sab % inelastic_data(l) % e_out_pdf(j + 1) - - ! Find secondary energy (variable E) - frac = (p_l_j1 - p_l_j) / (E_l_j1 - E_l_j) - if (frac == ZERO) then - E = E_l_j + (r1 - c_j) / p_l_j - else - E = E_l_j + (sqrt(max(ZERO, p_l_j * p_l_j + & - TWO * frac * (r1 - c_j))) - p_l_j) / frac - end if - - ! Now interpolate between incident energy bins i and i + 1 - if (l == i) then - E = E_1 + (E - E_i_1) * (E_J - E_1) / (E_i_J - E_i_1) - else - E = E_1 + (E - E_i1_1) * (E_J - E_1) / (E_i1_J - E_i1_1) - end if - - ! Find angular distribution for closest outgoing energy bin - if (r1 - c_j < c_j1 - r1) then - j = j - else - j = j + 1 - end if - - ! Sample outgoing cosine bin - k = 1 + int(prn() * sab % n_inelastic_mu) - - ! Will use mu from the randomly chosen incoming and closest outgoing - ! energy bins - mu = sab % inelastic_data(l) % mu(k, j) + ! Outgoing energy is same as incoming energy -- no need to do anything else - call fatal_error("Invalid secondary energy mode on S(a,b) table " & - &// trim(sab % name)) - end if ! (inelastic secondary energy treatment) - end if ! (elastic or inelastic) + ! Perform inelastic calculations + + ! Get index and interpolation factor for inelastic grid + if (E < sab % inelastic_e_in(1)) then + i = 1 + f = ZERO + else + i = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, E) + f = (E - sab%inelastic_e_in(i)) / & + (sab%inelastic_e_in(i+1) - sab%inelastic_e_in(i)) + end if + + ! Now that we have an incoming energy bin, we need to determine the + ! outgoing energy bin. This will depend on the "secondary energy + ! mode". If the mode is 0, then the outgoing energy bin is chosen from a + ! set of equally-likely bins. If the mode is 1, then the first + ! two and last two bins are skewed to have lower probabilities than the + ! other bins (0.1 for the first and last bins and 0.4 for the second and + ! second to last bins, relative to a normal bin probability of 1). + ! Finally, if the mode is 2, then a continuous distribution (with + ! accompanying PDF and CDF is utilized) + + if ((sab_tables(i_sab) % secondary_mode == SAB_SECONDARY_EQUAL) .or. & + (sab_tables(i_sab) % secondary_mode == SAB_SECONDARY_SKEWED)) then + if (sab_tables(i_sab) % secondary_mode == SAB_SECONDARY_EQUAL) then + ! All bins equally likely + + j = 1 + int(prn() * sab % n_inelastic_e_out) + elseif (sab_tables(i_sab) % secondary_mode == SAB_SECONDARY_SKEWED) then + ! Distribution skewed away from edge points + + ! Determine number of outgoing energy and angle bins + n_energy_out = sab % n_inelastic_e_out + + r = prn() * (n_energy_out - 3) + if (r > ONE) then + ! equally likely N-4 middle bins + j = int(r) + 2 + elseif (r > 0.6_8) then + ! second to last bin has relative probability of 0.4 + j = n_energy_out - 1 + elseif (r > HALF) then + ! last bin has relative probability of 0.1 + j = n_energy_out + elseif (r > 0.1_8) then + ! second bin has relative probability of 0.4 + j = 2 + else + ! first bin has relative probability of 0.1 + j = 1 + end if + end if + + ! Determine outgoing energy corresponding to E_in(i) and E_in(i+1) + E_ij = sab % inelastic_e_out(j,i) + E_i1j = sab % inelastic_e_out(j,i+1) + + ! Outgoing energy + E = (1 - f)*E_ij + f*E_i1j + + ! Sample outgoing cosine bin + k = 1 + int(prn() * sab % n_inelastic_mu) + + ! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1) + mu_ijk = sab % inelastic_mu(k,j,i) + mu_i1jk = sab % inelastic_mu(k,j,i+1) + + ! Cosine of angle between incoming and outgoing neutron + mu = (1 - f)*mu_ijk + f*mu_i1jk + + else if (sab_tables(i_sab) % secondary_mode == SAB_SECONDARY_CONT) then + ! Continuous secondary energy - this is to be similar to + ! Law 61 interpolation on outgoing energy + + ! Sample between ith and (i+1)th bin + r = prn() + if (f > r) then + l = i + 1 + else + l = i + end if + + ! Determine endpoints on grid i + n_energy_out = sab % inelastic_data(i) % n_e_out + E_i_1 = sab % inelastic_data(i) % e_out(1) + E_i_J = sab % inelastic_data(i) % e_out(n_energy_out) + + ! Determine endpoints on grid i + 1 + n_energy_out = sab % inelastic_data(i + 1) % n_e_out + E_i1_1 = sab % inelastic_data(i + 1) % e_out(1) + E_i1_J = sab % inelastic_data(i + 1) % e_out(n_energy_out) + + E_1 = E_i_1 + f * (E_i1_1 - E_i_1) + E_J = E_i_J + f * (E_i1_J - E_i_J) + + ! Determine outgoing energy bin + ! (First reset n_energy_out to the right value) + n_energy_out = sab % inelastic_data(l) % n_e_out + r1 = prn() + c_j = sab % inelastic_data(l) % e_out_cdf(1) + do j = 1, n_energy_out - 1 + c_j1 = sab % inelastic_data(l) % e_out_cdf(j + 1) + if (r1 < c_j1) exit + c_j = c_j1 + end do + + ! check to make sure k is <= n_energy_out - 1 + j = min(j, n_energy_out - 1) + + ! Get the data to interpolate between + E_l_j = sab % inelastic_data(l) % e_out(j) + p_l_j = sab % inelastic_data(l) % e_out_pdf(j) + + ! Next part assumes linear-linear interpolation in standard + E_l_j1 = sab % inelastic_data(l) % e_out(j + 1) + p_l_j1 = sab % inelastic_data(l) % e_out_pdf(j + 1) + + ! Find secondary energy (variable E) + frac = (p_l_j1 - p_l_j) / (E_l_j1 - E_l_j) + if (frac == ZERO) then + E = E_l_j + (r1 - c_j) / p_l_j + else + E = E_l_j + (sqrt(max(ZERO, p_l_j * p_l_j + & + TWO * frac * (r1 - c_j))) - p_l_j) / frac + end if + + ! Now interpolate between incident energy bins i and i + 1 + if (l == i) then + E = E_1 + (E - E_i_1) * (E_J - E_1) / (E_i_J - E_i_1) + else + E = E_1 + (E - E_i1_1) * (E_J - E_1) / (E_i1_J - E_i1_1) + end if + + ! Find angular distribution for closest outgoing energy bin + if (r1 - c_j < c_j1 - r1) then + j = j + else + j = j + 1 + end if + + ! Sample outgoing cosine bin + k = 1 + int(prn() * sab % n_inelastic_mu) + + ! Will use mu from the randomly chosen incoming and closest outgoing + ! energy bins + mu = sab % inelastic_data(l) % mu(k, j) + + else + call fatal_error("Invalid secondary energy mode on S(a,b) table " & + // trim(sab_tables(i_sab) % name)) + end if ! (inelastic secondary energy treatment) + end if ! (elastic or inelastic) + end associate ! Because of floating-point roundoff, it may be possible for mu to be ! outside of the range [-1,1). In these cases, we just set mu to exactly @@ -741,19 +759,19 @@ contains ! implemented here. !=============================================================================== - subroutine sample_target_velocity(nuc, v_target, E, uvw, v_neut, wgt, xs_eff) + subroutine sample_target_velocity(nuc, v_target, E, uvw, v_neut, wgt, xs_eff, kT) type(Nuclide), intent(in) :: nuc ! target nuclide at temperature T real(8), intent(out) :: v_target(3) ! target velocity - real(8), intent(in) :: v_neut(3) ! neutron velocity real(8), intent(in) :: E ! particle energy real(8), intent(in) :: uvw(3) ! direction cosines + real(8), intent(in) :: v_neut(3) ! neutron velocity real(8), intent(inout) :: wgt ! particle weight + real(8), intent(in) :: xs_eff ! effective elastic xs at temperature T + real(8), intent(in) :: kT ! equilibrium temperature of target in MeV real(8) :: awr ! target/neutron mass ratio - real(8) :: kT ! equilibrium temperature of target in MeV real(8) :: E_rel ! trial relative energy real(8) :: xs_0K ! 0K xs at E_rel - real(8) :: xs_eff ! effective elastic xs at temperature T real(8) :: wcf ! weight correction factor real(8) :: E_red ! reduced energy (same as used by Cullen in SIGMA1) real(8) :: E_low ! lowest practical relative energy @@ -782,7 +800,6 @@ contains character(80) :: sampling_scheme ! method of target velocity sampling - kT = nuc % kT awr = nuc % awr ! check if nuclide is a resonant scatterer @@ -817,12 +834,12 @@ contains case ('cxs') ! sample target velocity with the constant cross section (cxs) approx. - call sample_cxs_target_velocity(nuc, v_target, E, uvw) + call sample_cxs_target_velocity(nuc, v_target, E, uvw, kT) case ('wcm') ! sample target velocity with the constant cross section (cxs) approx. - call sample_cxs_target_velocity(nuc, v_target, E, uvw) + call sample_cxs_target_velocity(nuc, v_target, E, uvw, kT) ! adjust weight as prescribed by the weight correction method (wcm) E_rel = dot_product((v_neut - v_target), (v_neut - v_target)) @@ -874,7 +891,7 @@ contains do ! sample target velocity with the constant cross section (cxs) approx. - call sample_cxs_target_velocity(nuc, v_target, E, uvw) + call sample_cxs_target_velocity(nuc, v_target, E, uvw, kT) ! perform Doppler broadening rejection correction (dbrc) E_rel = dot_product((v_neut - v_target), (v_neut - v_target)) @@ -986,13 +1003,13 @@ contains ! can be found in FRA-TM-123. !=============================================================================== - subroutine sample_cxs_target_velocity(nuc, v_target, E, uvw) + subroutine sample_cxs_target_velocity(nuc, v_target, E, uvw, kT) type(Nuclide), intent(in) :: nuc ! target nuclide at temperature real(8), intent(out) :: v_target(3) real(8), intent(in) :: E real(8), intent(in) :: uvw(3) + real(8), intent(in) :: kT ! equilibrium temperature of target in MeV - real(8) :: kT ! equilibrium temperature of target in MeV real(8) :: awr ! target/neutron mass ratio real(8) :: alpha ! probability of sampling f2 over f1 real(8) :: mu ! cosine of angle between neutron and target vel @@ -1004,7 +1021,6 @@ contains real(8) :: beta_vt_sq ! (beta * speed of target)^2 real(8) :: vt ! speed of target - kT = nuc % kT awr = nuc % awr beta_vn = sqrt(awr * E / kT) diff --git a/src/reaction_header.F90 b/src/reaction_header.F90 index c0501d43e..a7896e9ff 100644 --- a/src/reaction_header.F90 +++ b/src/reaction_header.F90 @@ -7,6 +7,7 @@ module reaction_header use hdf5_interface, only: read_attribute, open_group, close_group, & open_dataset, read_dataset, close_dataset, get_shape use product_header, only: ReactionProduct + use stl_vector, only: VectorInt use string, only: to_str, starts_with implicit none @@ -16,12 +17,16 @@ module reaction_header ! distributions for a single reaction in a continuous-energy ACE-format table !=============================================================================== + type TemperatureXS + integer :: threshold ! Energy grid index of threshold + real(8), allocatable :: value(:) ! Cross section values + end type TemperatureXS + type Reaction integer :: MT ! ENDF MT value real(8) :: Q_value ! Reaction Q value - integer :: threshold ! Energy grid index of threshold logical :: scatter_in_cm ! scattering system in center-of-mass? - real(8), allocatable :: sigma(:) ! Cross section values + type(TemperatureXS), allocatable :: xs(:) type(ReactionProduct), allocatable :: products(:) contains procedure :: from_hdf5 => reaction_from_hdf5 @@ -29,10 +34,10 @@ module reaction_header contains - subroutine reaction_from_hdf5(this, group_id, temperature) + subroutine reaction_from_hdf5(this, group_id, temperatures) class(Reaction), intent(inout) :: this integer(HID_T), intent(in) :: group_id - character(6), intent(in) :: temperature + type(VectorInt), intent(in) :: temperatures integer :: i integer :: cm @@ -42,11 +47,12 @@ contains integer :: n_links integer :: hdf5_err integer(HID_T) :: pgroup - integer(HID_T) :: xs, xs_group + integer(HID_T) :: xs, temp_group integer(SIZE_T) :: name_len integer(HSIZE_T) :: dims(1) integer(HSIZE_T) :: j character(MAX_WORD_LEN) :: name + character(MAX_WORD_LEN) :: temp_str ! temperature dataset name, e.g. '294K' call read_attribute(this % Q_value, group_id, 'Q_value') call read_attribute(this % MT, group_id, 'mt') @@ -54,14 +60,18 @@ contains this % scatter_in_cm = (cm == 1) ! Read cross section and threshold_idx data - xs_group = open_group(group_id, temperature) - xs = open_dataset(xs_group, 'xs') - call read_attribute(this % threshold, xs, 'threshold_idx') - call get_shape(xs, dims) - allocate(this % sigma(dims(1))) - call read_dataset(this % sigma, xs) - call close_dataset(xs) - call close_group(xs_group) + allocate(this % xs(temperatures % size())) + do i = 1, temperatures % size() + temp_str = trim(to_str(temperatures % data(i))) // "K" + temp_group = open_group(group_id, temp_str) + xs = open_dataset(temp_group, 'xs') + call read_attribute(this % xs(i) % threshold, xs, 'threshold_idx') + call get_shape(xs, dims) + allocate(this % xs(i) % value(dims(1))) + call read_dataset(this % xs(i) % value, xs) + call close_dataset(xs) + call close_group(temp_group) + end do ! Determine number of products call h5gget_info_f(group_id, storage_type, n_links, max_corder, hdf5_err) diff --git a/src/relaxng/materials.rnc b/src/relaxng/materials.rnc index 4d1b217db..c5f4efd6f 100644 --- a/src/relaxng/materials.rnc +++ b/src/relaxng/materials.rnc @@ -4,7 +4,7 @@ element materials { (element name { xsd:string { maxLength="52" } } | attribute name { xsd:string { maxLength="52" } })? & - element temperature { xsd:string { maxLength = "6" } }? & + element temperature { xsd:double }? & element density { (element value { xsd:double } | attribute value { xsd:double })? & @@ -41,7 +41,5 @@ element materials { element sab { (element name { xsd:string } | attribute name { xsd:string }) }* - }+ & - - element default_temperature { xsd:string { maxLength = "6" } }? + }+ } diff --git a/src/relaxng/materials.rng b/src/relaxng/materials.rng index 20b7b86ec..3c92dc94a 100644 --- a/src/relaxng/materials.rng +++ b/src/relaxng/materials.rng @@ -1,197 +1,186 @@ - - - - + + + + + + + + + + + + - - + + + 52 + - - + + + 52 + - + + + + + + + + + + + + + + + + + + - + - 52 + 10 - + - 52 + 10 - - - - - 6 - - - - + + + + + + + + + + + + - - + + + data + iso-in-lab + - - + + + data + iso-in-lab + - - - 10 - - - - - 10 - - + + + + + + + + + + + + + + + + - - - - - - - - - - - - - - - - data - iso-in-lab - - - - - data - iso-in-lab - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + - + + 2 + - + + 2 + - - - - - + - - - 2 - + + + data + iso-in-lab + - - - 2 - + + + data + iso-in-lab + - - - - - data - iso-in-lab - - - - - data - iso-in-lab - - - - - - - - - - - - - - - - - - - - - - - - - - - + - - - - - - + + + + + + + + + + + + + + + + - - - - - - - - - 6 - - - - + + + + + + + + + + + + + + + + + + diff --git a/src/relaxng/settings.rnc b/src/relaxng/settings.rnc index c554dfa75..70550a40f 100644 --- a/src/relaxng/settings.rnc +++ b/src/relaxng/settings.rnc @@ -128,6 +128,12 @@ element settings { element survival_biasing { xsd:boolean }? & + element temperature_default { xsd:double }? & + + element temperature_method { xsd:string }? & + + element temperature_tolerance { xsd:double }? & + element threads { xsd:positiveInteger }? & element trace { list { xsd:positiveInteger+ } }? & diff --git a/src/relaxng/settings.rng b/src/relaxng/settings.rng index 0719bcbc6..246c78e68 100644 --- a/src/relaxng/settings.rng +++ b/src/relaxng/settings.rng @@ -565,6 +565,21 @@ + + + + + + + + + + + + + + + diff --git a/src/sab_header.F90 b/src/sab_header.F90 index 5facd3fc3..c8619efa3 100644 --- a/src/sab_header.F90 +++ b/src/sab_header.F90 @@ -2,14 +2,17 @@ module sab_header use, intrinsic :: ISO_FORTRAN_ENV + use algorithm, only: find, sort use constants + use dict_header, only: DictIntInt use distribution_univariate, only: Tabular - use error, only: warning + use error, only: warning, fatal_error use hdf5, only: HID_T, HSIZE_T, SIZE_T use h5lt, only: h5ltpath_valid_f, h5iget_name_f use hdf5_interface, only: read_attribute, get_shape, open_group, close_group, & open_dataset, read_dataset, close_dataset, get_datasets use secondary_correlated, only: CorrelatedAngleEnergy + use stl_vector, only: VectorInt, VectorReal use string, only: to_str, str_to_int implicit none @@ -33,13 +36,7 @@ module sab_header ! of light isotopes such as water, graphite, Be, etc !=============================================================================== - type SAlphaBeta - character(100) :: name ! name of table, e.g. lwtr.10t - real(8) :: awr ! weight of nucleus in neutron masses - real(8) :: kT ! temperature in MeV (k*T) - integer :: n_zaid ! Number of valid zaids - integer, allocatable :: zaid(:) ! List of valid Z and A identifiers, e.g. 6012 - + type SabData ! threshold for S(a,b) treatment (usually ~4 eV) real(8) :: threshold_inelastic real(8) :: threshold_elastic = ZERO @@ -48,7 +45,6 @@ module sab_header integer :: n_inelastic_e_in ! # of incoming E for inelastic integer :: n_inelastic_e_out ! # of outgoing E for inelastic integer :: n_inelastic_mu ! # of outgoing angles for inelastic - integer :: secondary_mode ! secondary mode (equal/skewed/continuous) real(8), allocatable :: inelastic_e_in(:) real(8), allocatable :: inelastic_sigma(:) ! The following are used only if secondary_mode is 0 or 1 @@ -67,101 +63,35 @@ module sab_header real(8), allocatable :: elastic_e_in(:) real(8), allocatable :: elastic_P(:) real(8), allocatable :: elastic_mu(:,:) + end type SabData + + type SAlphaBeta + character(100) :: name ! name of table, e.g. lwtr.10t + real(8) :: awr ! weight of nucleus in neutron masses + real(8), allocatable :: kTs(:) ! temperatures in MeV (k*T) + character(10), allocatable :: nuclides(:) ! List of valid nuclides + integer :: secondary_mode ! secondary mode (equal/skewed/continuous) + + ! cross sections and distributions at each temperature + type(SabData), allocatable :: data(:) contains - procedure :: print => salphabeta_print procedure :: from_hdf5 => salphabeta_from_hdf5 end type SAlphaBeta contains -!=============================================================================== -! PRINT_SAB_TABLE displays information about a S(a,b) table containing data -! describing thermal scattering from bound materials such as hydrogen in water. -!=============================================================================== - - subroutine salphabeta_print(this, unit) - class(SAlphaBeta), intent(in) :: this - integer, intent(in), optional :: unit - - integer :: size_sab ! memory used by S(a,b) table - integer :: unit_ ! unit to write to - integer :: i ! Loop counter for parsing through this % zaid - integer :: char_count ! Counter for the number of characters on a line - - ! set default unit for writing information - if (present(unit)) then - unit_ = unit - else - unit_ = OUTPUT_UNIT - end if - - ! Basic S(a,b) table information - write(unit_,*) 'S(a,b) Table ' // trim(this % name) - write(unit_,'(A)',advance="no") ' zaids = ' - ! Initialize the counter based on the above string - char_count = 11 - do i = 1, this % n_zaid - ! Deal with a line thats too long - if (char_count >= 73) then ! 73 = 80 - (5 ZAID chars + 1 space + 1 comma) - ! End the line - write(unit_,*) "" - ! Add 11 leading blanks - write(unit_,'(A)', advance="no") " " - ! reset the counter to 11 - char_count = 11 - end if - if (i < this % n_zaid) then - ! Include a comma - write(unit_,'(A)',advance="no") trim(to_str(this % zaid(i))) // ", " - char_count = char_count + len(trim(to_str(this % zaid(i)))) + 2 - else - ! Don't include a comma, since we are all done - write(unit_,'(A)',advance="no") trim(to_str(this % zaid(i))) - end if - - end do - write(unit_,*) "" ! Move to next line - write(unit_,*) ' awr = ' // trim(to_str(this % awr)) - write(unit_,*) ' kT = ' // trim(to_str(this % kT)) - - ! Inelastic data - write(unit_,*) ' # of Incoming Energies (Inelastic) = ' // & - trim(to_str(this % n_inelastic_e_in)) - write(unit_,*) ' # of Outgoing Energies (Inelastic) = ' // & - trim(to_str(this % n_inelastic_e_out)) - write(unit_,*) ' # of Outgoing Angles (Inelastic) = ' // & - trim(to_str(this % n_inelastic_mu)) - write(unit_,*) ' Threshold for Inelastic = ' // & - trim(to_str(this % threshold_inelastic)) - - ! Elastic data - if (this % n_elastic_e_in > 0) then - write(unit_,*) ' # of Incoming Energies (Elastic) = ' // & - trim(to_str(this % n_elastic_e_in)) - write(unit_,*) ' # of Outgoing Angles (Elastic) = ' // & - trim(to_str(this % n_elastic_mu)) - write(unit_,*) ' Threshold for Elastic = ' // & - trim(to_str(this % threshold_elastic)) - end if - - ! Determine memory used by S(a,b) table and write out - size_sab = 8 * (this % n_inelastic_e_in * (2 + this % n_inelastic_e_out * & - (1 + this % n_inelastic_mu)) + this % n_elastic_e_in * & - (2 + this % n_elastic_mu)) - write(unit_,*) ' Memory Used = ' // trim(to_str(size_sab)) // ' bytes' - - ! Blank line at end - write(unit_,*) - - end subroutine salphabeta_print - - subroutine salphabeta_from_hdf5(this, group_id, temperature) + subroutine salphabeta_from_hdf5(this, group_id, temperature, tolerance) class(SAlphaBeta), intent(inout) :: this integer(HID_T), intent(in) :: group_id - character(6), intent(in) :: temperature + type(VectorReal), intent(in) :: temperature ! list of temperatures + real(8), intent(in) :: tolerance integer :: i, j + integer :: n + integer :: t integer :: n_energy, n_energy_out, n_mu + integer :: i_closest + integer :: n_temperature integer :: hdf5_err integer(SIZE_T) :: name_len, name_file_len integer(HID_T) :: T_group @@ -175,10 +105,13 @@ contains character(20) :: type logical :: exists type(CorrelatedAngleEnergy) :: correlated_dist - character(MAX_FILE_LEN), allocatable :: temperatures(:) - integer, allocatable :: temperatures_integer(:) - character(6) :: my_temperature - integer :: temperature_integer + + character(MAX_WORD_LEN) :: temp_str + character(MAX_FILE_LEN), allocatable :: dset_names(:) + real(8), allocatable :: temps_available(:) ! temperatures available + real(8) :: temp_desired + real(8) :: temp_actual + type(VectorInt) :: temps_to_read ! Get name of table from group name_len = len(this % name) @@ -188,7 +121,7 @@ contains this % name = trim(this % name(2:)) call read_attribute(this % awr, group_id, 'atomic_weight_ratio') - call read_attribute(this % zaid, group_id, 'zaids') + call read_attribute(this % nuclides, group_id, 'nuclides') call read_attribute(type, group_id, 'secondary_mode') select case (type) case ('equal') @@ -198,168 +131,185 @@ contains case ('continuous') this % secondary_mode = SAB_SECONDARY_CONT end select - this % n_zaid = size(this % zaid) + + ! Read temperatures kT_group = open_group(group_id, 'kTs') - ! Before accessing the temperature data, see if the user-provied temperature - ! exists. We can find this out by looking at the datasets within kT_group - temperature_integer = & - str_to_int(temperature(1: len_trim(adjustl(temperature)) - 1)) - call get_datasets(kT_group, temperatures) - allocate(temperatures_integer(size(temperatures))) - do i = 1, size(temperatures) - temperatures_integer(i) = & - str_to_int(temperatures(i)(1: len_trim(adjustl(temperatures(i))) - 1)) + ! Determine temperatures available + call get_datasets(kT_group, dset_names) + allocate(temps_available(size(dset_names))) + do i = 1, size(dset_names) + ! Read temperature value + call read_dataset(temps_available(i), kT_group, trim(dset_names(i))) + temps_available(i) = temps_available(i) / K_BOLTZMANN end do - my_temperature = & - temperatures(minloc(abs(temperature_integer - temperatures_integer), & - dim=1)) - ! Now print a warning if there is no matching temperature and then use the - ! closest temperature - if (temperature /= my_temperature) then - if (temperature == '0K') then - call warning(trim(this % name) // " does not contain 0K data & - &needed for the resonance scattering options selected") + ! Determine actual temperatures to read + TEMP_LOOP: do i = 1, temperature % size() + temp_desired = temperature % data(i) + i_closest = minloc(abs(temps_available - temp_desired), dim=1) + temp_actual = temps_available(i_closest) + if (abs(temp_actual - temp_desired) < tolerance) then + if (find(temps_to_read, nint(temp_actual)) == -1) then + call temps_to_read % push_back(nint(temp_actual)) + end if else - call warning(trim(this % name) // " does not contain data at a & - &temperature of " // trim(temperature) // "; using the & - &nearest available temperature of " // trim(my_temperature)) + call fatal_error("Nuclear data library does not contain cross sections & + &for " // trim(this % name) // " at or near " // & + trim(to_str(nint(temp_desired))) // " K.") end if - end if + end do TEMP_LOOP - kT_dset = open_dataset(kT_group, my_temperature) - call read_dataset(this % kT, kT_dset) - call close_dataset(kT_dset) - call close_group(kT_group) + ! TODO: If using interpolation, add a block to add bounding temperatures for + ! each - ! Open my_temperature group - T_group = open_group(group_id, my_temperature) + ! Sort temperatures to read + call sort(temps_to_read) - ! Coherent elastic data - call h5ltpath_valid_f(T_group, 'elastic', .true., exists, hdf5_err) - if (exists) then - ! Read cross section data - elastic_group = open_group(T_group, 'elastic') - dset_id = open_dataset(elastic_group, 'xs') - call read_attribute(type, dset_id, 'type') - call get_shape(dset_id, dims2) - allocate(temp(dims2(1), dims2(2))) - call read_dataset(temp, dset_id) - call close_dataset(dset_id) + n_temperature = temps_to_read % size() + allocate(this % kTs(n_temperature)) + allocate(this % data(n_temperature)) - ! Set cross section data and type - this % n_elastic_e_in = int(dims2(1), 4) - allocate(this % elastic_e_in(this % n_elastic_e_in)) - allocate(this % elastic_P(this % n_elastic_e_in)) - this % elastic_e_in(:) = temp(:, 1) - this % elastic_P(:) = temp(:, 2) - select case (type) - case ('tab1') - this % elastic_mode = SAB_ELASTIC_DISCRETE - case ('bragg') - this % elastic_mode = SAB_ELASTIC_EXACT - end select - deallocate(temp) + do t = 1, n_temperature + ! Get temperature as a string + temp_str = trim(to_str(temps_to_read % data(t))) // "K" - ! Set elastic threshold - this % threshold_elastic = this % elastic_e_in(this % n_elastic_e_in) + ! Read exact temperature value + call read_dataset(this % kTs(t), kT_group, temp_str) - ! Read angle distribution - if (this % elastic_mode /= SAB_ELASTIC_EXACT) then - dset_id = open_dataset(elastic_group, 'mu_out') + ! Open group for temperature i + T_group = open_group(group_id, temp_str) + + ! Coherent elastic data + call h5ltpath_valid_f(T_group, 'elastic', .true., exists, hdf5_err) + if (exists) then + ! Read cross section data + elastic_group = open_group(T_group, 'elastic') + dset_id = open_dataset(elastic_group, 'xs') + call read_attribute(type, dset_id, 'type') call get_shape(dset_id, dims2) - this % n_elastic_mu = int(dims2(1), 4) - allocate(this % elastic_mu(dims2(1), dims2(2))) - call read_dataset(this % elastic_mu, dset_id) + allocate(temp(dims2(1), dims2(2))) + call read_dataset(temp, dset_id) call close_dataset(dset_id) - end if - call close_group(elastic_group) - end if + ! Set cross section data and type + this % data(t) % n_elastic_e_in = int(dims2(1), 4) + allocate(this % data(t) % elastic_e_in(this % data(t) % n_elastic_e_in)) + allocate(this % data(t) % elastic_P(this % data(t) % n_elastic_e_in)) + this % data(t) % elastic_e_in(:) = temp(:, 1) + this % data(t) % elastic_P(:) = temp(:, 2) + select case (type) + case ('tab1') + this % data(t) % elastic_mode = SAB_ELASTIC_DISCRETE + case ('bragg') + this % data(t) % elastic_mode = SAB_ELASTIC_EXACT + end select + deallocate(temp) - ! Inelastic data - call h5ltpath_valid_f(T_group, 'inelastic', .true., exists, hdf5_err) - if (exists) then - ! Read type of inelastic data - inelastic_group = open_group(T_group, 'inelastic') - - ! Read cross section data - dset_id = open_dataset(inelastic_group, 'xs') - call get_shape(dset_id, dims2) - allocate(temp(dims2(1), dims2(2))) - call read_dataset(temp, dset_id) - call close_dataset(dset_id) - - ! Set cross section data - this % n_inelastic_e_in = int(dims2(1), 4) - allocate(this % inelastic_e_in(this % n_inelastic_e_in)) - allocate(this % inelastic_sigma(this % n_inelastic_e_in)) - this % inelastic_e_in(:) = temp(:, 1) - this % inelastic_sigma(:) = temp(:, 2) - deallocate(temp) - - ! Set inelastic threshold - this % threshold_inelastic = this % inelastic_e_in(this % n_inelastic_e_in) - - if (this % secondary_mode /= SAB_SECONDARY_CONT) then - ! Read energy distribution - dset_id = open_dataset(inelastic_group, 'energy_out') - call get_shape(dset_id, dims2) - this % n_inelastic_e_out = int(dims2(1), 4) - allocate(this % inelastic_e_out(dims2(1), dims2(2))) - call read_dataset(this % inelastic_e_out, dset_id) - call close_dataset(dset_id) + ! Set elastic threshold + this % data(t) % threshold_elastic = this % data(t) % elastic_e_in(& + this % data(t) % n_elastic_e_in) ! Read angle distribution - dset_id = open_dataset(inelastic_group, 'mu_out') - call get_shape(dset_id, dims3) - this % n_inelastic_mu = int(dims3(1), 4) - allocate(this % inelastic_mu(dims3(1), dims3(2), dims3(3))) - call read_dataset(this % inelastic_mu, dset_id) - call close_dataset(dset_id) - else - ! Read correlated angle-energy distribution - call correlated_dist % from_hdf5(inelastic_group) + if (this % data(t) % elastic_mode /= SAB_ELASTIC_EXACT) then + dset_id = open_dataset(elastic_group, 'mu_out') + call get_shape(dset_id, dims2) + this % data(t) % n_elastic_mu = int(dims2(1), 4) + allocate(this % data(t) % elastic_mu(dims2(1), dims2(2))) + call read_dataset(this % data(t) % elastic_mu, dset_id) + call close_dataset(dset_id) + end if - ! Convert to S(a,b) native format - n_energy = size(correlated_dist % energy) - allocate(this % inelastic_data(n_energy)) - do i = 1, n_energy - associate (edist => correlated_dist % distribution(i)) - ! Get number of outgoing energies for incoming energy i - n_energy_out = size(edist % e_out) - this % inelastic_data(i) % n_e_out = n_energy_out - allocate(this % inelastic_data(i) % e_out(n_energy_out)) - allocate(this % inelastic_data(i) % e_out_pdf(n_energy_out)) - allocate(this % inelastic_data(i) % e_out_cdf(n_energy_out)) - - ! Copy outgoing energy distribution - this % inelastic_data(i) % e_out(:) = edist % e_out - this % inelastic_data(i) % e_out_pdf(:) = edist % p - this % inelastic_data(i) % e_out_cdf(:) = edist % c - - do j = 1, n_energy_out - select type (adist => edist % angle(j) % obj) - type is (Tabular) - ! On first pass, allocate space for angles - if (j == 1) then - n_mu = size(adist % x) - this % n_inelastic_mu = n_mu - allocate(this % inelastic_data(i) % mu(n_mu, n_energy_out)) - end if - - ! Copy outgoing angles - this % inelastic_data(i) % mu(:, j) = adist % x - end select - end do - end associate - end do + call close_group(elastic_group) end if - call close_group(inelastic_group) - end if - call close_group(T_group) + ! Inelastic data + call h5ltpath_valid_f(T_group, 'inelastic', .true., exists, hdf5_err) + if (exists) then + ! Read type of inelastic data + inelastic_group = open_group(T_group, 'inelastic') + + ! Read cross section data + dset_id = open_dataset(inelastic_group, 'xs') + call get_shape(dset_id, dims2) + allocate(temp(dims2(1), dims2(2))) + call read_dataset(temp, dset_id) + call close_dataset(dset_id) + + ! Set cross section data + this % data(t) % n_inelastic_e_in = int(dims2(1), 4) + allocate(this % data(t) % inelastic_e_in(this % data(t) % n_inelastic_e_in)) + allocate(this % data(t) % inelastic_sigma(this % data(t) % n_inelastic_e_in)) + this % data(t) % inelastic_e_in(:) = temp(:, 1) + this % data(t) % inelastic_sigma(:) = temp(:, 2) + deallocate(temp) + + ! Set inelastic threshold + this % data(t) % threshold_inelastic = this % data(t) % inelastic_e_in(& + this % data(t) % n_inelastic_e_in) + + if (this % secondary_mode /= SAB_SECONDARY_CONT) then + ! Read energy distribution + dset_id = open_dataset(inelastic_group, 'energy_out') + call get_shape(dset_id, dims2) + this % data(t) % n_inelastic_e_out = int(dims2(1), 4) + allocate(this % data(t) % inelastic_e_out(dims2(1), dims2(2))) + call read_dataset(this % data(t) % inelastic_e_out, dset_id) + call close_dataset(dset_id) + + ! Read angle distribution + dset_id = open_dataset(inelastic_group, 'mu_out') + call get_shape(dset_id, dims3) + this % data(t) % n_inelastic_mu = int(dims3(1), 4) + allocate(this % data(t) % inelastic_mu(dims3(1), dims3(2), dims3(3))) + call read_dataset(this % data(t) % inelastic_mu, dset_id) + call close_dataset(dset_id) + else + ! Read correlated angle-energy distribution + call correlated_dist % from_hdf5(inelastic_group) + + ! Convert to S(a,b) native format + n_energy = size(correlated_dist % energy) + allocate(this % data(t) % inelastic_data(n_energy)) + do i = 1, n_energy + associate (edist => correlated_dist % distribution(i)) + ! Get number of outgoing energies for incoming energy i + n_energy_out = size(edist % e_out) + this % data(t) % inelastic_data(i) % n_e_out = n_energy_out + allocate(this % data(t) % inelastic_data(i) % e_out(n_energy_out)) + allocate(this % data(t) % inelastic_data(i) % e_out_pdf(n_energy_out)) + allocate(this % data(t) % inelastic_data(i) % e_out_cdf(n_energy_out)) + + ! Copy outgoing energy distribution + this % data(t) % inelastic_data(i) % e_out(:) = edist % e_out + this % data(t) % inelastic_data(i) % e_out_pdf(:) = edist % p + this % data(t) % inelastic_data(i) % e_out_cdf(:) = edist % c + + do j = 1, n_energy_out + select type (adist => edist % angle(j) % obj) + type is (Tabular) + ! On first pass, allocate space for angles + if (j == 1) then + n_mu = size(adist % x) + this % data(t) % n_inelastic_mu = n_mu + allocate(this % data(t) % inelastic_data(i) % mu(& + n_mu, n_energy_out)) + end if + + ! Copy outgoing angles + this % data(t) % inelastic_data(i) % mu(:, j) = adist % x + end select + end do + end associate + end do + end if + + call close_group(inelastic_group) + end if + call close_group(T_group) + end do + + call close_group(kT_group) end subroutine salphabeta_from_hdf5 end module sab_header diff --git a/src/summary.F90 b/src/summary.F90 index 5bd440ff3..cc0c51758 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -115,36 +115,31 @@ contains integer :: i character(12), allocatable :: nucnames(:) real(8), allocatable :: awrs(:) - integer, allocatable :: zaids(:) ! Write useful data from nuclide objects nuclide_group = create_group(file_id, "nuclides") call write_dataset(nuclide_group, "n_nuclides_total", n_nuclides_total) - ! Build array of nuclide names, awrs, and zaids + ! Build array of nuclide names and awrs allocate(nucnames(n_nuclides_total)) allocate(awrs(n_nuclides_total)) - allocate(zaids(n_nuclides_total)) do i = 1, n_nuclides_total if (run_CE) then nucnames(i) = nuclides(i) % name awrs(i) = nuclides(i) % awr - zaids(i) = nuclides(i) % zaid else nucnames(i) = nuclides_MG(i) % obj % name awrs(i) = nuclides_MG(i) % obj % awr - zaids(i) = nuclides_MG(i) % obj % zaid end if end do - ! Write nuclide names, awrs and zaids + ! Write nuclide names and awrs call write_dataset(nuclide_group, "names", nucnames) call write_dataset(nuclide_group, "awrs", awrs) - call write_dataset(nuclide_group, "zaids", zaids) call close_group(nuclide_group) - deallocate(nucnames, awrs, zaids) + deallocate(nucnames, awrs) end subroutine write_nuclides @@ -529,9 +524,6 @@ contains ! Write name for this material call write_dataset(material_group, "name", m % name) - ! Write temperature for this material - call write_dataset(material_group, "temperature", m % temperature) - ! Write atom density with units call write_dataset(material_group, "atom_density", m % density) call write_attribute_string(material_group, "atom_density", "units", & diff --git a/src/tally.F90 b/src/tally.F90 index 3ee4702e2..3cee65b96 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -88,6 +88,7 @@ contains integer :: l ! loop index for nuclides in material integer :: m ! loop index for reactions integer :: q ! loop index for scoring bins + integer :: i_temp ! temperature index integer :: i_nuc ! index in nuclides array (from material) integer :: i_energy ! index in nuclide energy grid integer :: score_bin ! scoring bin, e.g. SCORE_FLUX @@ -887,16 +888,18 @@ contains if (i_nuclide > 0) then if (nuclides(i_nuclide)%reaction_index%has_key(score_bin)) then m = nuclides(i_nuclide)%reaction_index%get_key(score_bin) - associate (rxn => nuclides(i_nuclide) % reactions(m)) - ! Retrieve index on nuclide energy grid and interpolation - ! factor - i_energy = micro_xs(i_nuclide) % index_grid - f = micro_xs(i_nuclide) % interp_factor - if (i_energy >= rxn % threshold) then - score = ((ONE - f) * rxn % sigma(i_energy - & - rxn%threshold + 1) + f * rxn % sigma(i_energy - & - rxn%threshold + 2)) * atom_density * flux + ! Retrieve temperature and energy grid index and interpolation + ! factor + i_temp = micro_xs(i_nuclide) % index_temp + i_energy = micro_xs(i_nuclide) % index_grid + f = micro_xs(i_nuclide) % interp_factor + + associate (xs => nuclides(i_nuclide) % reactions(m) % xs(i_temp)) + if (i_energy >= xs % threshold) then + score = ((ONE - f) * xs % value(i_energy - & + xs % threshold + 1) + f * xs % value(i_energy - & + xs % threshold + 2)) * atom_density * flux end if end associate end if @@ -911,15 +914,18 @@ contains if (nuclides(i_nuc)%reaction_index%has_key(score_bin)) then m = nuclides(i_nuc)%reaction_index%get_key(score_bin) - associate (rxn => nuclides(i_nuc) % reactions(m)) - ! Retrieve index on nuclide energy grid and interpolation - ! factor - i_energy = micro_xs(i_nuc) % index_grid - f = micro_xs(i_nuc) % interp_factor - if (i_energy >= rxn % threshold) then - score = score + ((ONE - f) * rxn % sigma(i_energy - & - rxn%threshold + 1) + f * rxn % sigma(i_energy - & - rxn%threshold + 2)) * atom_density_ * flux + + ! Retrieve temperature and energy grid index and interpolation + ! factor + i_temp = micro_xs(i_nuc) % index_temp + i_energy = micro_xs(i_nuc) % index_grid + f = micro_xs(i_nuc) % interp_factor + + associate (xs => nuclides(i_nuc) % reactions(m) % xs(i_temp)) + if (i_energy >= xs % threshold) then + score = score + ((ONE - f) * xs % value(i_energy - & + xs % threshold + 1) + f * xs % value(i_energy - & + xs % threshold + 2)) * atom_density_ * flux end if end associate end if diff --git a/src/tracking.F90 b/src/tracking.F90 index 69fb78c35..f2613146e 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -84,9 +84,10 @@ contains ! Calculate microscopic and macroscopic cross sections if (run_CE) then - ! If the material is the same as the last material and the energy of the - ! particle hasn't changed, we don't need to lookup cross sections again. - if (p % material /= p % last_material) call calculate_xs(p) + ! If the material is the same as the last material and the temperature + ! hasn't changed, we don't need to lookup cross sections again. + if (p % material /= p % last_material .or. & + p % sqrtkT /= p % last_sqrtkT) call calculate_xs(p) else ! Since the MGXS can be angle dependent, this needs to be done ! After every collision for the MGXS mode diff --git a/tests/test_cmfd_feed/materials.xml b/tests/test_cmfd_feed/materials.xml index f773d83d0..70580e3a8 100644 --- a/tests/test_cmfd_feed/materials.xml +++ b/tests/test_cmfd_feed/materials.xml @@ -1,9 +1,7 @@ -294K - - + diff --git a/tests/test_cmfd_nofeed/materials.xml b/tests/test_cmfd_nofeed/materials.xml index f773d83d0..70580e3a8 100644 --- a/tests/test_cmfd_nofeed/materials.xml +++ b/tests/test_cmfd_nofeed/materials.xml @@ -1,9 +1,7 @@ -294K - - + diff --git a/tests/test_complex_cell/materials.xml b/tests/test_complex_cell/materials.xml index 857626d5a..6edf0a5f9 100644 --- a/tests/test_complex_cell/materials.xml +++ b/tests/test_complex_cell/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_confidence_intervals/materials.xml b/tests/test_confidence_intervals/materials.xml index e1124a303..0965c8783 100644 --- a/tests/test_confidence_intervals/materials.xml +++ b/tests/test_confidence_intervals/materials.xml @@ -2,7 +2,7 @@ - 294K + 294 diff --git a/tests/test_density/materials.xml b/tests/test_density/materials.xml index 53a7b0fcb..7b49233ed 100644 --- a/tests/test_density/materials.xml +++ b/tests/test_density/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_eigenvalue_genperbatch/materials.xml b/tests/test_eigenvalue_genperbatch/materials.xml index 37cfcdc6e..2472a7471 100644 --- a/tests/test_eigenvalue_genperbatch/materials.xml +++ b/tests/test_eigenvalue_genperbatch/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_eigenvalue_no_inactive/materials.xml b/tests/test_eigenvalue_no_inactive/materials.xml index 37cfcdc6e..2472a7471 100644 --- a/tests/test_eigenvalue_no_inactive/materials.xml +++ b/tests/test_eigenvalue_no_inactive/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_energy_grid/materials.xml b/tests/test_energy_grid/materials.xml index 2946a8d60..87a3b8fe1 100644 --- a/tests/test_energy_grid/materials.xml +++ b/tests/test_energy_grid/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_energy_grid/settings.xml b/tests/test_energy_grid/settings.xml index f925356a9..1e4b5937b 100644 --- a/tests/test_energy_grid/settings.xml +++ b/tests/test_energy_grid/settings.xml @@ -1,7 +1,7 @@ - nuclide + 20000 10 diff --git a/tests/test_energy_laws/materials.xml b/tests/test_energy_laws/materials.xml index e7acd32e2..e63f4018c 100644 --- a/tests/test_energy_laws/materials.xml +++ b/tests/test_energy_laws/materials.xml @@ -1,6 +1,5 @@ - 294K diff --git a/tests/test_entropy/materials.xml b/tests/test_entropy/materials.xml index 37cfcdc6e..2472a7471 100644 --- a/tests/test_entropy/materials.xml +++ b/tests/test_entropy/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_filter_distribcell/case-1/materials.xml b/tests/test_filter_distribcell/case-1/materials.xml index afd4b05b7..e7108b477 100644 --- a/tests/test_filter_distribcell/case-1/materials.xml +++ b/tests/test_filter_distribcell/case-1/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_filter_distribcell/case-2/materials.xml b/tests/test_filter_distribcell/case-2/materials.xml index afd4b05b7..794d410a3 100644 --- a/tests/test_filter_distribcell/case-2/materials.xml +++ b/tests/test_filter_distribcell/case-2/materials.xml @@ -1,9 +1,6 @@ - 294K - - diff --git a/tests/test_filter_distribcell/case-3/materials.xml b/tests/test_filter_distribcell/case-3/materials.xml index f03693e52..588912271 100644 --- a/tests/test_filter_distribcell/case-3/materials.xml +++ b/tests/test_filter_distribcell/case-3/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_filter_distribcell/case-4/materials.xml b/tests/test_filter_distribcell/case-4/materials.xml index 15cf0c54b..2eb744fe6 100644 --- a/tests/test_filter_distribcell/case-4/materials.xml +++ b/tests/test_filter_distribcell/case-4/materials.xml @@ -1,7 +1,5 @@ - 294K - diff --git a/tests/test_filter_mesh_2d/materials.xml b/tests/test_filter_mesh_2d/materials.xml index bead58f8b..8021f5f99 100644 --- a/tests/test_filter_mesh_2d/materials.xml +++ b/tests/test_filter_mesh_2d/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_filter_mesh_3d/materials.xml b/tests/test_filter_mesh_3d/materials.xml index bead58f8b..8021f5f99 100644 --- a/tests/test_filter_mesh_3d/materials.xml +++ b/tests/test_filter_mesh_3d/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_fixed_source/materials.xml b/tests/test_fixed_source/materials.xml index cfd0efd6c..6e4249da3 100644 --- a/tests/test_fixed_source/materials.xml +++ b/tests/test_fixed_source/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_fixed_source/settings.xml b/tests/test_fixed_source/settings.xml index 8a0ddb251..1e9b85d5a 100644 --- a/tests/test_fixed_source/settings.xml +++ b/tests/test_fixed_source/settings.xml @@ -6,6 +6,8 @@ 100 + 294 + diff --git a/tests/test_infinite_cell/materials.xml b/tests/test_infinite_cell/materials.xml index 8bd6c820b..6acd8df74 100644 --- a/tests/test_infinite_cell/materials.xml +++ b/tests/test_infinite_cell/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_lattice/materials.xml b/tests/test_lattice/materials.xml index 91772da79..971f5c548 100644 --- a/tests/test_lattice/materials.xml +++ b/tests/test_lattice/materials.xml @@ -10,8 +10,6 @@ =============================================================== --> - 294K - diff --git a/tests/test_lattice_hex/materials.xml b/tests/test_lattice_hex/materials.xml index 90b79b7b4..c7649fcf9 100644 --- a/tests/test_lattice_hex/materials.xml +++ b/tests/test_lattice_hex/materials.xml @@ -1,9 +1,7 @@ - 294K - - + @@ -11,7 +9,7 @@ - + @@ -19,7 +17,7 @@ - + @@ -29,7 +27,7 @@ - + diff --git a/tests/test_lattice_mixed/materials.xml b/tests/test_lattice_mixed/materials.xml index 90b79b7b4..c7649fcf9 100644 --- a/tests/test_lattice_mixed/materials.xml +++ b/tests/test_lattice_mixed/materials.xml @@ -1,9 +1,7 @@ - 294K - - + @@ -11,7 +9,7 @@ - + @@ -19,7 +17,7 @@ - + @@ -29,7 +27,7 @@ - + diff --git a/tests/test_lattice_multiple/materials.xml b/tests/test_lattice_multiple/materials.xml index bead58f8b..8021f5f99 100644 --- a/tests/test_lattice_multiple/materials.xml +++ b/tests/test_lattice_multiple/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_multipole/inputs_true.dat b/tests/test_multipole/inputs_true.dat index 5498e6873..930be9536 100644 --- a/tests/test_multipole/inputs_true.dat +++ b/tests/test_multipole/inputs_true.dat @@ -1 +1 @@ -2ad86dbb798ab68b45ca535fed3bd848625e0a2b79c68d11db91b14c74b8a48ed4d0129de7e6f6b6b7e9cb47a55d45aadca5ef535dc8a18cb881dddfc74d0bbb \ No newline at end of file +8462e17d102259b3a48a7e908bc75038a28a19d4a5e8bd38f26579e5baf3998bc2d05d1d3055ac1ed478d0c01bd64868d654919c2e6ca3fea86d79f03eda25ef \ No newline at end of file diff --git a/tests/test_multipole/results_true.dat b/tests/test_multipole/results_true.dat index d138fa16a..9555a2d53 100644 --- a/tests/test_multipole/results_true.dat +++ b/tests/test_multipole/results_true.dat @@ -1,5 +1,5 @@ k-combined: -1.457760E+00 1.119659E-02 +1.457760E+00 1.119656E-02 Cell ID = 11 Name = diff --git a/tests/test_multipole/test_multipole.py b/tests/test_multipole/test_multipole.py index 812f5b6d0..a44489aa5 100644 --- a/tests/test_multipole/test_multipole.py +++ b/tests/test_multipole/test_multipole.py @@ -7,7 +7,6 @@ import openmc from openmc.stats import Box from openmc.source import Source - class MultipoleTestHarness(PyAPITestHarness): def _build_inputs(self): #################### @@ -24,29 +23,20 @@ class MultipoleTestHarness(PyAPITestHarness): dense_fuel.add_nuclide('U235', 1.0) mats_file = openmc.Materials([moderator, dense_fuel]) - mats_file.default_temperature = '294K' mats_file.export_to_xml() - #################### # Geometry #################### - c1 = openmc.Cell(cell_id=1) - c1.fill = moderator - mod_univ = openmc.Universe(universe_id=1) - mod_univ.add_cell(c1) + c1 = openmc.Cell(cell_id=1, fill=moderator) + mod_univ = openmc.Universe(universe_id=1, cells=(c1,)) r0 = openmc.ZCylinder(R=0.3) - c11 = openmc.Cell(cell_id=11) - c11.region = -r0 - c11.fill = dense_fuel + c11 = openmc.Cell(cell_id=11, fill=dense_fuel, region=-r0) c11.temperature = [500, 0, 700, 800] - c12 = openmc.Cell(cell_id=12) - c12.region = +r0 - c12.fill = moderator - fuel_univ = openmc.Universe(universe_id=11) - fuel_univ.add_cells((c11, c12)) + c12 = openmc.Cell(cell_id=12, fill=moderator, region=+r0) + fuel_univ = openmc.Universe(universe_id=11, cells=(c11, c12)) lat = openmc.RectLattice(lattice_id=101) lat.dimension = [2, 2] @@ -61,17 +51,12 @@ class MultipoleTestHarness(PyAPITestHarness): y1 = openmc.YPlane(y0=3.0) for s in [x0, x1, y0, y1]: s.boundary_type = 'reflective' - c101 = openmc.Cell(cell_id=101) - c101.region = +x0 & -x1 & +y0 & -y1 - c101.fill = lat - root_univ = openmc.Universe(universe_id=0) - root_univ.add_cell(c101) + c101 = openmc.Cell(cell_id=101, fill=lat, region=+x0 & -x1 & +y0 & -y1) + root_univ = openmc.Universe(universe_id=0, cells=(c101,)) - geometry = openmc.Geometry() - geometry.root_universe = root_univ + geometry = openmc.Geometry(root_univ) geometry.export_to_xml() - #################### # Settings #################### @@ -82,10 +67,9 @@ class MultipoleTestHarness(PyAPITestHarness): sets_file.particles = 1000 sets_file.source = Source(space=Box([-1, -1, -1], [1, 1, 1])) sets_file.output = {'summary': True} - sets_file.use_windowed_multipole=True + sets_file.temperature = {'method': 'multipole'} sets_file.export_to_xml() - #################### # Plots #################### diff --git a/tests/test_natural_element/materials.xml b/tests/test_natural_element/materials.xml index 543b5351d..6568951f4 100644 --- a/tests/test_natural_element/materials.xml +++ b/tests/test_natural_element/materials.xml @@ -3,8 +3,6 @@ - 294K - diff --git a/tests/test_output/materials.xml b/tests/test_output/materials.xml index 37cfcdc6e..2472a7471 100644 --- a/tests/test_output/materials.xml +++ b/tests/test_output/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_output/settings.xml b/tests/test_output/settings.xml index fef5e2fcf..e2f3fc018 100644 --- a/tests/test_output/settings.xml +++ b/tests/test_output/settings.xml @@ -1,7 +1,7 @@ - + 10 diff --git a/tests/test_output/test_output.py b/tests/test_output/test_output.py index 8e36ead80..81f8f42c8 100644 --- a/tests/test_output/test_output.py +++ b/tests/test_output/test_output.py @@ -19,10 +19,6 @@ class OutputTestHarness(TestHarness): assert summary[0].endswith('h5'),\ 'Summary file is not a HDF5 file.' - # Check for the cross sections. - assert os.path.exists(os.path.join(os.getcwd(), 'cross_sections.out')),\ - 'Cross section output file does not exist.' - def _cleanup(self): TestHarness._cleanup(self) output = glob.glob(os.path.join(os.getcwd(), 'summary.*')) diff --git a/tests/test_particle_restart_eigval/materials.xml b/tests/test_particle_restart_eigval/materials.xml index 301efd2bd..3aa37fca6 100644 --- a/tests/test_particle_restart_eigval/materials.xml +++ b/tests/test_particle_restart_eigval/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_particle_restart_fixed/materials.xml b/tests/test_particle_restart_fixed/materials.xml index 95afc3999..f3851d7ef 100644 --- a/tests/test_particle_restart_fixed/materials.xml +++ b/tests/test_particle_restart_fixed/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_plot/materials.xml b/tests/test_plot/materials.xml index 0d5ee77ce..90b354267 100644 --- a/tests/test_plot/materials.xml +++ b/tests/test_plot/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_ptables_off/materials.xml b/tests/test_ptables_off/materials.xml index 37cfcdc6e..2472a7471 100644 --- a/tests/test_ptables_off/materials.xml +++ b/tests/test_ptables_off/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_quadric_surfaces/materials.xml b/tests/test_quadric_surfaces/materials.xml index baf523b34..f68768383 100644 --- a/tests/test_quadric_surfaces/materials.xml +++ b/tests/test_quadric_surfaces/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_reflective_plane/materials.xml b/tests/test_reflective_plane/materials.xml index 37cfcdc6e..2472a7471 100644 --- a/tests/test_reflective_plane/materials.xml +++ b/tests/test_reflective_plane/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_resonance_scattering/inputs_true.dat b/tests/test_resonance_scattering/inputs_true.dat index ba59e03f5..7b515cd1c 100644 --- a/tests/test_resonance_scattering/inputs_true.dat +++ b/tests/test_resonance_scattering/inputs_true.dat @@ -1 +1 @@ -dea135e820f605cf7a7f8d184330d2dd093cbcd2d53a71198819cb4dc7c455809dbb37eff79538d9b35d19932bd88b555ba833f32d01fa8d5289840b2250f72a \ No newline at end of file +15d4ce20d34fbafc757689f1a1f014c42b25efd1cac6bbefdd284c5c1af8bca264c75055fef719441655a61201d0ef098c82cacaa2aeea2b410d3006e983a942 \ No newline at end of file diff --git a/tests/test_resonance_scattering/test_resonance_scattering.py b/tests/test_resonance_scattering/test_resonance_scattering.py index 0e34a4a4c..3daf5a387 100644 --- a/tests/test_resonance_scattering/test_resonance_scattering.py +++ b/tests/test_resonance_scattering/test_resonance_scattering.py @@ -9,16 +9,21 @@ import openmc class ResonanceScatteringTestHarness(PyAPITestHarness): def _build_inputs(self): + # Nuclides + u238 = openmc.Nuclide('U238') + u235 = openmc.Nuclide('U235') + pu239 = openmc.Nuclide('Pu239') + h1 = openmc.Nuclide('H1') + # Materials mat = openmc.Material(material_id=1) mat.set_density('g/cc', 1.0) - mat.add_nuclide('U238', 1.0) - mat.add_nuclide('U235', 0.02) - mat.add_nuclide('Pu239', 0.02) - mat.add_nuclide('H1', 20.0) + mat.add_nuclide(u238, 1.0) + mat.add_nuclide(u235, 0.02) + mat.add_nuclide(pu239, 0.02) + mat.add_nuclide(h1, 20.0) mats_file = openmc.Materials([mat]) - mats_file.default_temperature = '294K' mats_file.export_to_xml() # Geometry @@ -37,29 +42,9 @@ class ResonanceScatteringTestHarness(PyAPITestHarness): geometry.export_to_xml() # Settings - nuclide = openmc.Nuclide('U238') - res_scatt_dbrc = openmc.ResonanceScattering() - res_scatt_dbrc.nuclide = nuclide - res_scatt_dbrc.nuclide_0K = nuclide # This is a bad idea! Just for tests - res_scatt_dbrc.method = 'DBRC' - res_scatt_dbrc.E_min = 1e-6 - res_scatt_dbrc.E_max = 210e-6 - - nuclide = openmc.Nuclide('U235') - res_scatt_wcm = openmc.ResonanceScattering() - res_scatt_wcm.nuclide = nuclide - res_scatt_wcm.nuclide_0K = nuclide - res_scatt_wcm.method = 'WCM' - res_scatt_wcm.E_min = 1e-6 - res_scatt_wcm.E_max = 210e-6 - - nuclide = openmc.Nuclide('Pu239') - res_scatt_ares = openmc.ResonanceScattering() - res_scatt_ares.nuclide = nuclide - res_scatt_ares.nuclide_0K = nuclide - res_scatt_ares.method = 'ARES' - res_scatt_ares.E_min = 1e-6 - res_scatt_ares.E_max = 210e-6 + res_scatt_dbrc = openmc.ResonanceScattering(u238, 'DBRC', 1e-6, 210e-6) + res_scatt_wcm = openmc.ResonanceScattering(u235, 'WCM', 1e-6, 210e-6) + res_scatt_ares = openmc.ResonanceScattering(pu239, 'ARES', 1e-6, 210e-6) sets_file = openmc.Settings() sets_file.batches = 10 diff --git a/tests/test_rotation/materials.xml b/tests/test_rotation/materials.xml index 37cfcdc6e..2472a7471 100644 --- a/tests/test_rotation/materials.xml +++ b/tests/test_rotation/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_salphabeta/materials.xml b/tests/test_salphabeta/materials.xml index 0d250e7b4..bfe0a6224 100644 --- a/tests/test_salphabeta/materials.xml +++ b/tests/test_salphabeta/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_score_current/materials.xml b/tests/test_score_current/materials.xml index bead58f8b..8021f5f99 100644 --- a/tests/test_score_current/materials.xml +++ b/tests/test_score_current/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_seed/materials.xml b/tests/test_seed/materials.xml index 37cfcdc6e..2472a7471 100644 --- a/tests/test_seed/materials.xml +++ b/tests/test_seed/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_source/inputs_true.dat b/tests/test_source/inputs_true.dat index 5744adeeb..0c8f18194 100644 --- a/tests/test_source/inputs_true.dat +++ b/tests/test_source/inputs_true.dat @@ -1 +1 @@ -29498faa9496b8eeeab79f6cb5a966cb03a54e9c3c7e9178e7cc132df575f1377aa780fa2684201c4becdf199f89c3e7c36a824a495e60359bc6fd30b02cad6f \ No newline at end of file +b791dc4a37d20599afe375d91a9c6da123d6579cdc3ae7093a4a464c82bce43f2349fc6828f0b5d9f6b11b8a1aa1031459cd19f2771277ef40d775fedb68c00f \ No newline at end of file diff --git a/tests/test_source/test_source.py b/tests/test_source/test_source.py index f83e39974..ce9012bc2 100644 --- a/tests/test_source/test_source.py +++ b/tests/test_source/test_source.py @@ -13,7 +13,7 @@ import openmc class SourceTestHarness(PyAPITestHarness): def _build_inputs(self): - mat1 = openmc.Material(material_id=1, temperature='294K') + mat1 = openmc.Material(material_id=1, temperature='294') mat1.set_density('g/cm3', 4.5) mat1.add_nuclide(openmc.Nuclide('U235'), 1.0) materials = openmc.Materials([mat1]) diff --git a/tests/test_source_file/materials.xml b/tests/test_source_file/materials.xml index 37cfcdc6e..2472a7471 100644 --- a/tests/test_source_file/materials.xml +++ b/tests/test_source_file/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_sourcepoint_batch/materials.xml b/tests/test_sourcepoint_batch/materials.xml index 37cfcdc6e..2472a7471 100644 --- a/tests/test_sourcepoint_batch/materials.xml +++ b/tests/test_sourcepoint_batch/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_sourcepoint_interval/materials.xml b/tests/test_sourcepoint_interval/materials.xml index 37cfcdc6e..2472a7471 100644 --- a/tests/test_sourcepoint_interval/materials.xml +++ b/tests/test_sourcepoint_interval/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_sourcepoint_latest/materials.xml b/tests/test_sourcepoint_latest/materials.xml index 37cfcdc6e..2472a7471 100644 --- a/tests/test_sourcepoint_latest/materials.xml +++ b/tests/test_sourcepoint_latest/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_sourcepoint_restart/materials.xml b/tests/test_sourcepoint_restart/materials.xml index 37cfcdc6e..2472a7471 100644 --- a/tests/test_sourcepoint_restart/materials.xml +++ b/tests/test_sourcepoint_restart/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_statepoint_batch/materials.xml b/tests/test_statepoint_batch/materials.xml index 37cfcdc6e..2472a7471 100644 --- a/tests/test_statepoint_batch/materials.xml +++ b/tests/test_statepoint_batch/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_statepoint_interval/materials.xml b/tests/test_statepoint_interval/materials.xml index 37cfcdc6e..2472a7471 100644 --- a/tests/test_statepoint_interval/materials.xml +++ b/tests/test_statepoint_interval/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_statepoint_restart/materials.xml b/tests/test_statepoint_restart/materials.xml index 37cfcdc6e..2472a7471 100644 --- a/tests/test_statepoint_restart/materials.xml +++ b/tests/test_statepoint_restart/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_statepoint_sourcesep/materials.xml b/tests/test_statepoint_sourcesep/materials.xml index 37cfcdc6e..2472a7471 100644 --- a/tests/test_statepoint_sourcesep/materials.xml +++ b/tests/test_statepoint_sourcesep/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_survival_biasing/materials.xml b/tests/test_survival_biasing/materials.xml index 03b2162b3..f271ddee2 100644 --- a/tests/test_survival_biasing/materials.xml +++ b/tests/test_survival_biasing/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_tally_assumesep/materials.xml b/tests/test_tally_assumesep/materials.xml index bead58f8b..8021f5f99 100644 --- a/tests/test_tally_assumesep/materials.xml +++ b/tests/test_tally_assumesep/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_tally_nuclides/materials.xml b/tests/test_tally_nuclides/materials.xml index 5ac4f6942..1f89c7df6 100644 --- a/tests/test_tally_nuclides/materials.xml +++ b/tests/test_tally_nuclides/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_trace/materials.xml b/tests/test_trace/materials.xml index 37cfcdc6e..2472a7471 100644 --- a/tests/test_trace/materials.xml +++ b/tests/test_trace/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_track_output/materials.xml b/tests/test_track_output/materials.xml index f31f61fc3..5dc9a6475 100644 --- a/tests/test_track_output/materials.xml +++ b/tests/test_track_output/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_translation/materials.xml b/tests/test_translation/materials.xml index 37cfcdc6e..2472a7471 100644 --- a/tests/test_translation/materials.xml +++ b/tests/test_translation/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_trigger_batch_interval/materials.xml b/tests/test_trigger_batch_interval/materials.xml index 5ac4f6942..1f89c7df6 100644 --- a/tests/test_trigger_batch_interval/materials.xml +++ b/tests/test_trigger_batch_interval/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_trigger_no_batch_interval/materials.xml b/tests/test_trigger_no_batch_interval/materials.xml index 5ac4f6942..1f89c7df6 100644 --- a/tests/test_trigger_no_batch_interval/materials.xml +++ b/tests/test_trigger_no_batch_interval/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_trigger_no_status/materials.xml b/tests/test_trigger_no_status/materials.xml index 5ac4f6942..1f89c7df6 100644 --- a/tests/test_trigger_no_status/materials.xml +++ b/tests/test_trigger_no_status/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_trigger_tallies/materials.xml b/tests/test_trigger_tallies/materials.xml index 5ac4f6942..1f89c7df6 100644 --- a/tests/test_trigger_tallies/materials.xml +++ b/tests/test_trigger_tallies/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_uniform_fs/materials.xml b/tests/test_uniform_fs/materials.xml index 37cfcdc6e..2472a7471 100644 --- a/tests/test_uniform_fs/materials.xml +++ b/tests/test_uniform_fs/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_union_energy_grids/geometry.xml b/tests/test_union_energy_grids/geometry.xml deleted file mode 100644 index bc56030e1..000000000 --- a/tests/test_union_energy_grids/geometry.xml +++ /dev/null @@ -1,8 +0,0 @@ - - - - - - - - diff --git a/tests/test_union_energy_grids/materials.xml b/tests/test_union_energy_grids/materials.xml deleted file mode 100644 index 2946a8d60..000000000 --- a/tests/test_union_energy_grids/materials.xml +++ /dev/null @@ -1,13 +0,0 @@ - - - - 294K - - - - - - - - - diff --git a/tests/test_union_energy_grids/results_true.dat b/tests/test_union_energy_grids/results_true.dat deleted file mode 100644 index 0a607592c..000000000 --- a/tests/test_union_energy_grids/results_true.dat +++ /dev/null @@ -1,2 +0,0 @@ -k-combined: -3.330789E-01 2.216495E-03 diff --git a/tests/test_union_energy_grids/settings.xml b/tests/test_union_energy_grids/settings.xml deleted file mode 100644 index 1eb22241c..000000000 --- a/tests/test_union_energy_grids/settings.xml +++ /dev/null @@ -1,18 +0,0 @@ - - - - union - - - 10 - 5 - 1000 - - - - - -4 -4 -4 4 4 4 - - - - diff --git a/tests/test_union_energy_grids/test_union_energy_grids.py b/tests/test_union_energy_grids/test_union_energy_grids.py deleted file mode 100644 index 2a595f3e6..000000000 --- a/tests/test_union_energy_grids/test_union_energy_grids.py +++ /dev/null @@ -1,11 +0,0 @@ -#!/usr/bin/env python - -import os -import sys -sys.path.insert(0, os.pardir) -from testing_harness import TestHarness - - -if __name__ == '__main__': - harness = TestHarness('statepoint.10.*') - harness.main() diff --git a/tests/test_universe/materials.xml b/tests/test_universe/materials.xml index 37cfcdc6e..2472a7471 100644 --- a/tests/test_universe/materials.xml +++ b/tests/test_universe/materials.xml @@ -1,8 +1,6 @@ - 294K - diff --git a/tests/test_void/materials.xml b/tests/test_void/materials.xml index 2a578e0ed..f70c3a40f 100644 --- a/tests/test_void/materials.xml +++ b/tests/test_void/materials.xml @@ -3,8 +3,6 @@ - 294K -