crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Address #768 comments

Browse source
This commit is contained in:
Sterling Harper 2016-12-05 16:51:55 -05:00
parent e8c9fe664a
commit 310fa134be
6 changed files with 60 additions and 44 deletions
Download patch file Download diff file Expand all files Collapse all files

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

@ -1639,8 +1639,9 @@ The ``<tally>`` 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

24
openmc/filter.py
View file

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

2
openmc/tallies.py
View file

@ -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:

62
src/tally_filter.F90
View file

@ -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

2
tests/test_filter_energyfun/inputs_true.dat
View file

@ -1 +1 @@
4dac56412e341cb64fa50d67e4dcf389b3b2cafb9bd8401eadc49654e050bad7988ca8239d5d9ee3e2953e1d36a2ca91a56557c7a34ccc72d520b41f462d9f1b
2a08e03d9e0f40486fa26ba00800d264f43906a38c851da48f0fd76bacc0b862e8f45f16c970b7c62db03d96e3acb1b58e74179f554a851ccaaf4492e97fc163

9
tests/test_filter_energyfun/test_filter_energyfun.py
View file

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