crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Major refactor to include cross sections at multiple temperatures within the

Browse source 
Nuclide and SAlphaData derived types.
This commit is contained in:
Paul Romano 2016-08-29 10:16:59 -05:00
parent 4448d17ee1
commit a95266ab22
121 changed files with 1670 additions and 2213 deletions
Download patch file Download diff file Expand all files Collapse all files

23
docs/source/io_formats/nuclear_data.rst
View file

@ -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: - **<TTT>K** (*double*) -- kT values (in MeV) for each Temperature
TTT (in Kelvin)
:Datasets:
- **<TTT>K** (*double*) -- kT values (in MeV) for each Temperature
TTT (in Kelvin)
**/<nuclide name>/reactions/reaction_<mt>/**
@ -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)
**/<nuclide name>/reactions/reaction_<mt>/product_<j>/**
@ -112,7 +116,7 @@ Thermal Neutron Scattering Data
**/<thermal name>/**
: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: - **<TTT>K** (*double*) -- kT values (in MeV) for each Temperature
TTT (in Kelvin)
:Datasets:
- **<TTT>K** (*double*) -- kT values (in MeV) for each Temperature
TTT (in Kelvin)
**/<thermal name>/elastic/<TTT>K/**

37
docs/source/methods/cross_sections.rst
View file

@ -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

73
docs/source/usersguide/input.rst
View file

@ -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:
``<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 <use_windowed_multipole> 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.
``<max_order>`` Element
---------------------------
@ -707,6 +709,47 @@ survival biasing, otherwise known as implicit capture or absorption.
*Default*: false
.. _temperature_default:
``<temperature_default>`` Element
---------------------------------
The ``<temperature_default>`` 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:
``<temperature_method>`` Element
--------------------------------
The ``<temperature_method>`` 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
<multipole_library>` must also be available).
*Default*: "nearest"
.. _temperature_tolerance:
``<temperature_tolerance>`` Element
-----------------------------------
The ``<temperature_tolerance>`` 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
``<threads>`` Element
---------------------
@ -1083,7 +1126,9 @@ Each ``<cell>`` 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 <temperature_default>` 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
<temperature_default>` 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
``<default_temperature>`` 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
``<temperature>`` element on every material, a ``<default_temperature>``
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
------------------------------------

1
examples/python/basic/build-xml.py
View file

@ -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()

1
examples/python/boxes/build-xml.py
View file

@ -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()

1
examples/python/lattice/hexagonal/build-xml.py
View file

@ -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()

1
examples/python/lattice/nested/build-xml.py
View file

@ -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()

1
examples/python/lattice/simple/build-xml.py
View file

@ -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()

1
examples/python/pincell/build-xml.py
View file

@ -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()

9
examples/python/pincell_multigroup/build-xml.py
View file

@ -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()

1
examples/python/reflective/build-xml.py
View file

@ -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()

2
examples/xml/basic/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="40">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
examples/xml/boxes/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
examples/xml/lattice/nested/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<!-- Definition of materials -->
<material id="1">
<density value="4.5" units="g/cc" />

2
examples/xml/lattice/simple/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<!-- Definition of materials -->
<material id="1">
<density value="4.5" units="g/cc" />

3
examples/xml/pincell/materials.xml
View file

@ -1,9 +1,6 @@
<?xml version="1.0"?>
<materials>
<!-- By default, use 294K cross sections -->
<default_temperature>294K</default_temperature>
<!--
Since O-18 is not present in ENDF/B-VII, it was necessary to combine the
atom densities for O-17 and O-18 in any materials containing Oxygen.

3
examples/xml/pincell_multigroup/materials.xml
View file

@ -1,8 +1,5 @@
<?xml version="1.0"?>
<materials>
<!-- Set default xs set to use 300K data -->
<default_xs>300K</default_xs>
<!-- UO2 -->
<material id="1">
<density units="macro" value="1.0" />

32
examples/xml/pincell_multigroup/mg_cross_sections.xml
View file

@ -11,8 +11,8 @@
-->
<xsdata>
<!-- Meta data for this data -->
<name>UO2.300K</name>
<alias>UO2.300K</alias>
<name>UO2</name>
<alias>UO2</alias>
<kT> 2.53E-8 </kT> <!-- in MeV -->
<order>0</order>
<fissionable>true</fissionable>
@ -67,8 +67,8 @@
<xsdata>
<!-- Meta data for this data -->
<name>MOX1.300K</name>
<alias>MOX1.300K</alias>
<name>MOX1</name>
<alias>MOX1</alias>
<kT> 2.53E-8 </kT> <!-- in MeV -->
<order>0</order>
<fissionable>true</fissionable>
@ -124,8 +124,8 @@
<xsdata>
<!-- Meta data for this data -->
<name>MOX2.300K</name>
<alias>MOX2.300K</alias>
<name>MOX2</name>
<alias>MOX2</alias>
<kT> 2.53E-8 </kT> <!-- in MeV -->
<order>0</order>
<fissionable>true</fissionable>
@ -180,8 +180,8 @@
<xsdata>
<!-- Meta data for this data -->
<name>MOX3.300K</name>
<alias>MOX3.300K</alias>
<name>MOX3</name>
<alias>MOX3</alias>
<kT> 2.53E-8 </kT> <!-- in MeV -->
<order>0</order>
<fissionable>true</fissionable>
@ -236,8 +236,8 @@
<xsdata>
<!-- Meta data for this data -->
<name>FC.300K</name>
<alias>FC.300K</alias>
<name>FC</name>
<alias>FC</alias>
<kT> 2.53E-8 </kT> <!-- in MeV -->
<order>0</order>
<fissionable>true</fissionable>
@ -286,8 +286,8 @@
<xsdata>
<!-- Meta data for this data -->
<name>GT.300K</name>
<alias>GT.300K</alias>
<name>GT</name>
<alias>GT</alias>
<kT> 2.53E-8 </kT> <!-- in MeV -->
<order>0</order>
<fissionable>false</fissionable>
@ -318,8 +318,8 @@
<xsdata>
<!-- Meta data for this data -->
<name>LWTR.300K</name>
<alias>LWTR.300K</alias>
<name>LWTR</name>
<alias>LWTR</alias>
<kT> 2.53E-8 </kT> <!-- in MeV -->
<order>0</order>
<fissionable>false</fissionable>
@ -351,8 +351,8 @@
<xsdata>
<!-- Meta data for this data -->
<name>CR.300K</name>
<alias>CR.300K</alias>
<name>CR</name>
<alias>CR</alias>
<kT> 2.53E-8 </kT> <!-- in MeV -->
<order>0</order>
<fissionable>false</fissionable>

2
examples/xml/reflective/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

9
openmc/data/neutron.py
View file

@ -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)

24
openmc/data/thermal.py
View file

@ -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

23
openmc/mgxs_library.py
View file

@ -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)

16
openmc/nuclide.py
View file

@ -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)

122
openmc/settings.py
View file

@ -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'))

4
openmc/summary.py
View file

@ -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

12
src/constants.F90
View file

@ -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

511
src/cross_section.F90
View file

@ -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

176
src/energy_grid.F90
View file

@ -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

52
src/geometry.F90
View file

@ -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

12
src/global.F90
View file

@ -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.
! ============================================================================

62
src/hdf5_interface.F90
View file

@ -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(:)

76
src/initialize.F90
View file

@ -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

291
src/input_xml.F90
View file

@ -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 <use_windowed_multipole> 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 <material>
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

3
src/material_header.F90
View file

@ -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

12
src/mgxs_header.F90
View file

@ -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

110
src/multipole.F90
View file

@ -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)

524
src/nuclide_header.F90
View file

@ -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

51
src/output.F90
View file

@ -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

2
src/particle_header.F90
View file

@ -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

454
src/physics.F90
View file

@ -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)

36
src/reaction_header.F90
View file

@ -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)

6
src/relaxng/materials.rnc
View file

@ -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" } }?
}+
}

307
src/relaxng/materials.rng
View file

@ -1,197 +1,186 @@
<?xml version="1.0" encoding="UTF-8"?>
<element name="materials" xmlns="http://relaxng.org/ns/structure/1.0" datatypeLibrary="http://www.w3.org/2001/XMLSchema-datatypes">
<interleave>
<oneOrMore>
<element name="material">
<interleave>
<oneOrMore>
<element name="material">
<interleave>
<choice>
<element name="id">
<data type="int"/>
</element>
<attribute name="id">
<data type="int"/>
</attribute>
</choice>
<optional>
<choice>
<element name="id">
<data type="int"/>
<element name="name">
<data type="string">
<param name="maxLength">52</param>
</data>
</element>
<attribute name="id">
<data type="int"/>
<attribute name="name">
<data type="string">
<param name="maxLength">52</param>
</data>
</attribute>
</choice>
<optional>
</optional>
<optional>
<element name="temperature">
<data type="double"/>
</element>
</optional>
<element name="density">
<interleave>
<optional>
<choice>
<element name="value">
<data type="double"/>
</element>
<attribute name="value">
<data type="double"/>
</attribute>
</choice>
</optional>
<choice>
<element name="name">
<element name="units">
<data type="string">
<param name="maxLength">52</param>
<param name="maxLength">10</param>
</data>
</element>
<attribute name="name">
<attribute name="units">
<data type="string">
<param name="maxLength">52</param>
<param name="maxLength">10</param>
</data>
</attribute>
</choice>
</optional>
<optional>
<element name="temperature">
<data type="string">
<param name="maxLength">6</param>
</data>
</element>
</optional>
<element name="density">
</interleave>
</element>
<zeroOrMore>
<element name="nuclide">
<interleave>
<choice>
<element name="name">
<data type="string"/>
</element>
<attribute name="name">
<data type="string"/>
</attribute>
</choice>
<optional>
<choice>
<element name="value">
<data type="double"/>
<element name="scattering">
<choice>
<value>data</value>
<value>iso-in-lab</value>
</choice>
</element>
<attribute name="value">
<data type="double"/>
<attribute name="scattering">
<choice>
<value>data</value>
<value>iso-in-lab</value>
</choice>
</attribute>
</choice>
</optional>
<choice>
<element name="units">
<data type="string">
<param name="maxLength">10</param>
</data>
</element>
<attribute name="units">
<data type="string">
<param name="maxLength">10</param>
</data>
</attribute>
<choice>
<element name="ao">
<data type="double"/>
</element>
<attribute name="ao">
<data type="double"/>
</attribute>
</choice>
<choice>
<element name="wo">
<data type="double"/>
</element>
<attribute name="wo">
<data type="double"/>
</attribute>
</choice>
</choice>
</interleave>
</element>
<zeroOrMore>
<element name="nuclide">
<interleave>
<choice>
<element name="name">
<data type="string"/>
</element>
<attribute name="name">
<data type="string"/>
</attribute>
</choice>
<optional>
<choice>
<element name="scattering">
<choice>
<value>data</value>
<value>iso-in-lab</value>
</choice>
</element>
<attribute name="scattering">
<choice>
<value>data</value>
<value>iso-in-lab</value>
</choice>
</attribute>
</choice>
</optional>
<choice>
<choice>
<element name="ao">
<data type="double"/>
</element>
<attribute name="ao">
<data type="double"/>
</attribute>
</choice>
<choice>
<element name="wo">
<data type="double"/>
</element>
<attribute name="wo">
<data type="double"/>
</attribute>
</choice>
</choice>
</interleave>
</element>
</zeroOrMore>
<zeroOrMore>
<element name="macroscopic">
</zeroOrMore>
<zeroOrMore>
<element name="macroscopic">
<choice>
<element name="name">
<data type="string"/>
</element>
<attribute name="name">
<data type="string"/>
</attribute>
</choice>
</element>
</zeroOrMore>
<zeroOrMore>
<element name="element">
<interleave>
<choice>
<element name="name">
<data type="string"/>
<data type="string">
<param name="maxLength">2</param>
</data>
</element>
<attribute name="name">
<data type="string"/>
<data type="string">
<param name="maxLength">2</param>
</data>
</attribute>
</choice>
</element>
</zeroOrMore>
<zeroOrMore>
<element name="element">
<interleave>
<optional>
<choice>
<element name="name">
<data type="string">
<param name="maxLength">2</param>
</data>
<element name="scattering">
<choice>
<value>data</value>
<value>iso-in-lab</value>
</choice>
</element>
<attribute name="name">
<data type="string">
<param name="maxLength">2</param>
</data>
<attribute name="scattering">
<choice>
<value>data</value>
<value>iso-in-lab</value>
</choice>
</attribute>
</choice>
<optional>
<choice>
<element name="scattering">
<choice>
<value>data</value>
<value>iso-in-lab</value>
</choice>
</element>
<attribute name="scattering">
<choice>
<value>data</value>
<value>iso-in-lab</value>
</choice>
</attribute>
</choice>
</optional>
<choice>
<choice>
<element name="ao">
<data type="double"/>
</element>
<attribute name="ao">
<data type="double"/>
</attribute>
</choice>
<choice>
<element name="wo">
<data type="double"/>
</element>
<attribute name="wo">
<data type="double"/>
</attribute>
</choice>
</choice>
</interleave>
</element>
</zeroOrMore>
<zeroOrMore>
<element name="sab">
</optional>
<choice>
<element name="name">
<data type="string"/>
</element>
<attribute name="name">
<data type="string"/>
</attribute>
<choice>
<element name="ao">
<data type="double"/>
</element>
<attribute name="ao">
<data type="double"/>
</attribute>
</choice>
<choice>
<element name="wo">
<data type="double"/>
</element>
<attribute name="wo">
<data type="double"/>
</attribute>
</choice>
</choice>
</element>
</zeroOrMore>
</interleave>
</element>
</oneOrMore>
<optional>
<element name="default_temperature">
<data type="string">
<param name="maxLength">6</param>
</data>
</element>
</optional>
</interleave>
</interleave>
</element>
</zeroOrMore>
<zeroOrMore>
<element name="sab">
<choice>
<element name="name">
<data type="string"/>
</element>
<attribute name="name">
<data type="string"/>
</attribute>
</choice>
</element>
</zeroOrMore>
</interleave>
</element>
</oneOrMore>
</element>

6
src/relaxng/settings.rnc
View file

@ -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+ } }? &

15
src/relaxng/settings.rng
View file

@ -565,6 +565,21 @@
<data type="boolean"/>
</element>
</optional>
<optional>
<element name="temperature_default">
<data type="double"/>
</element>
</optional>
<optional>
<element name="temperature_method">
<data type="string"/>
</element>
</optional>
<optional>
<element name="temperature_tolerance">
<data type="double"/>
</element>
</optional>
<optional>
<element name="threads">
<data type="positiveInteger"/>

428
src/sab_header.F90
View file

@ -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

14
src/summary.F90
View file

@ -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", &

42
src/tally.F90
View file

@ -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

7
src/tracking.F90
View file

@ -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

4
tests/test_cmfd_feed/materials.xml
View file

@ -1,9 +1,7 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<!-- Definition of materials -->
<!-- Definition of materials -->
<material id="1">
<density value="19" units="g/cc" />
<nuclide name="U235" wo="0.21" />

4
tests/test_cmfd_nofeed/materials.xml
View file

@ -1,9 +1,7 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<!-- Definition of materials -->
<!-- Definition of materials -->
<material id="1">
<density value="19" units="g/cc" />
<nuclide name="U235" wo="0.21" />

2
tests/test_complex_cell/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
tests/test_confidence_intervals/materials.xml
View file

@ -2,7 +2,7 @@
<materials>
<material id="1">
<temperature>294K</temperature>
<temperature>294</temperature>
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />
</material>

2
tests/test_density/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="0.1" units="atom/b-cm" />
<nuclide name="U235" ao="1.0" />

2
tests/test_eigenvalue_genperbatch/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
tests/test_eigenvalue_no_inactive/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
tests/test_energy_grid/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
tests/test_energy_grid/settings.xml
View file

@ -1,7 +1,7 @@
<?xml version="1.0"?>
<settings>
<energy_grid>nuclide</energy_grid>
<log_grid_bins>20000</log_grid_bins>
<eigenvalue>
<batches>10</batches>

1
tests/test_energy_laws/materials.xml
View file

@ -1,6 +1,5 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="20" units="g/cc" />
<nuclide name="U233" ao="1.0" />

2
tests/test_entropy/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
tests/test_filter_distribcell/case-1/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<!-- Definition of materials -->
<material id="1">
<density value="4.5" units="g/cc" />

3
tests/test_filter_distribcell/case-2/materials.xml
View file

@ -1,9 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<!-- Definition of materials -->
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
tests/test_filter_distribcell/case-3/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<!-- Fuel composition -->
<material id="1">
<density value="10.062" units="g/cm3" />

2
tests/test_filter_distribcell/case-4/materials.xml
View file

@ -1,7 +1,5 @@
<?xml version='1.0' encoding='utf-8'?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1" name="fuel">
<density units="g/cc" value="4.5" />
<nuclide ao="1.0" name="U235" />

2
tests/test_filter_mesh_2d/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<!-- Fuel composition -->
<material id="1">
<density value="10.062" units="g/cm3" />

2
tests/test_filter_mesh_3d/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<!-- Fuel composition -->
<material id="1">
<density value="10.062" units="g/cm3" />

2
tests/test_fixed_source/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="7.5" units="g/cc" />
<nuclide name="O16" ao="1.0" />

2
tests/test_fixed_source/settings.xml
View file

@ -6,6 +6,8 @@
<particles>100</particles>
</fixed_source>
<temperature_default>294</temperature_default>
<source>
<space type="point" parameters="0 0 0" />
</source>

2
tests/test_infinite_cell/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
tests/test_lattice/materials.xml
View file

@ -10,8 +10,6 @@
===============================================================
-->
<default_temperature>294K</default_temperature>
<!-- UO2-PuO2 Mixture -->
<material id="1">
<density value="9.54" units="g/cm3" />

10
tests/test_lattice_hex/materials.xml
View file

@ -1,9 +1,7 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<!-- Material 01: Water -->
<!-- Material 01: Water -->
<material id="01">
<density value="0.998" units="g/cc"/>
<nuclide name="H1" ao="2.0"/>
@ -11,7 +9,7 @@
<nuclide name="B10" ao="0.00085"/>
</material>
<!-- Material 13: Fuel (3.30%) -->
<!-- Material 13: Fuel (3.30%) -->
<material id="13">
<density value="10.371" units="g/cc"/>
<nuclide name="U235" ao="0.033408078"/>
@ -19,7 +17,7 @@
<nuclide name="O16" ao="2.0"/>
</material>
<!-- Material 21: Zircaloy -->
<!-- Material 21: Zircaloy -->
<material id="21">
<density value="6.44" units="g/cc" />
<nuclide name="Zr90" wo="0.501807"/>
@ -29,7 +27,7 @@
<nuclide name="Zr96" wo="0.029133"/>
</material>
<!-- Material 22: Steel -->
<!-- Material 22: Steel -->
<material id="22">
<density value="7.9" units="g/cc"/>
<nuclide name="Mn55" wo="2"/>

10
tests/test_lattice_mixed/materials.xml
View file

@ -1,9 +1,7 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<!-- Material 01: Water -->
<!-- Material 01: Water -->
<material id="01">
<density value="0.998" units="g/cc"/>
<nuclide name="H1" ao="2.0"/>
@ -11,7 +9,7 @@
<nuclide name="B10" ao="0.00085"/>
</material>
<!-- Material 13: Fuel (3.30%) -->
<!-- Material 13: Fuel (3.30%) -->
<material id="13">
<density value="10.371" units="g/cc"/>
<nuclide name="U235" ao="0.033408078"/>
@ -19,7 +17,7 @@
<nuclide name="O16" ao="2.0"/>
</material>
<!-- Material 21: Zircaloy -->
<!-- Material 21: Zircaloy -->
<material id="21">
<density value="6.44" units="g/cc" />
<nuclide name="Zr90" wo="0.501807"/>
@ -29,7 +27,7 @@
<nuclide name="Zr96" wo="0.029133"/>
</material>
<!-- Material 22: Steel -->
<!-- Material 22: Steel -->
<material id="22">
<density value="7.9" units="g/cc"/>
<nuclide name="Mn55" wo="2"/>

2
tests/test_lattice_multiple/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<!-- Fuel composition -->
<material id="1">
<density value="10.062" units="g/cm3" />

2
tests/test_multipole/inputs_true.dat
View file

@ -1 +1 @@
2ad86dbb798ab68b45ca535fed3bd848625e0a2b79c68d11db91b14c74b8a48ed4d0129de7e6f6b6b7e9cb47a55d45aadca5ef535dc8a18cb881dddfc74d0bbb
8462e17d102259b3a48a7e908bc75038a28a19d4a5e8bd38f26579e5baf3998bc2d05d1d3055ac1ed478d0c01bd64868d654919c2e6ca3fea86d79f03eda25ef

2
tests/test_multipole/results_true.dat
View file

@ -1,5 +1,5 @@
k-combined:
1.457760E+00 1.119659E-02
1.457760E+00 1.119656E-02
Cell
ID = 11
Name =

34
tests/test_multipole/test_multipole.py
View file

@ -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
####################

2
tests/test_natural_element/materials.xml
View file

@ -3,8 +3,6 @@
<!-- pu-met-fast-021 case 2 -->
<default_temperature>294K</default_temperature>
<!-- Plutonium -->
<material id="1">
<density value="18.5345" units="g/cm3" />

2
tests/test_output/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
tests/test_output/settings.xml
View file

@ -1,7 +1,7 @@
<?xml version="1.0"?>
<settings>
<output summary="true" cross_sections="true" />
<output summary="true" />
<eigenvalue>
<batches>10</batches>

4
tests/test_output/test_output.py
View file

@ -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.*'))

2
tests/test_particle_restart_eigval/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="20.0" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
tests/test_particle_restart_fixed/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="1.4" units="g/cc" />
<nuclide name="B10" ao="1.0" />

2
tests/test_plot/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
tests/test_ptables_off/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
tests/test_quadric_surfaces/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U238" ao="1.0" />

2
tests/test_reflective_plane/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
tests/test_resonance_scattering/inputs_true.dat
View file

@ -1 +1 @@
dea135e820f605cf7a7f8d184330d2dd093cbcd2d53a71198819cb4dc7c455809dbb37eff79538d9b35d19932bd88b555ba833f32d01fa8d5289840b2250f72a
15d4ce20d34fbafc757689f1a1f014c42b25efd1cac6bbefdd284c5c1af8bca264c75055fef719441655a61201d0ef098c82cacaa2aeea2b410d3006e983a942

41
tests/test_resonance_scattering/test_resonance_scattering.py
View file

@ -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

2
tests/test_rotation/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
tests/test_salphabeta/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
tests/test_score_current/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<!-- Fuel composition -->
<material id="1">
<density value="10.062" units="g/cm3" />

2
tests/test_seed/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
tests/test_source/inputs_true.dat
View file

@ -1 +1 @@
29498faa9496b8eeeab79f6cb5a966cb03a54e9c3c7e9178e7cc132df575f1377aa780fa2684201c4becdf199f89c3e7c36a824a495e60359bc6fd30b02cad6f
b791dc4a37d20599afe375d91a9c6da123d6579cdc3ae7093a4a464c82bce43f2349fc6828f0b5d9f6b11b8a1aa1031459cd19f2771277ef40d775fedb68c00f

2
tests/test_source/test_source.py
View file

@ -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])

2
tests/test_source_file/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
tests/test_sourcepoint_batch/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
tests/test_sourcepoint_interval/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
tests/test_sourcepoint_latest/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
tests/test_sourcepoint_restart/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
tests/test_statepoint_batch/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

Some files were not shown because too many files have changed in this diff Show more