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

Merge pull request #840 from paulromano/elastic-0k

Browse source 
Ability to add 0K elastic scattering from an ENDF file
This commit is contained in:
Sterling Harper 2017-04-03 12:36:44 -04:00 • committed by GitHub
commit 47deaa8c60
22 changed files with 646 additions and 571 deletions
Download patch file Download diff file Expand all files Collapse all files

2
.travis.yml
View file

@ -43,7 +43,7 @@ install: true
before_script:
- if [[ ! -e $HOME/nndc_hdf5/cross_sections.xml ]]; then
wget https://anl.box.com/shared/static/oy85y1i3gboho82ifdtjlu1q4ozhd7el.xz -O - | tar -C $HOME -xvJ;
wget https://anl.box.com/shared/static/a6sw2cep34wlz6b9i9jwiotaqoayxcxt.xz -O - | tar -C $HOME -xvJ;
fi
- export OPENMC_CROSS_SECTIONS=$HOME/nndc_hdf5/cross_sections.xml

3
docs/source/conf.py
View file

@ -25,7 +25,8 @@ except ImportError:
MOCK_MODULES = ['numpy', 'numpy.polynomial', 'numpy.polynomial.polynomial',
'h5py', 'pandas', 'uncertainties', 'openmoc']
'h5py', 'pandas', 'uncertainties', 'openmoc',
'openmc.data.reconstruct']
sys.modules.update((mod_name, MagicMock()) for mod_name in MOCK_MODULES)
import numpy as np

3
docs/source/pythonapi/index.rst
View file

@ -39,7 +39,6 @@ Simulation Settings
:template: myclass.rst
openmc.Source
openmc.ResonanceScattering
openmc.VolumeCalculation
openmc.Settings
@ -385,6 +384,8 @@ Core Functions
:template: myfunction.rst
openmc.data.atomic_mass
openmc.data.linearize
openmc.data.thin
openmc.data.write_compact_458_library
Angle-Energy Distributions

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

@ -361,54 +361,56 @@ or sub-elements and can be set to either "false" or "true".
``<resonance_scattering>`` Element
----------------------------------
The ``resonance_scattering`` element can contain one or more of the following
attributes or sub-elements:
The ``resonance_scattering`` element indicates to OpenMC that a method be used
to properly account for resonance elastic scattering (typically for nuclides
with Z > 40). This element can contain one or more of the following attributes
or sub-elements:
:scatterer:
An element with attributes/sub-elements called ``nuclide``, ``method``,
``E_min``, and ``E_max``. The ``nuclide`` attribute is the name, as given
by the ``name`` attribute within the ``nuclide`` sub-element of the
``material`` element in ``materials.xml``, of the nuclide to which a
resonance scattering treatment is to be applied.
The ``method`` attribute gives the type of resonance scattering treatment
that is to be applied to the ``nuclide``. Acceptable inputs - none of
which are case-sensitive - for the ``method`` attribute are ``ARES``,
``CXS``, ``WCM``, and ``DBRC``. Descriptions of each of these methods
are documented here_. The ``E_min`` attribute gives the minimum energy
above which the ``method`` is applied. The ``E_max`` attribute gives the
maximum energy below which the ``method`` is applied. One example would
be as follows:
:enable:
Indicates whether a resonance elastic scattering method should be turned
on. Accepts values of "true" or "false".
*Default*: If the ``<resonance_scattering>`` element is present, "true".
:method:
Which resonance elastic scattering method is to be applied: "ares"
(accelerated resonance elastic scattering), "dbrc" (Doppler broadening
rejection correction), or "wcm" (weight correction method). Descriptions of
each of these methods are documented here_.
.. _here: http://dx.doi.org/10.1016/j.anucene.2014.01.017
.. code-block:: xml
*Default*: "ares"
<resonance_scattering>
<scatterer>
<nuclide>U-238</nuclide>
<method>ARES</method>
<E_min>5.0e-6</E_min>
<E_max>40.0e-6</E_max>
</scatterer>
<scatterer>
<nuclide>Pu-239</nuclide>
<method>dbrc</method>
<E_min>0.01e-6</E_min>
<E_max>210.0e-6</E_max>
</scatterer>
</resonance_scattering>
:energy_min:
The energy in eV above which the resonance elastic scattering method should
be applied.
.. note:: If the ``resonance_scattering`` element is not given, the free gas,
constant cross section (``cxs``) scattering model, which has
historically been used by Monte Carlo codes to sample target
velocities, is used to treat the target motion of all nuclides. If
``resonance_scattering`` is present, the ``cxs`` method is applied
below ``E_min`` and the target-at-rest (asymptotic) kernel is used
above ``E_max``. An arbitrary number of ``scatterer`` elements may
be specified, each corresponding to a single nuclide at a single
temperature.
*Default*: 0.01 eV
*Defaults*: None (scatterer), ARES (method), 0.01 eV (E_min), 1.0 keV (E_max)
:energy_max:
The energy in eV below which the resonance elastic scattering method should
be applied.
*Default*: 1000.0 eV
:nuclides:
A list of nuclides to which the resonance elastic scattering method should
be applied.
*Default*: If ``<resonance_scattering>`` is present but the ``<nuclides>``
sub-element is not given, the method is applied to all nuclides with 0 K
elastic scattering data present.
.. note:: If the ``resonance_scattering`` element is not given, the free gas,
constant cross section scattering model, which has historically been
used by Monte Carlo codes to sample target velocities, is used to
treat the target motion of all nuclides. If
``resonance_scattering`` is present, the constant cross section
method is applied below ``energy_min`` and the target-at-rest
(asymptotic) kernel is used above ``energy_max``.
.. note:: This element is not used in the multi-group :ref:`energy_mode`.

1
openmc/data/__init__.py
View file

@ -28,3 +28,4 @@ from .library import *
from .fission_energy import *
from .resonance import *
from .multipole import *
from .grid import *

13
openmc/data/ace.py
View file

@ -24,7 +24,7 @@ from six import string_types
import numpy as np
from openmc.mixin import EqualityMixin
from openmc.data.endf import ENDF_FLOAT_RE
def ascii_to_binary(ascii_file, binary_file):
"""Convert an ACE file in ASCII format (type 1) to binary format (type 2).
@ -344,6 +344,17 @@ class Library(EqualityMixin):
datastr = '0.0 ' + ''.join(lines[12:12+n_lines])
xss = np.fromstring(datastr, sep=' ')
# When NJOY writes an ACE file, any values less than 1e-100 actually
# get written without the 'e'. Thus, what we do here is check
# whether the xss array is of the right size (if a number like
# 1.0-120 is encountered, np.fromstring won't capture any numbers
# after it). If it's too short, then we apply the ENDF float regular
# expression. We don't do this by default because it's expensive!
if xss.size != nxs[1] + 1:
datastr = ENDF_FLOAT_RE.sub(r'\1e\2', datastr)
xss = np.fromstring(datastr, sep=' ')
assert xss.size == nxs[1] + 1
table = Table(name, atomic_weight_ratio, temperature, pairs,
nxs, jxs, xss)
self.tables.append(table)

4
openmc/data/endf.py
View file

@ -46,7 +46,7 @@ SUM_RULES = {1: [2, 3],
106: list(range(750, 800)),
107: list(range(800, 850))}
_ENDF_FLOAT_RE = re.compile(r'([\s\-\+]?\d*\.\d+)([\+\-]\d+)')
ENDF_FLOAT_RE = re.compile(r'([\s\-\+]?\d*\.\d+)([\+\-]\d+)')
def float_endf(s):
@ -69,7 +69,7 @@ def float_endf(s):
The number
"""
return float(_ENDF_FLOAT_RE.sub(r'\1e\2', s))
return float(ENDF_FLOAT_RE.sub(r'\1e\2', s))
def get_text_record(file_obj):

