diff --git a/.travis.yml b/.travis.yml index e80ab4a9c..c1d3c86b2 100644 --- a/.travis.yml +++ b/.travis.yml @@ -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 diff --git a/docs/source/conf.py b/docs/source/conf.py index 89dd585d7..f6e342880 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -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 diff --git a/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst index 161b8564a..5ff2af2c2 100644 --- a/docs/source/pythonapi/index.rst +++ b/docs/source/pythonapi/index.rst @@ -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 diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index e1775d364..b4ba5350b 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -361,54 +361,56 @@ or sub-elements and can be set to either "false" or "true". ```` 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 ```` 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" - - - U-238 - ARES - 5.0e-6 - 40.0e-6 - - - Pu-239 - dbrc - 0.01e-6 - 210.0e-6 - - + :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 ```` is present but the ```` + 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`. diff --git a/openmc/data/__init__.py b/openmc/data/__init__.py index 9e7a44f0d..7158e9fe3 100644 --- a/openmc/data/__init__.py +++ b/openmc/data/__init__.py @@ -28,3 +28,4 @@ from .library import * from .fission_energy import * from .resonance import * from .multipole import * +from .grid import * diff --git a/openmc/data/ace.py b/openmc/data/ace.py index d7ec97edc..9827df5b9 100644 --- a/openmc/data/ace.py +++ b/openmc/data/ace.py @@ -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) diff --git a/openmc/data/endf.py b/openmc/data/endf.py index 34553ad2d..292d73ae3 100644 --- a/openmc/data/endf.py +++ b/openmc/data/endf.py @@ -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): diff --git a/openmc/data/grid.py b/openmc/data/grid.py new file mode 100644 index 000000000..f08bac999 --- /dev/null +++ b/openmc/data/grid.py @@ -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)] diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index 581adc417..3449c7d58 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -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] diff --git a/openmc/settings.py b/openmc/settings.py index ff8c8a00f..3a5bd63a2 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -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 diff --git a/src/constants.F90 b/src/constants.F90 index 5952aeaba..c44896b31 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -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 diff --git a/src/cross_section.F90 b/src/cross_section.F90 index b471ae54f..19e44538c 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -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 diff --git a/src/global.F90 b/src/global.F90 index c48367214..780f323f3 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -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) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 163dcf901..9bf78324d 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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 diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 571dd9023..ed4a9a0a6 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -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) diff --git a/src/physics.F90 b/src/physics.F90 index fceca20a9..a946c161b 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -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 diff --git a/src/relaxng/settings.rnc b/src/relaxng/settings.rnc index 1cac11e54..2043feb16 100644 --- a/src/relaxng/settings.rnc +++ b/src/relaxng/settings.rnc @@ -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+ } })? }? } diff --git a/src/relaxng/settings.rng b/src/relaxng/settings.rng index 8db4efcbc..f0c7ddbbf 100644 --- a/src/relaxng/settings.rng +++ b/src/relaxng/settings.rng @@ -715,54 +715,66 @@ - - - - - - - 12 - - - - - 12 - - - - - - - 16 - - - - - 16 - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/summary.F90 b/src/summary.F90 index 648042f68..9c952fe70 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -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(:) diff --git a/tests/test_resonance_scattering/inputs_true.dat b/tests/test_resonance_scattering/inputs_true.dat index 2801f115d..745fd30c3 100644 --- a/tests/test_resonance_scattering/inputs_true.dat +++ b/tests/test_resonance_scattering/inputs_true.dat @@ -25,23 +25,10 @@ - - U238 - DBRC - 1.0 - 210.0 - - - U235 - WCM - 1.0 - 210.0 - - - Pu239 - ARES - 1.0 - 210.0 - + true + ares + 1.0 + 210.0 + U238 U235 Pu239 diff --git a/tests/test_resonance_scattering/results_true.dat b/tests/test_resonance_scattering/results_true.dat index e7056e4fa..2dfadb7dd 100644 --- a/tests/test_resonance_scattering/results_true.dat +++ b/tests/test_resonance_scattering/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.440556E+00 6.383274E-02 +1.457296E+00 1.246018E-02 diff --git a/tests/test_resonance_scattering/test_resonance_scattering.py b/tests/test_resonance_scattering/test_resonance_scattering.py index 6689b17d5..94c9f5c3d 100644 --- a/tests/test_resonance_scattering/test_resonance_scattering.py +++ b/tests/test_resonance_scattering/test_resonance_scattering.py @@ -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()