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