111
openmc/data/grid.py Normal file
View file

@ -0,0 +1,111 @@
import numpy as np
def linearize(x, f, tolerance=0.001):
"""Return a tabulated representation of a function of one variable.
Parameters
----------
x : Iterable of float
Initial x values at which the function should be evaluated
f : Callable
Function of a single variable
tolerance : float
Tolerance on the interpolation error
Returns
-------
numpy.ndarray
Tabulated values of the independent variable
numpy.ndarray
Tabulated values of the dependent variable
"""
# Initialize output arrays
x_out = []
y_out = []
# Initialize stack
x_stack = [x[0]]
y_stack = [f(x[0])]
for i in range(x.shape[0] - 1):
x_stack.insert(0, x[i + 1])
y_stack.insert(0, f(x[i + 1]))
while True:
x_high, x_low = x_stack[-2:]
y_high, y_low = y_stack[-2:]
x_mid = 0.5*(x_low + x_high)
y_mid = f(x_mid)
y_interp = y_low + (y_high - y_low)/(x_high - x_low)*(x_mid - x_low)
error = abs((y_interp - y_mid)/y_mid)
if error > tolerance:
x_stack.insert(-1, x_mid)
y_stack.insert(-1, y_mid)
else:
x_out.append(x_stack.pop())
y_out.append(y_stack.pop())
if len(x_stack) == 1:
break
x_out.append(x_stack.pop())
y_out.append(y_stack.pop())
return np.array(x_out), np.array(y_out)
def thin(x, y, tolerance=0.001):
"""Check for (x,y) points that can be removed.
Parameters
----------
x : numpy.ndarray
Independent variable
y : numpy.ndarray
Dependent variable
tolerance : float
Tolerance on interpolation error
Returns
-------
numpy.ndarray
Tabulated values of the independent variable
numpy.ndarray
Tabulated values of the dependent variable
"""
# Initialize output arrays
x_out = x.copy()
y_out = y.copy()
N = x.shape[0]
i_left = 0
i_right = 2
while i_left < N - 2 and i_right < N:
m = (y[i_right] - y[i_left])/(x[i_right] - x[i_left])
for i in range(i_left + 1, i_right):
# Determine error in interpolated point
y_interp = y[i_left] + m*(x[i] - x[i_left])
if abs(y[i]) > 0.:
error = abs((y_interp - y[i])/y[i])
else:
error = 2*tolerance
if error > tolerance:
for i_remove in range(i_left + 1, i_right - 1):
x_out[i_remove] = np.nan
y_out[i_remove] = np.nan
i_left = i_right - 1
i_right = i_left + 1
break
i_right += 1
for i_remove in range(i_left + 1, i_right - 1):
x_out[i_remove] = np.nan
y_out[i_remove] = np.nan
return x_out[np.isfinite(x_out)], y_out[np.isfinite(y_out)]

118
openmc/data/neutron.py
View file

@ -2,6 +2,7 @@ from __future__ import division, unicode_literals
import sys
from collections import OrderedDict, Iterable, Mapping, MutableMapping
from itertools import chain
from math import log10
from numbers import Integral, Real
from warnings import warn
@ -12,17 +13,22 @@ import h5py
from . import HDF5_VERSION, HDF5_VERSION_MAJOR
from .ace import Table, get_table
from .data import ATOMIC_SYMBOL, K_BOLTZMANN, EV_PER_MEV
from .endf import Evaluation, SUM_RULES
from .fission_energy import FissionEnergyRelease
from .function import Tabulated1D, Sum, ResonancesWithBackground
from .endf import Evaluation, SUM_RULES
from .grid import linearize, thin
from .product import Product
from .reaction import Reaction, _get_photon_products_ace
from .resonance import Resonances, _RESOLVED
from . import resonance as res
from .urr import ProbabilityTables
import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
# Fractions of resonance widths used for reconstructing resonances
_RESONANCE_ENERGY_GRID = np.logspace(-3, 3, 61)
def _get_metadata(zaid, metastable_scheme='nndc'):
"""Return basic identifying data for a nuclide with a given ZAID.
@ -278,7 +284,7 @@ class IncidentNeutron(EqualityMixin):
@resonances.setter
def resonances(self, resonances):
cv.check_type('resonances', resonances, Resonances)
cv.check_type('resonances', resonances, res.Resonances)
self._resonances = resonances
@summed_reactions.setter
@ -343,6 +349,98 @@ class IncidentNeutron(EqualityMixin):
if strT in data.urr:
self.urr[strT] = data.urr[strT]
def add_elastic_0K_from_endf(self, filename, overwrite=False):
"""Append 0K elastic scattering cross section from an ENDF file.
Parameters
----------
filename : str
Path to ENDF file
overwrite : bool
If existing 0 K data is present, this flag can be used to indicate
that it should be overwritten. Otherwise, an exception will be
thrown.
Raises
------
ValueError
If 0 K data is already present and the `overwrite` parameter is
False.
"""
# Check for existing data
if '0K' in self.energy and not overwrite:
raise ValueError('0 K data already exists for this nuclide.')
data = type(self).from_endf(filename)
if data.resonances is not None:
x = []
y = []
for rr in data.resonances:
if isinstance(rr, res.RMatrixLimited):
raise TypeError('R-Matrix Limited not supported.')
elif isinstance(rr, res.Unresolved):
continue
# Get energies/widths for resonances
e_peak = rr.parameters['energy'].values
if isinstance(rr, res.MultiLevelBreitWigner):
gamma = rr.parameters['totalWidth'].values
elif isinstance(rr, res.ReichMoore):
df = rr.parameters
gamma = (df['neutronWidth'] +
df['captureWidth'] +
abs(df['fissionWidthA']) +
abs(df['fissionWidthB'])).values
# Determine peak energies and widths
e_min, e_max = rr.energy_min, rr.energy_max
in_range = (e_peak > e_min) & (e_peak < e_max)
e_peak = e_peak[in_range]
gamma = gamma[in_range]
# Get midpoints between resonances (use min/max energy of
# resolved region as absolute lower/upper bound)
e_mid = np.concatenate(
([e_min], (e_peak[1:] + e_peak[:-1])/2, [e_max]))
# Add grid around each resonance that includes the peak +/- the
# width times each value in _RESONANCE_ENERGY_GRID. Values are
# constrained so that points around one resonance don't overlap
# with points around another. This algorithm is from Fudge.
energies = []
for e, g, e_lower, e_upper in zip(e_peak, gamma, e_mid[:-1],
e_mid[1:]):
e_left = e - g*_RESONANCE_ENERGY_GRID
energies.append(e_left[e_left > e_lower][::-1])
e_right = e + g*_RESONANCE_ENERGY_GRID[1:]
energies.append(e_right[e_right < e_upper])
# Concatenate all points
energies = np.concatenate(energies)
# Create 1000 equal log-spaced energies over RRR, combine with
# resonance peaks and half-height energies
e_log = np.logspace(log10(e_min), log10(e_max), 1000)
energies = np.union1d(e_log, energies)
# Linearize and thin cross section
xi, yi = linearize(energies, data[2].xs['0K'])
xi, yi = thin(xi, yi)
# If there are multiple resolved resonance ranges (e.g. Pu239 in
# ENDF/B-VII.1), combine them
x = np.concatenate((x, xi))
y = np.concatenate((y, yi))
else:
energies = data[2].xs['0K'].x
x, y = linearize(energies, data[2].xs['0K'])
x, y = thin(x, y)
# Set 0K energy grid and elastic scattering cross section
self.energy['0K'] = x
self[2].xs['0K'] = Tabulated1D(x, y)
def get_reaction_components(self, mt):
"""Determine what reactions make up summed reaction.
@ -411,12 +509,22 @@ class IncidentNeutron(EqualityMixin):
for temperature in self.temperatures:
eg.create_dataset(temperature, data=self.energy[temperature])
# Write 0K energy grid if needed
if '0K' in self.energy and '0K' not in eg:
eg.create_dataset('0K', data=self.energy['0K'])
# Write reaction data
rxs_group = g.create_group('reactions')
for rx in self.reactions.values():
rx_group = rxs_group.create_group('reaction_{:03}'.format(rx.mt))
rx.to_hdf5(rx_group)
# Write 0K elastic scattering if needed
if '0K' in rx.xs and '0K' not in rx_group:
group = rx_group.create_group('0K')
dset = group.create_dataset('xs', data=rx.xs['0K'].y)
dset.attrs['threshold_idx'] = 1
# Write total nu data if available
if len(rx.derived_products) > 0 and 'total_nu' not in g:
tgroup = g.create_group('total_nu')
@ -673,7 +781,7 @@ class IncidentNeutron(EqualityMixin):
atomic_weight_ratio, temperature)
if (2, 151) in ev.section:
data.resonances = Resonances.from_endf(ev)
data.resonances = res.Resonances.from_endf(ev)
# Read each reaction
for mf, mt, nc, mod in ev.reaction_list:
@ -682,7 +790,7 @@ class IncidentNeutron(EqualityMixin):
# Replace cross sections for elastic, capture, fission
try:
if any(isinstance(r, _RESOLVED) for r in data.resonances):
if any(isinstance(r, res._RESOLVED) for r in data.resonances):
for mt in (2, 102, 18):
if mt in data.reactions:
rx = data.reactions[mt]

