From 24c00acad403ab60a7328cba932f0eac5c691a20 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 19 Mar 2017 18:17:39 -0500 Subject: [PATCH 01/12] Allow 0K elastic scattering to be added to IncidentNeutron from ENDF files --- .travis.yml | 2 +- docs/source/pythonapi/index.rst | 2 + openmc/data/__init__.py | 1 + openmc/data/grid.py | 111 ++++++++++++++++++ openmc/data/neutron.py | 56 ++++++++- openmc/settings.py | 24 ++-- src/cross_section.F90 | 8 +- src/input_xml.F90 | 68 +++-------- src/nuclide_header.F90 | 36 +++++- src/physics.F90 | 27 +++-- src/summary.F90 | 2 +- .../results_true.dat | 2 +- .../test_resonance_scattering.py | 38 ++---- 13 files changed, 263 insertions(+), 114 deletions(-) create mode 100644 openmc/data/grid.py 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/pythonapi/index.rst b/docs/source/pythonapi/index.rst index c23bf1921..44e91abc2 100644 --- a/docs/source/pythonapi/index.rst +++ b/docs/source/pythonapi/index.rst @@ -371,6 +371,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/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/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..e24058ca8 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,12 +13,13 @@ 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 .resonance import RMatrixLimited, Unresolved, Resonances, _RESOLVED from .urr import ProbabilityTables import openmc.checkvalue as cv from openmc.mixin import EqualityMixin @@ -343,6 +345,46 @@ class IncidentNeutron(EqualityMixin): if strT in data.urr: self.urr[strT] = data.urr[strT] + def add_elastic_0K_from_endf(self, filename): + """Append 0K elastic scattering cross section from an ENDF file. + + Parameters + ---------- + filename : str + Path to ENDF file + + """ + data = type(self).from_endf(filename) + if data.resonances is not None: + x = [] + y = [] + for rr in data.resonances: + if isinstance(rr, RMatrixLimited): + raise TypeError('R-Matrix Limited not supported.') + elif isinstance(rr, Unresolved): + continue + + # Create 1000 equal log-spaced energies over RRR + e_min, e_max = rr.energy_min, rr.energy_max + energies = np.logspace(log10(e_min), log10(e_max), 1000) + + # 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 +453,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') diff --git a/openmc/settings.py b/openmc/settings.py index ff8c8a00f..e1b452a6a 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -1120,33 +1120,33 @@ class ResonanceScattering(object): Parameters ---------- - nuclide : openmc.Nuclide + nuclide : openmc.Nuclide or str 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. + The minimum energy in eV 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. + The maximum energy in eV below which the specified method is applied. + By default, the asymptotic target-at-rest model is applied above E_max. Attributes ---------- - nuclide : openmc.Nuclide + nuclide : openmc.Nuclide or str 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. + The minimum energy in eV 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. + The maximum energy in eV below which the specified method is applied. + By default, the asymptotic target-at-rest model is applied above E_max. """ @@ -1178,7 +1178,9 @@ class ResonanceScattering(object): @nuclide.setter def nuclide(self, nuc): - cv.check_type('nuclide', nuc, Nuclide) + cv.check_type('nuclide', nuc, (Nuclide,) + string_types) + if isinstance(nuc, string_types): + nuc = Nuclide(nuc) self._nuclide = nuc @method.setter 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/input_xml.F90 b/src/input_xml.F90 index 163dcf901..4866caf00 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -5150,8 +5150,8 @@ 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 (treat_res_scat) & + call assign_0K_elastic_scattering(nuclides(i_nuclide)) ! Determine if minimum/maximum energy for this nuclide is greater/less ! than the previous @@ -5280,25 +5280,14 @@ 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 @@ -5308,51 +5297,28 @@ contains 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))) - 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..147091390 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -65,9 +65,7 @@ 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 @@ -201,6 +199,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 +308,35 @@ 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 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, temp_str) + 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 +354,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..777be8197 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -828,6 +828,7 @@ contains 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 @@ -887,21 +888,22 @@ contains ! 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 @@ -941,21 +943,22 @@ contains ! 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 CDF since we're not exactly at the energy indices 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/results_true.dat b/tests/test_resonance_scattering/results_true.dat index e7056e4fa..28bd3e9cb 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.423955E+00 1.459996E-02 diff --git a/tests/test_resonance_scattering/test_resonance_scattering.py b/tests/test_resonance_scattering/test_resonance_scattering.py index 6689b17d5..0b80c5da7 100644 --- a/tests/test_resonance_scattering/test_resonance_scattering.py +++ b/tests/test_resonance_scattering/test_resonance_scattering.py @@ -9,42 +9,28 @@ 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) + 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) sets_file = openmc.Settings() sets_file.batches = 10 @@ -58,5 +44,5 @@ class ResonanceScatteringTestHarness(PyAPITestHarness): if __name__ == '__main__': - harness = ResonanceScatteringTestHarness('statepoint.10.*') + harness = ResonanceScatteringTestHarness('statepoint.10.h5') harness.main() From 8c7ca7afea9f4fd02dec48330a45fccae660c34f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 22 Mar 2017 09:37:18 -0500 Subject: [PATCH 02/12] Use resonance energies when performing reconstruction for 0 K elastic --- openmc/data/neutron.py | 21 +++++++++++++++++-- .../results_true.dat | 2 +- 2 files changed, 20 insertions(+), 3 deletions(-) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index e24058ca8..84969f6e2 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -345,15 +345,29 @@ class IncidentNeutron(EqualityMixin): if strT in data.urr: self.urr[strT] = data.urr[strT] - def add_elastic_0K_from_endf(self, filename): + 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 = [] @@ -366,7 +380,10 @@ class IncidentNeutron(EqualityMixin): # Create 1000 equal log-spaced energies over RRR e_min, e_max = rr.energy_min, rr.energy_max - energies = np.logspace(log10(e_min), log10(e_max), 1000) + e_log = np.logspace(log10(e_min), log10(e_max), 1000) + e_res = rr.parameters['energy'].values + in_range = (e_res > e_min) & (e_res < e_max) + energies = np.union1d(e_log, e_res[in_range]) # Linearize and thin cross section xi, yi = linearize(energies, data[2].xs['0K']) diff --git a/tests/test_resonance_scattering/results_true.dat b/tests/test_resonance_scattering/results_true.dat index 28bd3e9cb..f548f2f16 100644 --- a/tests/test_resonance_scattering/results_true.dat +++ b/tests/test_resonance_scattering/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.423955E+00 1.459996E-02 +1.399343E+00 1.149372E-01 From 474e11e87633b11875d09c68924f8828a35014a5 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 22 Mar 2017 10:59:36 -0500 Subject: [PATCH 03/12] Change how resonance scattering options are specified --- docs/source/conf.py | 3 +- docs/source/pythonapi/index.rst | 1 - docs/source/usersguide/input.rst | 84 +++++---- openmc/settings.py | 175 ++++++------------ src/constants.F90 | 7 + src/global.F90 | 20 +- src/input_xml.F90 | 119 ++++++------ src/nuclide_header.F90 | 15 -- src/physics.F90 | 28 ++- src/relaxng/settings.rnc | 16 +- src/relaxng/settings.rng | 108 ++++++----- .../test_resonance_scattering/inputs_true.dat | 23 +-- .../results_true.dat | 2 +- .../test_resonance_scattering.py | 27 +-- 14 files changed, 265 insertions(+), 363 deletions(-) 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 44e91abc2..3299458f1 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 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/settings.py b/openmc/settings.py index e1b452a6a..8f7338cb3 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, float) + cv.check_greater_than(name, value, 0) + elif key == 'energy_max': + name = 'resonance scattering minimum energy' + cv.check_type(name, value, float) + 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,112 +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 or str - 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 in eV above which the specified method is applied. - By default, CXS will be used below E_min. - E_max : float - The maximum energy in eV below which the specified method is applied. - By default, the asymptotic target-at-rest model is applied above E_max. - - Attributes - ---------- - nuclide : openmc.Nuclide or str - 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 in eV above which the specified method is applied. - By default, CXS will be used below E_min. - E_max : float - The maximum energy in eV 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,) + string_types) - if isinstance(nuc, string_types): - nuc = Nuclide(nuc) - 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/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 4866caf00..432990322 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,53 @@ 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 < ZERO) then + call fatal_error("Upper resonance scattering energy bound is negative") + 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 +2139,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 +5141,7 @@ contains call file_close(file_id) ! Assign resonant scattering data - if (treat_res_scat) & - call assign_0K_elastic_scattering(nuclides(i_nuclide)) + 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 @@ -5289,13 +5279,10 @@ contains integer :: i, j real(8) :: xs_cdf_sum - 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 ! Build CDF for 0K elastic scattering xs_cdf_sum = ZERO diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 147091390..b8a54f303 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -65,12 +65,9 @@ module nuclide_header ! Resonance scattering info logical :: resonant = .false. ! resonant scatterer? - character(16) :: scheme ! target velocity sampling scheme 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? @@ -104,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 diff --git a/src/physics.F90 b/src/physics.F90 index 777be8197..2a4afeaed 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -832,24 +832,24 @@ contains 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 @@ -858,19 +858,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) @@ -881,7 +880,7 @@ contains wcf = xs_0K / xs_eff wgt = wcf * wgt - case ('dbrc') + case (RES_SCAT_DBRC) E_red = sqrt((awr * E) / kT) E_low = (((E_red - FOUR)**2) * kT) / awr E_up = (((E_red + FOUR)**2) * kT) / awr @@ -936,7 +935,7 @@ contains if (.not. reject) exit end do - case ('ares') + case (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 @@ -1025,9 +1024,6 @@ contains if (.not. reject) exit end do - - case default - call fatal_error("Not a recognized resonance scattering treatment!") 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/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 f548f2f16..622e12cca 100644 --- a/tests/test_resonance_scattering/results_true.dat +++ b/tests/test_resonance_scattering/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.399343E+00 1.149372E-01 +1.439342E+00 1.224486E-02 diff --git a/tests/test_resonance_scattering/test_resonance_scattering.py b/tests/test_resonance_scattering/test_resonance_scattering.py index 0b80c5da7..94c9f5c3d 100644 --- a/tests/test_resonance_scattering/test_resonance_scattering.py +++ b/tests/test_resonance_scattering/test_resonance_scattering.py @@ -27,20 +27,23 @@ class ResonanceScatteringTestHarness(PyAPITestHarness): 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__': From 6f4fb29d657a9d0f0c39f68924f6a14dd8ac5a15 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 22 Mar 2017 12:51:11 -0500 Subject: [PATCH 04/12] Handle numbers in ACE files less than 1e-100 correctly --- openmc/data/ace.py | 13 ++++++++++++- openmc/data/endf.py | 4 ++-- 2 files changed, 14 insertions(+), 3 deletions(-) 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): From f49e0f9cf18494098c02eba48a0e23fc2b01e077 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 24 Mar 2017 14:53:25 -0500 Subject: [PATCH 05/12] When reconstructing resonances, add energies at half-height --- openmc/data/neutron.py | 36 ++++++++++++++++++++++++++---------- 1 file changed, 26 insertions(+), 10 deletions(-) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index 84969f6e2..67de69ff1 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 functools import reduce from math import log10 from numbers import Integral, Real from warnings import warn @@ -19,7 +20,7 @@ from .function import Tabulated1D, Sum, ResonancesWithBackground from .grid import linearize, thin from .product import Product from .reaction import Reaction, _get_photon_products_ace -from .resonance import RMatrixLimited, Unresolved, Resonances, _RESOLVED +from . import resonance as res from .urr import ProbabilityTables import openmc.checkvalue as cv from openmc.mixin import EqualityMixin @@ -280,7 +281,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 @@ -373,17 +374,32 @@ class IncidentNeutron(EqualityMixin): x = [] y = [] for rr in data.resonances: - if isinstance(rr, RMatrixLimited): + if isinstance(rr, res.RMatrixLimited): raise TypeError('R-Matrix Limited not supported.') - elif isinstance(rr, Unresolved): + elif isinstance(rr, res.Unresolved): continue - # Create 1000 equal log-spaced energies over RRR + # Get energies/widths for resonances + e_peak = rr.parameters['energy'].values + if isinstance(rr, res.MultiLevelBreitWigner): + width = rr.parameters['totalWidth'].values + elif isinstance(rr, res.ReichMoore): + df = rr.parameters + width = (df['neutronWidth'] + df['captureWidth'] + + df['fissionWidthA'] + df['fissionWidthB']).values + + # Determine energies at half-height 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] + e_half_left = e_peak - width[in_range]/2 + e_half_right = e_peak + width[in_range]/2 + + # 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) - e_res = rr.parameters['energy'].values - in_range = (e_res > e_min) & (e_res < e_max) - energies = np.union1d(e_log, e_res[in_range]) + energies = reduce(np.union1d, (e_log, e_peak, e_half_left, + e_half_right)) # Linearize and thin cross section xi, yi = linearize(energies, data[2].xs['0K']) @@ -742,7 +758,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: @@ -751,7 +767,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] From 60e4d6c539dc63ecf6137fb86867f513fe64f643 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 25 Mar 2017 16:58:37 -0500 Subject: [PATCH 06/12] Add energy checks suggested by @walshjon --- src/input_xml.F90 | 14 ++++++++++++-- tests/test_resonance_scattering/results_true.dat | 2 +- 2 files changed, 13 insertions(+), 3 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 432990322..514244dde 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -885,8 +885,9 @@ contains 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 < ZERO) then - call fatal_error("Upper resonance scattering energy bound is negative") + 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 @@ -5301,6 +5302,15 @@ contains end do end associate + ! Check to make sure resonance scattering upper energy is below URR + ! minimum, if present + if (nuc % urr_present) then + if (res_scat_energy_max > nuc % urr_data(1) % energy(1)) then + call fatal_error("Resonance scattering maximum energy is above the & + &bottom of the unresolved resonance region.") + end if + end if + exit end if end do diff --git a/tests/test_resonance_scattering/results_true.dat b/tests/test_resonance_scattering/results_true.dat index 622e12cca..0db0e4a95 100644 --- a/tests/test_resonance_scattering/results_true.dat +++ b/tests/test_resonance_scattering/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.439342E+00 1.224486E-02 +1.441041E+00 8.832991E-03 From cf7faabc82824d5a5d09c793e4401d2cb0667c59 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 27 Mar 2017 07:46:31 -0500 Subject: [PATCH 07/12] Respond to comments from @smharper on #840 --- openmc/settings.py | 4 ++-- src/input_xml.F90 | 9 --------- src/nuclide_header.F90 | 3 ++- 3 files changed, 4 insertions(+), 12 deletions(-) diff --git a/openmc/settings.py b/openmc/settings.py index 8f7338cb3..3a5bd63a2 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -798,11 +798,11 @@ class Settings(object): _RES_SCAT_METHODS) elif key == 'energy_min': name = 'resonance scattering minimum energy' - cv.check_type(name, value, float) + 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, float) + cv.check_type(name, value, Real) cv.check_greater_than(name, value, 0) elif key == 'nuclides': cv.check_type('resonance scattering nuclides', value, diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 514244dde..c42820bed 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -5302,15 +5302,6 @@ contains end do end associate - ! Check to make sure resonance scattering upper energy is below URR - ! minimum, if present - if (nuc % urr_present) then - if (res_scat_energy_max > nuc % urr_data(1) % energy(1)) then - call fatal_error("Resonance scattering maximum energy is above the & - &bottom of the unresolved resonance region.") - end if - end if - exit end if end do diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index b8a54f303..ed4a9a0a6 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -306,6 +306,7 @@ module nuclide_header ! 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))) @@ -315,7 +316,7 @@ module nuclide_header ! Check for 0K energy grid if (object_exists(energy_group, '0K')) then - energy_dset = open_dataset(energy_group, temp_str) + 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) From f15f740abf63a90b9760506ccba72859bd44e82d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 27 Mar 2017 22:37:18 -0500 Subject: [PATCH 08/12] Make sure all energies are positive when reconstructing resonances Apparently E_peak - gamma/2 can go negative (e.g. U235) --- openmc/data/neutron.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index 67de69ff1..b0528b475 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -395,6 +395,9 @@ class IncidentNeutron(EqualityMixin): e_half_left = e_peak - width[in_range]/2 e_half_right = e_peak + width[in_range]/2 + # Make sure all values are positive + e_half_left = e_half_left[e_half_left > 0.] + # 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) From 05440f5f2a48aa78f4fd8e3e732fff8167ce9e7b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 27 Mar 2017 22:38:56 -0500 Subject: [PATCH 09/12] Handle res scat case when i_E_low == i_E_up (assume cxs) --- src/physics.F90 | 197 +++++++++++++++++++++--------------------------- 1 file changed, 84 insertions(+), 113 deletions(-) diff --git a/src/physics.F90 b/src/physics.F90 index 2a4afeaed..d33ed1877 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -830,8 +830,6 @@ contains integer :: i_E_rel ! index to trial relative energy integer :: n_grid ! number of energies on 0K grid - logical :: reject ! resample if true - integer :: sampling_method ! method of target velocity sampling awr = nuc % awr @@ -880,7 +878,7 @@ contains wcf = xs_0K / xs_eff wgt = wcf * wgt - case (RES_SCAT_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 @@ -905,125 +903,98 @@ contains 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 (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(n_grid)) then - i_E_low = n_grid - 1 else - i_E_low = binary_search(nuc % energy_0K, n_grid, 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(n_grid)) then - i_E_up = n_grid - 1 - else - i_E_up = binary_search(nuc % energy_0K, n_grid, 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_dbrc = xs_0K / xs_max + if (prn() < R_dbrc) 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 - - reject = .true. - - do - - ! 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_speed) 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 - - ! 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. + 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 + 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 - end if - if (.not. reject) exit - end do + ! 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)) + + ! 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 + + ARES_REJECT_LOOP: do + ! 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_speed) 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 + + ! 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) + exit ARES_REJECT_LOOP + end if + end if + end do ARES_REJECT_LOOP + end if + end if end select end subroutine sample_target_velocity From a77c8eee4181fcae84d04133d4622f508b67ec27 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 29 Mar 2017 11:02:42 -0500 Subject: [PATCH 10/12] Extend xs_cdf to 0 to prevent potential out-of-bounds error --- src/input_xml.F90 | 3 ++- src/physics.F90 | 16 +++++----------- tests/test_resonance_scattering/results_true.dat | 2 +- 3 files changed, 8 insertions(+), 13 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index c42820bed..9bf78324d 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -5287,7 +5287,8 @@ contains ! 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 associate (E => nuc % energy_0K, xs => nuc % elastic_0K) do j = 1, size(E) - 1 diff --git a/src/physics.F90 b/src/physics.F90 index d33ed1877..2a1970a67 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -938,17 +938,11 @@ contains 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 - 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 + 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 ! cdf value at upper bound attainable energy m = (nuc % xs_cdf(i_E_up) - nuc % xs_cdf(i_E_up - 1)) & diff --git a/tests/test_resonance_scattering/results_true.dat b/tests/test_resonance_scattering/results_true.dat index 0db0e4a95..d1437a42a 100644 --- a/tests/test_resonance_scattering/results_true.dat +++ b/tests/test_resonance_scattering/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.441041E+00 8.832991E-03 +1.436406E+00 1.833812E-02 From d17ee5d54869e58eca7bc49abb11afd98e1f0417 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 29 Mar 2017 11:30:10 -0500 Subject: [PATCH 11/12] Simplify logic for Maxwellian sampling in ARES --- src/physics.F90 | 27 ++++++++------------------- 1 file changed, 8 insertions(+), 19 deletions(-) diff --git a/src/physics.F90 b/src/physics.F90 index 2a1970a67..a946c161b 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -809,20 +809,16 @@ 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 @@ -931,8 +927,8 @@ contains ! 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) exit DBRC_REJECT_LOOP + R = xs_0K / xs_max + if (prn() < R) exit DBRC_REJECT_LOOP end do DBRC_REJECT_LOOP elseif (sampling_method == RES_SCAT_ARES) then @@ -950,20 +946,13 @@ contains cdf_up = nuc % xs_cdf(i_E_up - 1) & + m * (E_up - nuc % energy_0K(i_E_up)) - ! 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 - ARES_REJECT_LOOP: do ! 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 + xi = prn() + E_t = 16.0_8 * kT * xi**2 + R = FOUR * xi * exp(ONE - E_t/kT) - if (prn() < R_speed) then + 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), & From a76d46de1822460be2532c717d69a290e55f720b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 30 Mar 2017 09:08:38 -0500 Subject: [PATCH 12/12] Adopt Fudge method for reconstructing resonances --- openmc/data/neutron.py | 42 ++++++++++++++----- .../results_true.dat | 2 +- 2 files changed, 32 insertions(+), 12 deletions(-) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index b0528b475..3449c7d58 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -2,7 +2,6 @@ from __future__ import division, unicode_literals import sys from collections import OrderedDict, Iterable, Mapping, MutableMapping from itertools import chain -from functools import reduce from math import log10 from numbers import Integral, Real from warnings import warn @@ -26,6 +25,10 @@ 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. @@ -382,27 +385,44 @@ class IncidentNeutron(EqualityMixin): # Get energies/widths for resonances e_peak = rr.parameters['energy'].values if isinstance(rr, res.MultiLevelBreitWigner): - width = rr.parameters['totalWidth'].values + gamma = rr.parameters['totalWidth'].values elif isinstance(rr, res.ReichMoore): df = rr.parameters - width = (df['neutronWidth'] + df['captureWidth'] + - df['fissionWidthA'] + df['fissionWidthB']).values + gamma = (df['neutronWidth'] + + df['captureWidth'] + + abs(df['fissionWidthA']) + + abs(df['fissionWidthB'])).values - # Determine energies at half-height + # 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] - e_half_left = e_peak - width[in_range]/2 - e_half_right = e_peak + width[in_range]/2 + gamma = gamma[in_range] - # Make sure all values are positive - e_half_left = e_half_left[e_half_left > 0.] + # 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 = reduce(np.union1d, (e_log, e_peak, e_half_left, - e_half_right)) + energies = np.union1d(e_log, energies) # Linearize and thin cross section xi, yi = linearize(energies, data[2].xs['0K']) diff --git a/tests/test_resonance_scattering/results_true.dat b/tests/test_resonance_scattering/results_true.dat index d1437a42a..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.436406E+00 1.833812E-02 +1.457296E+00 1.246018E-02