From 310fa134be0881f33aa6c91aa14584049dd62591 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 5 Dec 2016 16:51:55 -0500 Subject: [PATCH] Address #768 comments --- docs/source/usersguide/input.rst | 5 +- openmc/filter.py | 24 ++++--- openmc/tallies.py | 2 +- src/tally_filter.F90 | 62 ++++++++++++------- tests/test_filter_energyfun/inputs_true.dat | 2 +- .../test_filter_energyfun.py | 9 +-- 6 files changed, 60 insertions(+), 44 deletions(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index ea407e725b..09eab32621 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1639,8 +1639,9 @@ The ```` element accepts the following sub-elements: :energy: ``energyfunction`` filters multiply tally scores by an arbitrary function. The function is described by a piecewise linear-linear set of - (energy, y) values. This entry specifies the energy values. (Only used - for ``energyfunction`` filters) + (energy, y) values. This entry specifies the energy values. The function + will be evaluated as zero outside of the bounds of this energy grid. + (Only used for ``energyfunction`` filters) :y: ``energyfunction`` filters multiply tally scores by an arbitrary diff --git a/openmc/filter.py b/openmc/filter.py index ac6e961cf3..6bd7e26239 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1,6 +1,7 @@ from abc import ABCMeta from collections import Iterable, OrderedDict import copy +import hashlib from numbers import Real, Integral from xml.etree import ElementTree as ET @@ -158,8 +159,7 @@ class Filter(object): bins = np.array(bins) # If the bin is 0D numpy array, promote to 1D. - if bins.shape == (): - bins.shape = (1,) + bins = np.atleast_1d(bins) # Check the bin values. self.check_bins(bins) @@ -194,7 +194,7 @@ class Filter(object): pass - def to_xml(self): + def to_xml_element(self): """Return XML Element representing the Filter.""" element = ET.Element('filter') element.set('type', self.short_name.lower()) @@ -1616,7 +1616,8 @@ class EnergyFunctionFilter(Filter): """Multiplies tally scores by an arbitrary function of incident energy. The arbitrary function is described by a piecewise linear-linear - interpolation of energy and y values. + interpolation of energy and y values. Values outside of the given energy + range will be evaluated as zero. Parameters ---------- @@ -1723,8 +1724,7 @@ class EnergyFunctionFilter(Filter): energy = np.array(energy) # If the grid is 0D numpy array, promote to 1D. - if energy.shape == (): - energy.shape = (1,) + energy = np.atleast_1d(energy) # Make sure the values are Real and positive. cv.check_type('filter energy grid', energy, Iterable, Real) @@ -1739,8 +1739,7 @@ class EnergyFunctionFilter(Filter): y = np.array(y) # If the array is 0D, promote to 1D. - if y.shape == (): - y.shape = (1,) + y = np.atleast_1d(y) # Make sure the values are Real. cv.check_type('filter interpolant values', y, Iterable, Real) @@ -1751,7 +1750,7 @@ class EnergyFunctionFilter(Filter): def bins(self, bins): raise RuntimeError('EnergyFunctionFilters have no bins.') - def to_xml(self): + def to_xml_element(self): """Return XML Element representing the Filter.""" element = ET.Element('filter') element.set('type', self.short_name.lower()) @@ -1813,13 +1812,12 @@ class EnergyFunctionFilter(Filter): # and fill it with a hash of the __repr__. We want a hash that is # reproducible after restarting the interpreter so we'll use hashlib.md5 # rather than the intrinsic hash(). - from hashlib import md5 - hash_fun = md5() + hash_fun = hashlib.md5() hash_fun.update(repr(self).encode('utf-8')) out = hash_fun.hexdigest() - # The full 16 bytes make for a really wide column. Just 7 bytes of the - # digest are probably sufficient. + # The full 16 bytes make for a really wide column. Just 7 bytes (14 + # hex characters) of the digest are probably sufficient. out = out[:14] filter_bins = np.repeat(out, self.stride) diff --git a/openmc/tallies.py b/openmc/tallies.py index 35fb321751..a9a1390dfa 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -1036,7 +1036,7 @@ class Tally(object): # Optional Tally filters for self_filter in self.filters: - element.append(self_filter.to_xml()) + element.append(self_filter.to_xml_element()) # Optional Nuclides if len(self.nuclides) > 0: diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index f2680f4a9e..4b372defe0 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -1,7 +1,7 @@ module tally_filter use algorithm, only: binary_search - use constants, only: ONE, NO_BIN_FOUND, FP_PRECISION + use constants, only: ONE, NO_BIN_FOUND, FP_PRECISION, ERROR_REAL use dict_header, only: DictIntInt use geometry_header, only: BASE_UNIVERSE, RectLattice, HexLattice use global @@ -243,6 +243,8 @@ contains logical :: end_in_mesh ! ending coordinates inside mesh? type(RegularMesh), pointer :: m + weight = ERROR_REAL + ! Get a pointer to the mesh. m => meshes(this % mesh) @@ -251,10 +253,10 @@ contains ! valid mesh bin. if (current_bin == NO_BIN_FOUND) then call get_mesh_bin(m, p % coord(1) % xyz, next_bin) + weight = ONE else next_bin = NO_BIN_FOUND end if - weight = ONE else ! A track can span multiple mesh bins so we need to handle a lot of @@ -485,13 +487,14 @@ contains ! Starting one coordinate level deeper, find the next bin. next_bin = NO_BIN_FOUND + weight = ERROR_REAL do i = start, p % n_coord if (this % map % has_key(p % coord(i) % universe)) then next_bin = this % map % get_key(p % coord(i) % universe) + weight = ONE exit end if end do - weight = ONE end subroutine get_next_bin_universe subroutine to_statepoint_universe(this, filter_group) @@ -546,12 +549,13 @@ contains real(8), intent(out) :: weight next_bin = NO_BIN_FOUND + weight = ERROR_REAL if (current_bin == NO_BIN_FOUND) then if (this % map % has_key(p % material)) then next_bin = this % map % get_key(p % material) + weight = ONE end if end if - weight = ONE end subroutine get_next_bin_material subroutine to_statepoint_material(this, filter_group) @@ -621,13 +625,14 @@ contains ! Starting one coordinate level deeper, find the next bin. next_bin = NO_BIN_FOUND + weight = ERROR_REAL do i = start, p % n_coord if (this % map % has_key(p % coord(i) % cell)) then next_bin = this % map % get_key(p % coord(i) % cell) + weight = ONE exit end if end do - weight = ONE end subroutine get_next_bin_cell subroutine to_statepoint_cell(this, filter_group) @@ -681,10 +686,8 @@ contains integer, intent(out) :: next_bin real(8), intent(out) :: weight - logical :: cell_found integer :: distribcell_index, offset, i - cell_found = .false. if (current_bin == NO_BIN_FOUND) then distribcell_index = cells(this % cell) % distribcell_index offset = 0 @@ -707,15 +710,13 @@ contains end if if (this % cell == p % coord(i) % cell) then next_bin = offset + 1 - cell_found = .true. - exit + weight = ONE + return end if end do - else - next_bin = NO_BIN_FOUND end if - if (.not. cell_found) next_bin = NO_BIN_FOUND - weight = ONE + next_bin = NO_BIN_FOUND + weight = ERROR_REAL end subroutine get_next_bin_distribcell subroutine to_statepoint_distribcell(this, filter_group) @@ -770,12 +771,13 @@ contains real(8), intent(out) :: weight next_bin = NO_BIN_FOUND + weight = ERROR_REAL if (current_bin == NO_BIN_FOUND) then if (this % map % has_key(p % cell_born)) then next_bin = this % map % get_key(p % cell_born) + weight = ONE end if end if - weight = ONE end subroutine get_next_bin_cellborn subroutine to_statepoint_cellborn(this, filter_group) @@ -832,15 +834,16 @@ contains integer :: i next_bin = NO_BIN_FOUND + weight = ERROR_REAL if (current_bin == NO_BIN_FOUND) then do i = 1, this % n_bins if (p % surface == this % surfaces(i)) then next_bin = i + weight = ONE exit end if end do end if - weight = ONE end subroutine get_next_bin_surface subroutine to_statepoint_surface(this, filter_group) @@ -905,6 +908,7 @@ contains ! Tallies are ordered in increasing groups, group indices ! however are the opposite, so switch next_bin = num_energy_groups - next_bin + 1 + weight = ONE else ! Make sure the correct energy is used. @@ -917,16 +921,18 @@ contains ! Check if energy of the particle is within energy bins. if (E < this % bins(1) .or. E > this % bins(n + 1)) then next_bin = NO_BIN_FOUND + weight = ERROR_REAL else ! Search to find incoming energy bin. next_bin = binary_search(this % bins, n + 1, E) + weight = ONE end if end if else next_bin = NO_BIN_FOUND + weight = ERROR_REAL end if - weight = ONE end subroutine get_next_bin_energy subroutine to_statepoint_energy(this, filter_group) @@ -974,21 +980,24 @@ contains ! Tallies are ordered in increasing groups, group indices ! however are the opposite, so switch next_bin = num_energy_groups - next_bin + 1 + weight = ONE else ! Check if energy of the particle is within energy bins. if (p % E < this % bins(1) .or. p % E > this % bins(n + 1)) then next_bin = NO_BIN_FOUND + weight = ERROR_REAL else ! Search to find incoming energy bin. next_bin = binary_search(this % bins, n + 1, p % E) + weight = ONE end if end if else next_bin = NO_BIN_FOUND + weight = ERROR_REAL end if - weight = ONE end subroutine get_next_bin_energyout subroutine to_statepoint_energyout(this, filter_group) @@ -1026,10 +1035,11 @@ contains if (current_bin == NO_BIN_FOUND) then next_bin = 1 + weight = ONE else next_bin = NO_BIN_FOUND + weight = ERROR_REAL end if - weight = ONE end subroutine get_next_bin_dg subroutine to_statepoint_dg(this, filter_group) @@ -1068,15 +1078,17 @@ contains ! Check if energy of the particle is within energy bins. if (p % mu < this % bins(1) .or. p % mu > this % bins(n + 1)) then next_bin = NO_BIN_FOUND + weight = ERROR_REAL else ! Search to find incoming energy bin. next_bin = binary_search(this % bins, n + 1, p % mu) + weight = ONE end if else next_bin = NO_BIN_FOUND + weight = ERROR_REAL end if - weight = ONE end subroutine get_next_bin_mu subroutine to_statepoint_mu(this, filter_group) @@ -1129,15 +1141,17 @@ contains ! Check if particle is within polar angle bins. if (theta < this % bins(1) .or. theta > this % bins(n + 1)) then next_bin = NO_BIN_FOUND + weight = ERROR_REAL else ! Search to find polar angle bin. next_bin = binary_search(this % bins, n + 1, theta) + weight = ONE end if else next_bin = NO_BIN_FOUND + weight = ERROR_REAL end if - weight = ONE end subroutine get_next_bin_polar subroutine to_statepoint_polar(this, filter_group) @@ -1190,15 +1204,17 @@ contains ! Check if particle is within azimuthal angle bins. if (phi < this % bins(1) .or. phi > this % bins(n + 1)) then next_bin = NO_BIN_FOUND + weight = ERROR_REAL else ! Search to find azimuthal angle bin. next_bin = binary_search(this % bins, n + 1, phi) + weight = ONE end if else next_bin = NO_BIN_FOUND + weight = ERROR_REAL end if - weight = ONE end subroutine get_next_bin_azimuthal subroutine to_statepoint_azimuthal(this, filter_group) @@ -1253,7 +1269,7 @@ contains ! Check if energy of the particle is within energy bins. if (E < this % energy(1) .or. E > this % energy(n)) then next_bin = NO_BIN_FOUND - weight = ONE + weight = ERROR_REAL else ! Search to find incoming energy bin. @@ -1270,7 +1286,7 @@ contains else next_bin = NO_BIN_FOUND - weight = ONE + weight = ERROR_REAL end if end select end subroutine get_next_bin_energyfunction diff --git a/tests/test_filter_energyfun/inputs_true.dat b/tests/test_filter_energyfun/inputs_true.dat index 26f723392d..c9961071ae 100644 --- a/tests/test_filter_energyfun/inputs_true.dat +++ b/tests/test_filter_energyfun/inputs_true.dat @@ -1 +1 @@ -4dac56412e341cb64fa50d67e4dcf389b3b2cafb9bd8401eadc49654e050bad7988ca8239d5d9ee3e2953e1d36a2ca91a56557c7a34ccc72d520b41f462d9f1b \ No newline at end of file +2a08e03d9e0f40486fa26ba00800d264f43906a38c851da48f0fd76bacc0b862e8f45f16c970b7c62db03d96e3acb1b58e74179f554a851ccaaf4492e97fc163 \ No newline at end of file diff --git a/tests/test_filter_energyfun/test_filter_energyfun.py b/tests/test_filter_energyfun/test_filter_energyfun.py index 34141b038e..36dffb8718 100644 --- a/tests/test_filter_energyfun/test_filter_energyfun.py +++ b/tests/test_filter_energyfun/test_filter_energyfun.py @@ -18,15 +18,16 @@ class FilterEnergyFunHarness(PyAPITestHarness): # Add Am241 to the fuel. self._input_set.materials[1].add_nuclide('Am241', 1e-7) - # Make an EnergyFunctionFilter for the Am242m / Am242 branching ratio. + # Make an EnergyFunctionFilter for the Am242m / Am242 branching ratio + # (from ENDF/B-VII.1 data). x = [1e-5, 3.69e-1, 1e3, 1e5, 6e5, 1e6, 2e6, 4e6, 3e7] y = [0.1, 0.1, 0.1333, 0.158, 0.18467, 0.25618, 0.4297, 0.48, 0.48] filt = openmc.EnergyFunctionFilter(x, y) # Make tallies. - tallies = [openmc.Tally() for i in range(2)] + tallies = [openmc.Tally(), openmc.Tally()] for t in tallies: - t.scores = ['102'] + t.scores = ['(n,gamma)'] t.nuclides = ['Am241'] tallies[1].filters = [filt] self._input_set.tallies = openmc.Tallies(tallies) @@ -52,5 +53,5 @@ class FilterEnergyFunHarness(PyAPITestHarness): if __name__ == '__main__': - harness = FilterEnergyFunHarness('statepoint.10.*', True) + harness = FilterEnergyFunHarness('statepoint.10.h5', True) harness.main()