173
openmc/settings.py
View file

@ -11,8 +11,8 @@ from openmc.clean_xml import clean_xml_indentation
import openmc.checkvalue as cv
from openmc import Nuclide, VolumeCalculation, Source, Mesh
_RUN_MODES = ['eigenvalue', 'fixed source', 'plot', 'volume',
'particle restart']
_RUN_MODES = ['eigenvalue', 'fixed source', 'plot', 'volume', 'particle restart']
_RES_SCAT_METHODS = ['dbrc', 'wcm', 'ares']
class Settings(object):
@ -80,8 +80,18 @@ class Settings(object):
Number of particles per generation
ptables : bool
Determine whether probability tables are used.
resonance_scattering : ResonanceScattering or iterable of ResonanceScattering
The elastic scattering model to use for resonant isotopes
resonance_scattering : dict
Settings for resonance elastic scattering. Accepted keys are 'enable'
(bool), 'method' (str), 'energy_min' (float), 'energy_max' (float), and
'nuclides' (list). The 'method' can be set to 'dbrc' (Doppler broadening
rejection correction), 'wcm' (weight correction method), and 'ares'
(accelerated resonance elastic scattering). If not specified, 'ares' is
the default method. The 'energy_min' and 'energy_max' values indicate
the minimum and maximum energies above and below which the resonance
elastic scattering method is to be applied. The 'nuclides' list
indicates what nuclides the method should be applied to. In its absence,
the method will be applied to all nuclides with 0 K elastic scattering
data present.
run_cmfd : bool
Indicate if coarse mesh finite difference acceleration is to be used
run_mode : {'eigenvalue', 'fixed source', 'plot', 'volume', 'particle restart'}
@ -215,8 +225,7 @@ class Settings(object):
self._dd_allow_leakage = False
self._dd_count_interactions = False
self._resonance_scattering = cv.CheckedList(
ResonanceScattering, 'resonance scattering models')
self._resonance_scattering = {}
self._volume_calculations = cv.CheckedList(
VolumeCalculation, 'volume calculations')
@ -778,10 +787,27 @@ class Settings(object):
@resonance_scattering.setter
def resonance_scattering(self, res):
if not isinstance(res, MutableSequence):
res = [res]
self._resonance_scattering = cv.CheckedList(
ResonanceScattering, 'resonance scattering models', res)
cv.check_type('resonance scattering settings', res, Mapping)
keys = ('enable', 'method', 'energy_min', 'energy_max', 'nuclides')
for key, value in res.items():
cv.check_value('resonance scattering dictionary key', key, keys)
if key == 'enable':
cv.check_type('resonance scattering enable', value, bool)
elif key == 'method':
cv.check_value('resonance scattering method', value,
_RES_SCAT_METHODS)
elif key == 'energy_min':
name = 'resonance scattering minimum energy'
cv.check_type(name, value, Real)
cv.check_greater_than(name, value, 0)
elif key == 'energy_max':
name = 'resonance scattering minimum energy'
cv.check_type(name, value, Real)
cv.check_greater_than(name, value, 0)
elif key == 'nuclides':
cv.check_type('resonance scattering nuclides', value,
Iterable, string_types)
self._resonance_scattering = res
@volume_calculations.setter
def volume_calculations(self, vol_calcs):
@ -1049,10 +1075,24 @@ class Settings(object):
subelement.text = str(self._dd_count_interactions).lower()
def _create_resonance_scattering_subelement(self, root):
if len(self.resonance_scattering) > 0:
res = self.resonance_scattering
if res:
elem = ET.SubElement(root, 'resonance_scattering')
for r in self.resonance_scattering:
elem.append(r.to_xml_element())
if 'enable' in res:
subelem = ET.SubElement(elem, 'enable')
subelem.text = str(res['enable']).lower()
if 'method' in res:
subelem = ET.SubElement(elem, 'method')
subelem.text = res['method']
if 'energy_min' in res:
subelem = ET.SubElement(elem, 'energy_min')
subelem.text = str(res['energy_min'])
if 'energy_max' in res:
subelem = ET.SubElement(elem, 'energy_max')
subelem.text = str(res['energy_max'])
if 'nuclides' in res:
subelem = ET.SubElement(elem, 'nuclides')
subelem.text = ' '.join(res['nuclides'])
def _create_create_fission_neutrons_subelement(self, root):
if self._create_fission_neutrons is not None:
@ -1113,110 +1153,3 @@ class Settings(object):
# Write the XML Tree to the settings.xml file
tree = ET.ElementTree(root_element)
tree.write(path, xml_declaration=True, encoding='utf-8', method="xml")
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.
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.
"""
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 method(self):
return self._method
@property
def E_min(self):
return self._E_min
@property
def E_max(self):
return self._E_max
@nuclide.setter
def nuclide(self, nuc):
cv.check_type('nuclide', nuc, Nuclide)
self._nuclide = nuc
@method.setter
def method(self, m):
cv.check_value('method', m, ('ARES', 'CXS', 'DBRC', 'WCM'))
self._method = m
@E_min.setter
def E_min(self, E):
cv.check_type('E_min', E, Real)
cv.check_greater_than('E_min', E, 0, True)
self._E_min = E
@E_max.setter
def E_max(self, E):
cv.check_type('E_max', E, Real)
cv.check_greater_than('E_max', E, 0, True)
self._E_max = E
def to_xml_element(self):
"""Return XML representation of the resonance scattering model
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing resonance scattering model
"""
scatterer = ET.Element("scatterer")
subelement = ET.SubElement(scatterer, 'nuclide')
subelement.text = self.nuclide.name
if self.method is not None:
subelement = ET.SubElement(scatterer, 'method')
subelement.text = self.method
if self.E_min is not None:
subelement = ET.SubElement(scatterer, 'E_min')
subelement.text = str(self.E_min)
if self.E_max is not None:
subelement = ET.SubElement(scatterer, 'E_max')
subelement.text = str(self.E_max)
return scatterer

7
src/constants.F90
View file

@ -269,6 +269,13 @@ module constants
TEMPERATURE_NEAREST = 1, &
TEMPERATURE_INTERPOLATION = 2
! Resonance elastic scattering methods
integer, parameter :: &
RES_SCAT_ARES = 1, &
RES_SCAT_DBRC = 2, &
RES_SCAT_WCM = 3, &
RES_SCAT_CXS = 4
! ============================================================================
! TALLY-RELATED CONSTANTS

8
src/cross_section.F90
View file

@ -837,15 +837,17 @@ contains
real(8) :: xs_out ! 0K xs at trial energy
integer :: i_grid ! index on nuclide energy grid
integer :: n_grid
real(8) :: f ! interp factor on nuclide energy grid
! Determine index on nuclide energy grid
n_grid = size(nuc % energy_0K)
if (E < nuc % energy_0K(1)) then
i_grid = 1
elseif (E > nuc % energy_0K(nuc % n_grid_0K)) then
i_grid = nuc % n_grid_0K - 1
elseif (E > nuc % energy_0K(n_grid)) then
i_grid = n_grid - 1
else
i_grid = binary_search(nuc % energy_0K, nuc % n_grid_0K, E)
i_grid = binary_search(nuc % energy_0K, n_grid, E)
end if
! check for rare case where two energy points are the same

20
src/global.F90
View file

@ -429,9 +429,11 @@ module global
! ============================================================================
! RESONANCE SCATTERING VARIABLES
logical :: treat_res_scat = .false. ! is resonance scattering treated?
integer :: n_res_scatterers_total = 0 ! total number of resonant scatterers
type(Nuclide0K), allocatable, target :: nuclides_0K(:) ! 0K nuclides info
logical :: res_scat_on = .false. ! is resonance scattering treated?
integer :: res_scat_method = RES_SCAT_ARES ! resonance scattering method
real(8) :: res_scat_energy_min = 0.01_8
real(8) :: res_scat_energy_max = 1000.0_8
character(10), allocatable :: res_scat_nuclides(:)
!$omp threadprivate(micro_xs, material_xs, fission_bank, n_bank, &
!$omp& trace, thread_id, current_work, matching_bins, &
@ -468,17 +470,11 @@ contains
deallocate(nuclides)
end if
if (allocated(nuclides_0K)) then
deallocate(nuclides_0K)
end if
if (allocated(res_scat_nuclides)) deallocate(res_scat_nuclides)
if (allocated(nuclides_MG)) then
deallocate(nuclides_MG)
end if
if (allocated(nuclides_MG)) deallocate(nuclides_MG)
if (allocated(macro_xs)) then
deallocate(macro_xs)
end if
if (allocated(macro_xs)) deallocate(macro_xs)
if (allocated(sab_tables)) deallocate(sab_tables)
if (allocated(micro_xs)) deallocate(micro_xs)

187
src/input_xml.F90
View file

@ -94,11 +94,9 @@ contains
type(XMLNode) :: node_sp
type(XMLNode) :: node_output
type(XMLNode) :: node_res_scat
type(XMLNode) :: node_scatterer
type(XMLNode) :: node_trigger
type(XMLNode) :: node_vol
type(XMLNode) :: node_tab_leg
type(XMLNode), allocatable :: node_scat_list(:)
type(XMLNode), allocatable :: node_source_list(:)
type(XMLNode), allocatable :: node_vol_list(:)
@ -851,59 +849,54 @@ contains
! Resonance scattering parameters
if (check_for_node(root, "resonance_scattering")) then
node_res_scat = root % child("resonance_scattering")
call get_node_list(node_res_scat, "scatterer", node_scat_list)
! check that a nuclide is specified
if (size(node_scat_list) >= 1) then
treat_res_scat = .true.
n_res_scatterers_total = size(node_scat_list)
! store 0K info for resonant scatterers
allocate(nuclides_0K(n_res_scatterers_total))
do i = 1, n_res_scatterers_total
node_scatterer = node_scat_list(i)
! check to make sure a nuclide is specified
if (.not. check_for_node(node_scatterer, "nuclide")) then
call fatal_error("No nuclide specified for scatterer " &
// trim(to_str(i)) // " in settings.xml file!")
end if
call get_node_value(node_scatterer, "nuclide", &
nuclides_0K(i) % nuclide)
if (check_for_node(node_scatterer, "method")) then
call get_node_value(node_scatterer, "method", &
nuclides_0K(i) % scheme)
end if
if (check_for_node(node_scatterer, "E_min")) then
call get_node_value(node_scatterer, "E_min", &
nuclides_0K(i) % E_min)
end if
! check that E_min is non-negative
if (nuclides_0K(i) % E_min < ZERO) then
call fatal_error("Lower resonance scattering energy bound is &
&negative")
end if
if (check_for_node(node_scatterer, "E_max")) then
call get_node_value(node_scatterer, "E_max", &
nuclides_0K(i) % E_max)
end if
! check that E_max is not less than E_min
if (nuclides_0K(i) % E_max < nuclides_0K(i) % E_min) then
call fatal_error("Lower resonance scattering energy bound exceeds &
&upper")
end if
nuclides_0K(i) % nuclide = trim(nuclides_0K(i) % nuclide)
nuclides_0K(i) % scheme = to_lower(trim(nuclides_0K(i) % scheme))
end do
! See if resonance scattering is enabled
if (check_for_node(node_res_scat, "enable")) then
call get_node_value(node_res_scat, "enable", res_scat_on)
else
call fatal_error("No resonant scatterers are specified within the &
&resonance_scattering element in settings.xml")
res_scat_on = .true.
end if
! Determine what method is used
if (check_for_node(node_res_scat, "method")) then
call get_node_value(node_res_scat, "method", temp_str)
select case(to_lower(temp_str))
case ('ares')
res_scat_method = RES_SCAT_ARES
case ('dbrc')
res_scat_method = RES_SCAT_DBRC
case ('wcm')
res_scat_method = RES_SCAT_WCM
case default
call fatal_error("Unrecognized resonance elastic scattering method: " &
// trim(temp_str) // ".")
end select
end if
! Minimum energy for resonance scattering
if (check_for_node(node_res_scat, "energy_min")) then
call get_node_value(node_res_scat, "energy_min", res_scat_energy_min)
end if
if (res_scat_energy_min < ZERO) then
call fatal_error("Lower resonance scattering energy bound is negative")
end if
! Maximum energy for resonance scattering
if (check_for_node(node_res_scat, "energy_max")) then
call get_node_value(node_res_scat, "energy_max", res_scat_energy_max)
end if
if (res_scat_energy_max < res_scat_energy_min) then
call fatal_error("Upper resonance scattering energy bound is below the &
&lower resonance scattering energy bound.")
end if
! Get nuclides that resonance scattering should be applied to
if (check_for_node(node_res_scat, "nuclides")) then
n = node_word_count(node_res_scat, "nuclides")
allocate(res_scat_nuclides(n))
if (n > 0) then
call get_node_array(node_res_scat, "nuclides", res_scat_nuclides)
end if
end if
end if
@ -2147,12 +2140,11 @@ contains
end do
! Check that 0K nuclides are listed in the cross_sections.xml file
if (allocated(nuclides_0K)) then
do i = 1, size(nuclides_0K)
if (.not. library_dict % has_key(to_lower(nuclides_0K(i) % nuclide))) then
if (allocated(res_scat_nuclides)) then
do i = 1, size(res_scat_nuclides)
if (.not. library_dict % has_key(to_lower(res_scat_nuclides(i)))) then
call fatal_error("Could not find resonant scatterer " &
// trim(nuclides_0K(i) % nuclide) &
// " in cross_sections.xml file!")
// trim(res_scat_nuclides(i)) // " in cross_sections.xml file!")
end if
end do
end if
@ -5150,8 +5142,7 @@ contains
call file_close(file_id)
! Assign resonant scattering data
if (treat_res_scat) call read_0K_elastic_scattering(&
nuclides(i_nuclide), libraries, library_dict)
if (res_scat_on) call assign_0K_elastic_scattering(nuclides(i_nuclide))
! Determine if minimum/maximum energy for this nuclide is greater/less
! than the previous
@ -5280,79 +5271,43 @@ contains
end subroutine assign_temperatures
!===============================================================================
! READ_0K_ELASTIC_SCATTERING
! ASSIGN_0K_ELASTIC_SCATTERING
!===============================================================================
subroutine read_0K_elastic_scattering(nuc, libraries, library_dict)
type(Nuclide), intent(inout) :: nuc
type(Library), intent(in) :: libraries(:)
type(DictCharInt), intent(inout) :: library_dict
subroutine assign_0K_elastic_scattering(nuc)
type(Nuclide), intent(inout) :: nuc
integer :: i, j
integer :: i_library
integer :: method
integer(HID_T) :: file_id
integer(HID_T) :: group_id
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
! Copy basic information from settings.xml
do i = 1, size(res_scat_nuclides)
if (nuc % name == res_scat_nuclides(i)) then
! Set nuclide to be resonant
nuc % resonant = .true.
nuc % scheme = trim(nuclides_0K(i) % scheme)
nuc % E_min = nuclides_0K(i) % E_min
nuc % E_max = nuclides_0K(i) % E_max
! Get index in libraries array
name = nuc % name
i_library = library_dict % get_key(to_lower(name))
call write_message('Reading ' // trim(name) // ' 0K data from ' // &
trim(libraries(i_library) % path), 6)
! Open file and make sure version matches
file_id = file_open(libraries(i_library) % path, 'r')
! Read nuclide data from HDF5
group_id = open_group(file_id, name)
method = TEMPERATURE_NEAREST
call resonant_nuc % from_hdf5(group_id, temperature, &
method, 1000.0_8, master)
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 % 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
xs_cdf_sum = ZERO
allocate(nuc % xs_cdf(size(nuc % energy_0K)))
allocate(nuc % xs_cdf(0:size(nuc % energy_0K)))
nuc % xs_cdf(0) = ZERO
do j = 1, size(nuc % energy_0K) - 1
! Negative cross sections result in a CDF that is not monotonically
! increasing. Set all negative xs values to ZERO.
if (nuc % elastic_0K(j) < ZERO) nuc % elastic_0K(j) = ZERO
associate (E => nuc % energy_0K, xs => nuc % elastic_0K)
do j = 1, size(E) - 1
! Negative cross sections result in a CDF that is not monotonically
! increasing. Set all negative xs values to zero.
if (xs(j) < ZERO) xs(j) = ZERO
! build xs cdf
xs_cdf_sum = xs_cdf_sum &
+ (sqrt(nuc % energy_0K(j)) * nuc % elastic_0K(j) &
+ sqrt(nuc % energy_0K(j+1)) * nuc % elastic_0K(j+1)) / TWO &
* (nuc % energy_0K(j+1) - nuc % energy_0K(j))
nuc % xs_cdf(j) = xs_cdf_sum
end do
! build xs cdf
xs_cdf_sum = xs_cdf_sum + (sqrt(E(j))*xs(j) + sqrt(E(j+1))*xs(j+1))&
/ TWO * (E(j+1) - E(j))
nuc % xs_cdf(j) = xs_cdf_sum
end do
end associate
exit
end if
end do
end subroutine read_0K_elastic_scattering
end subroutine assign_0K_elastic_scattering
!===============================================================================
! READ_MULTIPOLE_DATA checks for the existence of a multipole library in the

52
src/nuclide_header.F90
View file

@ -65,14 +65,9 @@ module nuclide_header
! Resonance scattering info
logical :: resonant = .false. ! resonant scatterer?
character(10) :: name_0K = '' ! name of 0K nuclide, e.g. 92235.00c
character(16) :: scheme ! target velocity sampling scheme
integer :: n_grid_0K ! number of 0K energy grid points
real(8), allocatable :: energy_0K(:) ! energy grid for 0K xs
real(8), allocatable :: elastic_0K(:) ! Microscopic elastic cross section
real(8), allocatable :: xs_cdf(:) ! CDF of v_rel times cross section
real(8) :: E_min ! lower cutoff energy for res scattering
real(8) :: E_max ! upper cutoff energy for res scattering
! Fission information
logical :: has_partial_fission = .false. ! nuclide has partial fission reactions?
@ -106,18 +101,6 @@ module nuclide_header
procedure, private :: create_derived => nuclide_create_derived
end type Nuclide
!===============================================================================
! NUCLIDE0K temporarily contains all 0K cross section data and other parameters
! needed to treat resonance scattering before transferring them to Nuclide
!===============================================================================
type Nuclide0K
character(10) :: nuclide ! name of nuclide, e.g. U238
character(16) :: scheme = 'ares' ! target velocity sampling scheme
real(8) :: E_min = 0.01_8 ! lower cutoff energy for res scattering
real(8) :: E_max = 1000.0_8 ! upper cutoff energy for res scattering
end type Nuclide0K
!===============================================================================
! NUCLIDEMICROXS contains cached microscopic cross sections for a
! particular nuclide at the current energy
@ -201,6 +184,7 @@ module nuclide_header
integer(HID_T) :: kT_group
integer(HID_T) :: rxs_group
integer(HID_T) :: rx_group
integer(HID_T) :: xs, temp_group
integer(HID_T) :: total_nu
integer(HID_T) :: fer_group ! fission_energy_release group
integer(HID_T) :: fer_dset
@ -309,24 +293,36 @@ module nuclide_header
allocate(this % kTs(n_temperature))
allocate(this % grid(n_temperature))
! Get kT values
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))
end do
call close_group(kT_group)
! Read energy grid
energy_group = open_group(group_id, 'energy')
! Read energy grid
energy_group = open_group(group_id, 'energy')
do i = 1, n_temperature
temp_str = trim(to_str(temps_to_read % data(i))) // "K"
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
call close_group(kT_group)
! Check for 0K energy grid
if (object_exists(energy_group, '0K')) then
energy_dset = open_dataset(energy_group, '0K')
call get_shape(energy_dset, dims)
allocate(this % energy_0K(int(dims(1), 4)))
call read_dataset(this % energy_0K, energy_dset)
call close_dataset(energy_dset)
end if
call close_group(energy_group)
! Get MT values based on group names
rxs_group = open_group(group_id, 'reactions')
@ -344,6 +340,20 @@ module nuclide_header
zero_padded(MTs % data(i), 3)))
call this % reactions(i) % from_hdf5(rx_group, temps_to_read)
! Check for 0K elastic scattering
if (this % reactions(i) % MT == 2) then
if (object_exists(rx_group, '0K')) then
temp_group = open_group(rx_group, '0K')
xs = open_dataset(temp_group, 'xs')
call get_shape(xs, dims)
allocate(this % elastic_0K(int(dims(1), 4)))
call read_dataset(this % elastic_0K, xs)
call close_dataset(xs)
call close_group(temp_group)
end if
end if
call close_group(rx_group)
end do
call close_group(rxs_group)

221
src/physics.F90
View file

@ -809,46 +809,41 @@ contains
real(8) :: E_red ! reduced energy (same as used by Cullen in SIGMA1)
real(8) :: E_low ! lowest practical relative energy
real(8) :: E_up ! highest practical relative energy
real(8) :: E_mode ! most probable Maxwellian energy
real(8) :: E_t_max ! highest practical target energy
real(8) :: E_t ! trial target energy
real(8) :: xs_max ! max 0K xs over practical relative energies
real(8) :: xs_low ! 0K xs at lowest practical relative energy
real(8) :: xs_up ! 0K xs at highest practical relative energy
real(8) :: m ! slope for interpolation
real(8) :: R_dbrc ! DBRC rejection criterion
real(8) :: R_speed ! target speed rejection criterion
real(8) :: xi ! pseudorandom number on [0,1)
real(8) :: R ! rejection criterion for DBRC / target speed
real(8) :: cdf_low ! xs cdf at lowest practical relative energy
real(8) :: cdf_up ! xs cdf at highest practical relative energy
real(8) :: cdf_rel ! trial xs cdf value
real(8) :: p_mode ! probability at most probable energy
real(8) :: p_t ! probability at trial target energy
real(8) :: mu ! cosine between neutron and target velocities
integer :: i_E_low ! 0K index to lowest practical relative energy
integer :: i_E_up ! 0K index to highest practical relative energy
integer :: i_E_rel ! index to trial relative energy
integer :: n_grid ! number of energies on 0K grid
logical :: reject ! resample if true
character(80) :: sampling_scheme ! method of target velocity sampling
integer :: sampling_method ! method of target velocity sampling
awr = nuc % awr
! check if nuclide is a resonant scatterer
if (nuc % resonant) then
! sampling scheme to use
sampling_scheme = nuc % scheme
! sampling method to use
sampling_method = res_scat_method
! upper resonance scattering energy bound (target is at rest above this E)
if (E > nuc % E_max) then
if (E > res_scat_energy_max) then
v_target = ZERO
return
! lower resonance scattering energy bound (should be no resonances below)
else if (E < nuc % E_min) then
sampling_scheme = 'cxs'
else if (E < res_scat_energy_min) then
sampling_method = RES_SCAT_CXS
end if
! otherwise, use free gas model
@ -857,19 +852,18 @@ contains
v_target = ZERO
return
else
sampling_scheme = 'cxs'
sampling_method = RES_SCAT_CXS
end if
end if
! use appropriate target velocity sampling method
select case (sampling_scheme)
case ('cxs')
select case (sampling_method)
case (RES_SCAT_CXS)
! sample target velocity with the constant cross section (cxs) approx.
call sample_cxs_target_velocity(nuc, v_target, E, uvw, kT)
case ('wcm')
case (RES_SCAT_WCM)
! sample target velocity with the constant cross section (cxs) approx.
call sample_cxs_target_velocity(nuc, v_target, E, uvw, kT)
@ -880,151 +874,110 @@ contains
wcf = xs_0K / xs_eff
wgt = wcf * wgt
case ('dbrc')
case (RES_SCAT_DBRC, RES_SCAT_ARES)
E_red = sqrt((awr * E) / kT)
E_low = (((E_red - FOUR)**2) * kT) / awr
E_up = (((E_red + FOUR)**2) * kT) / awr
! find lower and upper energy bound indices
! lower index
n_grid = size(nuc % energy_0K)
if (E_low < nuc % energy_0K(1)) then
i_E_low = 1
elseif (E_low > nuc % energy_0K(nuc % n_grid_0K)) then
i_E_low = nuc % n_grid_0K - 1
elseif (E_low > nuc % energy_0K(n_grid)) then
i_E_low = n_grid - 1
else
i_E_low = binary_search(nuc % energy_0K, nuc % n_grid_0K, E_low)
i_E_low = binary_search(nuc % energy_0K, n_grid, E_low)
end if
! upper index
if (E_up < nuc % energy_0K(1)) then
i_E_up = 1
elseif (E_up > nuc % energy_0K(nuc % n_grid_0K)) then
i_E_up = nuc % n_grid_0K - 1
elseif (E_up > nuc % energy_0K(n_grid)) then
i_E_up = n_grid - 1
else
i_E_up = binary_search(nuc % energy_0K, nuc % n_grid_0K, E_up)
i_E_up = binary_search(nuc % energy_0K, n_grid, E_up)
end if
! interpolate xs since we're not exactly at the energy indices
xs_low = nuc % elastic_0K(i_E_low)
m = (nuc % elastic_0K(i_E_low + 1) - xs_low) &
& / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
xs_low = xs_low + m * (E_low - nuc % energy_0K(i_E_low))
xs_up = nuc % elastic_0K(i_E_up)
m = (nuc % elastic_0K(i_E_up + 1) - xs_up) &
& / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up))
xs_up = xs_up + m * (E_up - nuc % energy_0K(i_E_up))
! get max 0K xs value over range of practical relative energies
xs_max = max(xs_low, &
& maxval(nuc % elastic_0K(i_E_low + 1 : i_E_up - 1)), xs_up)
reject = .true.
! sample target velocities until one is accepted by the DBRC
do
! sample target velocity with the constant cross section (cxs) approx.
if (i_E_up == i_E_low) then
! Handle degenerate case -- if the upper/lower bounds occur for the same
! index, then using cxs is probably a good approximation
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))
xs_0K = elastic_xs_0K(E_rel, nuc)
R_dbrc = xs_0K / xs_max
if (prn() < R_dbrc) reject = .false.
if (.not. reject) exit
end do
case ('ares')
E_red = sqrt((awr * E) / kT)
E_low = (((E_red - FOUR)**2) * kT) / awr
E_up = (((E_red + FOUR)**2) * kT) / awr
! find lower and upper energy bound indices
! lower index
if (E_low < nuc % energy_0K(1)) then
i_E_low = 1
elseif (E_low > nuc % energy_0K(nuc % n_grid_0K)) then
i_E_low = nuc % n_grid_0K - 1
else
i_E_low = binary_search(nuc % energy_0K, nuc % n_grid_0K, E_low)
end if
if (sampling_method == RES_SCAT_DBRC) then
! interpolate xs since we're not exactly at the energy indices
xs_low = nuc % elastic_0K(i_E_low)
m = (nuc % elastic_0K(i_E_low + 1) - xs_low) &
/ (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
xs_low = xs_low + m * (E_low - nuc % energy_0K(i_E_low))
xs_up = nuc % elastic_0K(i_E_up)
m = (nuc % elastic_0K(i_E_up + 1) - xs_up) &
/ (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up))
xs_up = xs_up + m * (E_up - nuc % energy_0K(i_E_up))
! upper index
if (E_up < nuc % energy_0K(1)) then
i_E_up = 1
elseif (E_up > nuc % energy_0K(nuc % n_grid_0K)) then
i_E_up = nuc % n_grid_0K - 1
else
i_E_up = binary_search(nuc % energy_0K, nuc % n_grid_0K, E_up)
end if
! get max 0K xs value over range of practical relative energies
xs_max = max(xs_low, &
maxval(nuc % elastic_0K(i_E_low + 1 : i_E_up)), xs_up)
! interpolate xs CDF since we're not exactly at the energy indices
! cdf value at lower bound attainable energy
if (i_E_low > 1) then
m = (nuc % xs_cdf(i_E_low) - nuc % xs_cdf(i_E_low - 1)) &
& / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
cdf_low = nuc % xs_cdf(i_E_low - 1) &
& + m * (E_low - nuc % energy_0K(i_E_low))
else
m = nuc % xs_cdf(i_E_low) &
& / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
cdf_low = m * (E_low - nuc % energy_0K(i_E_low))
if (E_low <= nuc % energy_0K(1)) cdf_low = ZERO
end if
DBRC_REJECT_LOOP: do
! sample target velocity with the constant cross section (cxs) approx.
call sample_cxs_target_velocity(nuc, v_target, E, uvw, kT)
! cdf value at upper bound attainable energy
m = (nuc % xs_cdf(i_E_up) - nuc % xs_cdf(i_E_up - 1)) &
& / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up))
cdf_up = nuc % xs_cdf(i_E_up - 1) &
& + m * (E_up - nuc % energy_0K(i_E_up))
! perform Doppler broadening rejection correction (dbrc)
E_rel = dot_product((v_neut - v_target), (v_neut - v_target))
xs_0K = elastic_xs_0K(E_rel, nuc)
R = xs_0K / xs_max
if (prn() < R) exit DBRC_REJECT_LOOP
end do DBRC_REJECT_LOOP
! values used to sample the Maxwellian
E_mode = kT
p_mode = TWO * sqrt(E_mode / pi) * sqrt((ONE / kT)**3) &
& * exp(-E_mode / kT)
E_t_max = 16.0_8 * E_mode
elseif (sampling_method == RES_SCAT_ARES) then
! interpolate xs CDF since we're not exactly at the energy indices
! cdf value at lower bound attainable energy
m = (nuc % xs_cdf(i_E_low) - nuc % xs_cdf(i_E_low - 1)) &
/ (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
cdf_low = nuc % xs_cdf(i_E_low - 1) &
+ m * (E_low - nuc % energy_0K(i_E_low))
if (E_low <= nuc % energy_0K(1)) cdf_low = ZERO
reject = .true.
! cdf value at upper bound attainable energy
m = (nuc % xs_cdf(i_E_up) - nuc % xs_cdf(i_E_up - 1)) &
/ (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up))
cdf_up = nuc % xs_cdf(i_E_up - 1) &
+ m * (E_up - nuc % energy_0K(i_E_up))
do
ARES_REJECT_LOOP: do
! perform Maxwellian rejection sampling
xi = prn()
E_t = 16.0_8 * kT * xi**2
R = FOUR * xi * exp(ONE - E_t/kT)
! perform Maxwellian rejection sampling
E_t = E_t_max * prn()**2
p_t = TWO * sqrt(E_t / pi) * sqrt((ONE / kT)**3) &
& * exp(-E_t / kT)
R_speed = p_t / p_mode
if (prn() < R) then
! sample a relative energy using the xs cdf
cdf_rel = cdf_low + prn() * (cdf_up - cdf_low)
i_E_rel = binary_search(nuc % xs_cdf(i_E_low-1:i_E_up), &
i_E_up - i_E_low + 2, cdf_rel)
E_rel = nuc % energy_0K(i_E_low + i_E_rel - 1)
m = (nuc % xs_cdf(i_E_low + i_E_rel - 1) &
- nuc % xs_cdf(i_E_low + i_E_rel - 2)) &
/ (nuc % energy_0K(i_E_low + i_E_rel) &
- nuc % energy_0K(i_E_low + i_E_rel - 1))
E_rel = E_rel + (cdf_rel - nuc % xs_cdf(i_E_low + i_E_rel - 2)) / m
if (prn() < R_speed) then
! perform rejection sampling on cosine between
! neutron and target velocities
mu = (E_t + awr * (E - E_rel)) / (TWO * sqrt(awr * E * E_t))
! sample a relative energy using the xs cdf
cdf_rel = cdf_low + prn() * (cdf_up - cdf_low)
i_E_rel = binary_search(nuc % xs_cdf(i_E_low-1:i_E_up), &
& i_E_up - i_E_low + 2, cdf_rel)
E_rel = nuc % energy_0K(i_E_low + i_E_rel - 1)
m = (nuc % xs_cdf(i_E_low + i_E_rel - 1) &
& - nuc % xs_cdf(i_E_low + i_E_rel - 2)) &
& / (nuc % energy_0K(i_E_low + i_E_rel) &
& - nuc % energy_0K(i_E_low + i_E_rel - 1))
E_rel = E_rel + (cdf_rel - nuc % xs_cdf(i_E_low + i_E_rel - 2)) / m
! perform rejection sampling on cosine between
! neutron and target velocities
mu = (E_t + awr * (E - E_rel)) / (TWO * sqrt(awr * E * E_t))
if (abs(mu) < ONE) then
! set and accept target velocity
E_t = E_t / awr
v_target = sqrt(E_t) * rotate_angle(uvw, mu)
reject = .false.
end if
if (abs(mu) < ONE) then
! set and accept target velocity
E_t = E_t / awr
v_target = sqrt(E_t) * rotate_angle(uvw, mu)
exit ARES_REJECT_LOOP
end if
end if
end do ARES_REJECT_LOOP
end if
if (.not. reject) exit
end do
case default
call fatal_error("Not a recognized resonance scattering treatment!")
end if
end select
end subroutine sample_target_velocity

16
src/relaxng/settings.rnc
View file

@ -158,15 +158,11 @@ element settings {
}? &
element resonance_scattering {
element scatterer {
(element nuclide { xsd:string { maxLength = "12" } } |
attribute nuclide { xsd:string { maxLength = "12" } }) &
(element method { xsd:string { maxLength = "16" } } |
attribute method { xsd:string { maxLength = "16" } }) &
(element E_min { xsd:double } |
attribute E_min { xsd:double }) &
(element E_max { xsd:double } |
attribute E_max { xsd:double })?
}*
(element enable { xsd:boolean } | attribute enable { xsd:boolean })? &
(element method { xsd:string } | attribute method { xsd:string })? &
(element energy_min { xsd:double } | attribute energy_min { xsd:double })? &
(element energy_max { xsd:double } | attribute energy_max { xsd:double })? &
(element nuclides { list { xsd:string+ } } |
attribute nuclides { list { xsd:string+ } })?
}?
}

108
src/relaxng/settings.rng
View file

@ -715,54 +715,66 @@
</optional>
<optional>
<element name="resonance_scattering">
<zeroOrMore>
<element name="scatterer">
<interleave>
<choice>
<element name="nuclide">
<data type="string">
<param name="maxLength">12</param>
</data>
</element>
<attribute name="nuclide">
<data type="string">
<param name="maxLength">12</param>
</data>
</attribute>
</choice>
<choice>
<element name="method">
<data type="string">
<param name="maxLength">16</param>
</data>
</element>
<attribute name="method">
<data type="string">
<param name="maxLength">16</param>
</data>
</attribute>
</choice>
<choice>
<element name="E_min">
<data type="double"/>
</element>
<attribute name="E_min">
<data type="double"/>
</attribute>
</choice>
<optional>
<choice>
<element name="E_max">
<data type="double"/>
</element>
<attribute name="E_max">
<data type="double"/>
</attribute>
</choice>
</optional>
</interleave>
</element>
</zeroOrMore>
<interleave>
<optional>
<choice>
<element name="enable">
<data type="boolean"/>
</element>
<attribute name="enable">
<data type="boolean"/>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="method">
<data type="string"/>
</element>
<attribute name="method">
<data type="string"/>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="energy_min">
<data type="double"/>
</element>
<attribute name="energy_min">
<data type="double"/>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="energy_max">
<data type="double"/>
</element>
<attribute name="energy_max">
<data type="double"/>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="nuclides">
<list>
<oneOrMore>
<data type="string"/>
</oneOrMore>
</list>
</element>
<attribute name="nuclides">
<list>
<oneOrMore>
<data type="string"/>
</oneOrMore>
</list>
</attribute>
</choice>
</optional>
</interleave>
</element>
</optional>
</interleave>

2
src/summary.F90
View file

@ -111,7 +111,7 @@ contains
subroutine write_geometry(file_id)
integer(HID_T), intent(in) :: file_id
integer :: i, j, k, m, offset
integer :: i, j, k, m
integer, allocatable :: lattice_universes(:,:,:)
integer, allocatable :: cell_materials(:)
integer, allocatable :: cell_ids(:)

23
tests/test_resonance_scattering/inputs_true.dat
View file

@ -25,23 +25,10 @@
</space>
</source>
<resonance_scattering>
<scatterer>
<nuclide>U238</nuclide>
<method>DBRC</method>
<E_min>1.0</E_min>
<E_max>210.0</E_max>
</scatterer>
<scatterer>
<nuclide>U235</nuclide>
<method>WCM</method>
<E_min>1.0</E_min>
<E_max>210.0</E_max>
</scatterer>
<scatterer>
<nuclide>Pu239</nuclide>
<method>ARES</method>
<E_min>1.0</E_min>
<E_max>210.0</E_max>
</scatterer>
<enable>true</enable>
<method>ares</method>
<energy_min>1.0</energy_min>
<energy_max>210.0</energy_max>
<nuclides>U238 U235 Pu239</nuclides>
</resonance_scattering>
</settings>

2
tests/test_resonance_scattering/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
1.440556E+00 6.383274E-02
1.457296E+00 1.246018E-02

59
tests/test_resonance_scattering/test_resonance_scattering.py
View file

@ -9,54 +9,43 @@ 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.export_to_xml()
# Geometry
dumb_surface = openmc.XPlane(x0=100)
dumb_surface.boundary_type = 'reflective'
c1 = openmc.Cell(cell_id=1)
c1.fill = mat
c1.region = -dumb_surface
root_univ = openmc.Universe(universe_id=0)
root_univ.add_cell(c1)
geometry = openmc.Geometry()
geometry.root_universe = root_univ
dumb_surface = openmc.XPlane(x0=100, boundary_type='reflective')
c1 = openmc.Cell(cell_id=1, fill=mat, region=-dumb_surface)
root_univ = openmc.Universe(universe_id=0, cells=[c1])
geometry = openmc.Geometry(root_univ)
geometry.export_to_xml()
# Settings
res_scatt_dbrc = openmc.ResonanceScattering(u238, 'DBRC', 1.0, 210.0)
res_scatt_wcm = openmc.ResonanceScattering(u235, 'WCM', 1.0, 210.0)
res_scatt_ares = openmc.ResonanceScattering(pu239, 'ARES', 1.0, 210.0)
# Resonance elastic scattering settings
res_scat_settings = {
'enable': True,
'energy_min': 1.0,
'energy_max': 210.0,
'method': 'ares',
'nuclides': ['U238', 'U235', 'Pu239']
}
sets_file = openmc.Settings()
sets_file.batches = 10
sets_file.inactive = 5
sets_file.particles = 1000
sets_file.source = openmc.source.Source(
settings = openmc.Settings()
settings.batches = 10
settings.inactive = 5
settings.particles = 1000
settings.source = openmc.source.Source(
space=openmc.stats.Box([-4, -4, -4], [4, 4, 4]))
sets_file.resonance_scattering = [res_scatt_dbrc, res_scatt_wcm,
res_scatt_ares]
sets_file.export_to_xml()
settings.resonance_scattering = res_scat_settings
settings.export_to_xml()
if __name__ == '__main__':
harness = ResonanceScatteringTestHarness('statepoint.10.*')
harness = ResonanceScatteringTestHarness('statepoint.10.h5')
harness.main()