Merge pull request #538 from smharper/distribcell

Move distribcell map indices from filters to cells
This commit is contained in:
Will Boyd 2016-01-04 19:18:16 -05:00
commit ffbc0b8a53
15 changed files with 109 additions and 115 deletions

View file

@ -185,10 +185,6 @@ if run_mode == 'k-eigenvalue':
Type of the j-th filter. Can be 'universe', 'material', 'cell', 'cellborn',
'surface', 'mesh', 'energy', 'energyout', or 'distribcell'.
**/tallies/tally <uid>/filter <j>/offset** (*int*)
Filter offset (used for distribcell filter).
**/tallies/tally <uid>/filter <j>/n_bins** (*int*)
Number of bins for the j-th filter.

View file

@ -121,6 +121,10 @@ The current revision of the summary file format is 1.
Region specification for the cell.
**/geometry/cells/cell <uid>/distribcell_index** (*int*)
Index of this cell in distribcell filter arrays.
**/geometry/surfaces/surface <uid>/index** (*int*)
Index in surfaces array used internally in OpenMC.

View file

@ -42,8 +42,6 @@ class Filter(object):
The number of filter bins
mesh : Mesh or None
A Mesh object for 'mesh' type filters.
offset : Integral
A value used to index tally bins for 'distribcell' tallies.
stride : Integral
The number of filter, nuclide and score bins within each of this
filter's bins.
@ -57,7 +55,6 @@ class Filter(object):
self._num_bins = 0
self._bins = None
self._mesh = None
self._offset = -1
self._stride = None
if type is not None:
@ -93,7 +90,6 @@ class Filter(object):
clone._bins = copy.deepcopy(self.bins, memo)
clone._num_bins = self.num_bins
clone._mesh = copy.deepcopy(self.mesh, memo)
clone._offset = self.offset
clone._stride = self.stride
memo[id(self)] = clone
@ -108,7 +104,6 @@ class Filter(object):
string = 'Filter\n'
string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self.type)
string += '{0: <16}{1}{2}\n'.format('\tBins', '=\t', self.bins)
string += '{0: <16}{1}{2}\n'.format('\tOffset', '=\t', self.offset)
return string
@property
@ -134,10 +129,6 @@ class Filter(object):
def mesh(self):
return self._mesh
@property
def offset(self):
return self._offset
@property
def stride(self):
return self._stride
@ -226,11 +217,6 @@ class Filter(object):
self.type = 'mesh'
self.bins = self.mesh.id
@offset.setter
def offset(self, offset):
cv.check_type('filter offset', offset, Integral)
self._offset = offset
@stride.setter
def stride(self, stride):
cv.check_type('filter stride', stride, Integral)
@ -623,7 +609,7 @@ class Filter(object):
# If this region is in Cell corresponding to the
# distribcell filter bin, store it in dictionary
if cell_id == self.bins[0]:
offset = openmc_geometry.get_offset(path, self.offset)
offset = openmc_geometry.get_cell_instance(path)
offsets_to_coords[offset] = coords
# Each distribcell offset is a DataFrame bin

View file

@ -42,10 +42,10 @@ class Geometry(object):
self._root_universe = root_universe
def get_offset(self, path, filter_offset):
"""Returns the corresponding location in the results array for a given path and
filter number. This is primarily intended to post-processing result when
a distribcell filter is used.
def get_cell_instance(self, path):
"""Return the instance number for the final cell in a geometry path.
The instance is an index into tally distribcell filter arrays.
Parameters
----------
@ -55,24 +55,31 @@ class Geometry(object):
lattice passed through. For the case of the lattice, a tuple should
be provided to indicate which coordinates in the lattice should be
entered. This should be in the form: (lat_id, i_x, i_y, i_z)
filter_offset : int
An integer that specifies which offset map the filter is using
Returns
-------
offset : int
Location in the results array for the path and filter
instance : int
Index in tally results array for distribcell filters
"""
# Find the distribcell index of the cell.
cells = self.get_all_cells()
if path[-1] in cells:
distribcell_index = cells[path[-1]].distribcell_index
else:
raise RuntimeError('Could not find cell {} specified in a \
distribcell filter'.format(path[-1]))
# Return memoize'd offset if possible
if (path, filter_offset) in self._offsets:
offset = self._offsets[(path, filter_offset)]
if (path, distribcell_index) in self._offsets:
offset = self._offsets[(path, distribcell_index)]
# Begin recursive call to compute offset starting with the base Universe
else:
offset = self._root_universe.get_offset(path, filter_offset)
self._offsets[(path, filter_offset)] = offset
offset = self._root_universe.get_cell_instance(path,
distribcell_index)
self._offsets[(path, distribcell_index)] = offset
# Return the final offset
return offset

View file

@ -360,9 +360,6 @@ class StatePoint(object):
# Read the Filter type
filter_type = self._f['{0}{1}/type'.format(subbase, j)].value.decode()
# Read the Filter offset
offset = self._f['{0}{1}/offset'.format(subbase, j)].value
n_bins = self._f['{0}{1}/n_bins'.format(subbase, j)].value
# Read the bin values
@ -370,7 +367,6 @@ class StatePoint(object):
# Create Filter object
filter = openmc.Filter(filter_type, bins)
filter.offset = offset
filter.num_bins = n_bins
if filter_type == 'mesh':

View file

@ -272,6 +272,11 @@ class Summary(object):
cell.region = Region.from_expression(
region, {s.id: s for s in self.surfaces.values()})
# Get the distribcell index
ind = self._f['geometry/cells'][key]['distribcell_index'].value
if ind != 0:
cell.distribcell_index = ind
# Add the Cell to the global dictionary of all Cells
self.cells[index] = cell

View file

@ -63,6 +63,8 @@ class Cell(object):
that is used to translate (shift) the universe.
offsets : ndarray
Array of offsets used for distributed cell searches
distribcell_index : int
Index of this cell in distribcell arrays
"""
@ -76,6 +78,7 @@ class Cell(object):
self._rotation = None
self._translation = None
self._offsets = None
self._distribcell_index = None
def __eq__(self, other):
if not isinstance(other, Cell):
@ -122,6 +125,8 @@ class Cell(object):
string += '{0: <16}{1}{2}\n'.format('\tTranslation', '=\t',
self._translation)
string += '{0: <16}{1}{2}\n'.format('\tOffset', '=\t', self._offsets)
string += '{0: <16}{1}{2}\n'.format('\tDistribcell index', '=\t',
self._distribcell_index)
return string
@ -164,6 +169,10 @@ class Cell(object):
def offsets(self):
return self._offsets
@property
def distribcell_index(self):
return self._distribcell_index
@id.setter
def id(self, cell_id):
if cell_id is None:
@ -231,6 +240,11 @@ class Cell(object):
cv.check_type('cell region', region, Region)
self._region = region
@distribcell_index.setter
def distribcell_index(self, ind):
cv.check_type('distribcell index', ind, Integral)
self._distribcell_index = ind
def add_surface(self, surface, halfspace):
"""Add a half-space to the list of half-spaces whose intersection defines the
cell.
@ -271,7 +285,7 @@ class Cell(object):
else:
self.region = Intersection(self.region, region)
def get_offset(self, path, filter_offset):
def get_cell_instance(self, path, distribcell_index):
# Get the current element and remove it from the list
cell_id = path[0]
path = path[1:]
@ -282,12 +296,12 @@ class Cell(object):
# If the Cell is filled by a Universe
elif self._type == 'fill':
offset = self._offsets[filter_offset-1]
offset += self._fill.get_offset(path, filter_offset)
offset = self.offsets[distribcell_index-1]
offset += self.fill.get_cell_instance(path, distribcell_index)
# If the Cell is filled by a Lattice
else:
offset = self._fill.get_offset(path, filter_offset)
offset = self.fill.get_cell_instance(path, distribcell_index)
return offset
@ -591,7 +605,7 @@ class Universe(object):
self._cells.clear()
def get_offset(self, path, filter_offset):
def get_cell_instance(self, path, distribcell_index):
# Get the current element and remove it from the list
path = path[1:]
@ -599,7 +613,7 @@ class Universe(object):
cell_id = path[0]
# Make a recursive call to the Cell within this Universe
offset = self._cells[cell_id].get_offset(path, filter_offset)
offset = self.cells[cell_id].get_cell_instance(path, distribcell_index)
# Return the offset computed at all nested Universe levels
return offset
@ -1059,21 +1073,22 @@ class RectLattice(Lattice):
cv.check_greater_than('lattice pitch', dim, 0.0)
self._pitch = pitch
def get_offset(self, path, filter_offset):
def get_cell_instance(self, path, distribcell_index):
# Get the current element and remove it from the list
i = path[0]
path = path[1:]
# For 2D Lattices
if len(self._dimension) == 2:
offset = self._offsets[i[1]-1, i[2]-1, 0, filter_offset-1]
offset += self._universes[i[1]][i[2]].get_offset(path, filter_offset)
offset = self._offsets[i[1]-1, i[2]-1, 0, distribcell_index-1]
offset += self._universes[i[1]][i[2]].get_cell_instance(path,
distribcell_index)
# For 3D Lattices
else:
offset = self._offsets[i[1]-1, i[2]-1, i[3]-1, filter_offset-1]
offset += self._universes[i[1]-1][i[2]-1][i[3]-1].get_offset(path,
filter_offset)
offset = self._offsets[i[1]-1, i[2]-1, i[3]-1, distribcell_index-1]
offset += self._universes[i[1]-1][i[2]-1][i[3]-1].get_cell_instance(
path, distribcell_index)
return offset

View file

@ -131,10 +131,13 @@ module geometry_header
integer, allocatable :: offset (:) ! Distribcell offset for tally
! counter
integer, allocatable :: region(:) ! Definition of spatial region as
! Boolean expression of half-spaces
! Boolean expression of half-spaces
integer, allocatable :: rpn(:) ! Reverse Polish notation for region
! expression
logical :: simple ! Is the region simple (intersections only)
! expression
logical :: simple ! Is the region simple (intersections
! only)
integer :: distribcell_index ! Index corresponding to this cell in
! distribcell arrays
! Rotation matrix and translation vector
real(8), allocatable :: translation(:)

View file

@ -1050,12 +1050,12 @@ contains
do i = 1, n_tallies
t => tallies(i)
do j = 1, t%n_filters
filter => t%filters(j)
do j = 1, t % n_filters
filter => t % filters(j)
if (filter%type == FILTER_DISTRIBCELL) then
if (.not. cell_list%contains(filter%int_bins(1))) then
call cell_list%add(filter%int_bins(1))
if (filter % type == FILTER_DISTRIBCELL) then
if (.not. cell_list % contains(filter % int_bins(1))) then
call cell_list % add(filter % int_bins(1))
end if
end if
@ -1066,8 +1066,8 @@ contains
! to determine the number of offset tables to allocate
do i = 1, n_universes
univ => universes(i)
do j = 1, univ%n_cells
if (cell_list%contains(univ%cells(j))) then
do j = 1, univ % n_cells
if (cell_list % contains(univ % cells(j))) then
n_maps = n_maps + 1
end if
end do
@ -1086,32 +1086,14 @@ contains
found(:,:) = .false.
k = 1
! Search through universes for distributed cells and assign each one a
! unique distribcell array index.
do i = 1, n_universes
univ => universes(i)
do j = 1, univ%n_cells
if (cell_list%contains(univ%cells(j))) then
! Loop over all tallies
do l = 1, n_tallies
t => tallies(l)
do m = 1, t%n_filters
filter => t%filters(m)
! Loop over only distribcell filters
! If filter points to cell we just found, set offset index
if (filter%type == FILTER_DISTRIBCELL) then
if (filter%int_bins(1) == univ%cells(j)) then
filter%offset = k
end if
end if
end do
end do
univ_list(k) = univ%id
do j = 1, univ % n_cells
if (cell_list % contains(univ % cells(j))) then
cells(univ % cells(j)) % distribcell_index = k
univ_list(k) = univ % id
k = k + 1
end if
end do
@ -1119,26 +1101,26 @@ contains
! Allocate the offset tables for lattices
do i = 1, n_lattices
lat => lattices(i)%obj
lat => lattices(i) % obj
select type(lat)
type is (RectLattice)
allocate(lat%offset(n_maps, lat%n_cells(1), lat%n_cells(2), &
lat%n_cells(3)))
allocate(lat % offset(n_maps, lat % n_cells(1), lat % n_cells(2), &
lat % n_cells(3)))
type is (HexLattice)
allocate(lat%offset(n_maps, 2 * lat%n_rings - 1, &
2 * lat%n_rings - 1, lat%n_axial))
allocate(lat % offset(n_maps, 2 * lat % n_rings - 1, &
2 * lat % n_rings - 1, lat % n_axial))
end select
lat%offset(:, :, :, :) = 0
lat % offset(:, :, :, :) = 0
end do
! Allocate offset table for fill cells
do i = 1, n_cells
if (cells(i)%material == NONE) then
allocate(cells(i)%offset(n_maps))
if (cells(i) % material == NONE) then
allocate(cells(i) % offset(n_maps))
end if
end do

View file

@ -1049,8 +1049,9 @@ contains
do i = 1, n_cells
c => cells(i)
! Initialize the number of cell instances - this is a base case for distribcells
! Initialize distribcell instances and distribcell index
c % instances = 0
c % distribcell_index = NONE
! Get pointer to i-th cell node
call get_list_item(node_cell_list, i, node_cell)

View file

@ -1378,8 +1378,7 @@ contains
label = ''
univ => universes(BASE_UNIVERSE)
offset = 0
call find_offset(t % filters(i_filter) % offset, &
t % filters(i_filter) % int_bins(1), &
call find_offset(t % filters(i_filter) % int_bins(1), &
univ, bin-1, offset, label)
case (FILTER_SURFACE)
i = t % filters(i_filter) % int_bins(bin)
@ -1413,15 +1412,15 @@ contains
! with the given offset
!===============================================================================
recursive subroutine find_offset(map, goal, univ, final, offset, path)
recursive subroutine find_offset(goal, univ, final, offset, path)
integer, intent(in) :: map ! Index in maps vector
integer, intent(in) :: goal ! The target cell ID
integer, intent(in) :: goal ! The target cell index
type(Universe), intent(in) :: univ ! Universe to begin search
integer, intent(in) :: final ! Target offset
integer, intent(inout) :: offset ! Current offset
character(*), intent(inout) :: path ! Path to offset
integer :: map ! Index in maps vector
integer :: i, j ! Index over cells
integer :: k, l, m ! Indices in lattice
integer :: old_k, old_l, old_m ! Previous indices in lattice
@ -1436,6 +1435,9 @@ contains
type(Universe), pointer :: next_univ ! Next universe to loop through
class(Lattice), pointer :: lat ! Pointer to current lattice
! Get the distribcell index for this cell
map = cells(goal) % distribcell_index
n = univ % n_cells
! Write to the geometry stack
@ -1447,17 +1449,13 @@ contains
! Look through all cells in this universe
do i = 1, n
cell_index = univ % cells(i)
c => cells(cell_index)
! If the cell ID matches the goal and the offset matches final,
! write to the geometry stack
if (cell_dict % get_key(c % id) == goal .AND. offset == final) then
path = trim(path) // "->" // to_str(c%id)
! If the cell matches the goal and the offset matches final, write to the
! geometry stack
if (univ % cells(i) == goal .AND. offset == final) then
c => cells(univ % cells(i))
path = trim(path) // "->" // to_str(c % id)
return
end if
end do
! Find the fill cell or lattice cell that we need to enter
@ -1537,7 +1535,7 @@ contains
offset = c % offset(map) + offset
next_univ => universes(c % fill)
call find_offset(map, goal, next_univ, final, offset, path)
call find_offset(goal, next_univ, final, offset, path)
return
! ====================================================================
@ -1577,7 +1575,7 @@ contains
path = trim(path) // "(" // trim(to_str(k)) // &
"," // trim(to_str(l)) // "," // &
trim(to_str(m)) // ")"
call find_offset(map, goal, next_univ, final, offset, path)
call find_offset(goal, next_univ, final, offset, path)
return
else
old_m = m
@ -1593,7 +1591,7 @@ contains
path = trim(path) // "(" // trim(to_str(old_k)) // &
"," // trim(to_str(old_l)) // "," // &
trim(to_str(old_m)) // ")"
call find_offset(map, goal, next_univ, final, offset, path)
call find_offset(goal, next_univ, final, offset, path)
return
end if
@ -1638,8 +1636,7 @@ contains
trim(to_str(k - lat % n_rings)) // "," // &
trim(to_str(l - lat % n_rings)) // "," // &
trim(to_str(m)) // ")"
call find_offset(map, goal, next_univ, final, offset, &
path)
call find_offset(goal, next_univ, final, offset, path)
return
else
old_m = m
@ -1656,7 +1653,7 @@ contains
trim(to_str(old_k - lat % n_rings)) // "," // &
trim(to_str(old_l - lat % n_rings)) // "," // &
trim(to_str(old_m)) // ")"
call find_offset(map, goal, next_univ, final, offset, path)
call find_offset(goal, next_univ, final, offset, path)
return
end if

View file

@ -268,7 +268,6 @@ contains
call write_dataset(filter_group, "type", "delayedgroup")
end select
call write_dataset(filter_group, "offset", tally%filters(j)%offset)
call write_dataset(filter_group, "n_bins", tally%filters(j)%n_bins)
if (tally % filters(j) % type == FILTER_ENERGYIN .or. &
tally % filters(j) % type == FILTER_ENERGYOUT .or. &

View file

@ -196,6 +196,8 @@ contains
end do
call write_dataset(cell_group, "region", adjustl(region_spec))
call write_dataset(cell_group, "distribcell_index", c % distribcell_index)
call close_group(cell_group)
end do CELL_LOOP
@ -548,7 +550,6 @@ contains
filter_group = create_group(tally_group, "filter " // trim(to_str(j)))
! Write number of bins for this filter
call write_dataset(filter_group, "offset", t%filters(j)%offset)
call write_dataset(filter_group, "n_bins", t%filters(j)%n_bins)
! Write filter bins

View file

@ -1712,6 +1712,7 @@ contains
integer :: j
integer :: n ! number of bins for single filter
integer :: offset ! offset for distribcell
integer :: distribcell_index ! index in distribcell arrays
real(8) :: E ! particle energy
real(8) :: theta, phi ! Polar and Azimuthal Angles, respectively
type(TallyObject), pointer :: t
@ -1756,12 +1757,14 @@ contains
case (FILTER_DISTRIBCELL)
! determine next distribcell bin
distribcell_index = cells(t % filters(i) % int_bins(1)) &
% distribcell_index
matching_bins(i) = NO_BIN_FOUND
offset = 0
do j = 1, p % n_coord
if (cells(p % coord(j) % cell) % type == CELL_FILL) then
offset = offset + cells(p % coord(j) % cell) % &
offset(t % filters(i) % offset)
offset(distribcell_index)
elseif(cells(p % coord(j) % cell) % type == CELL_LATTICE) then
if (lattices(p % coord(j + 1) % lattice) % obj &
% are_valid_indices([&
@ -1769,7 +1772,7 @@ contains
p % coord(j + 1) % lattice_y, &
p % coord(j + 1) % lattice_z])) then
offset = offset + lattices(p % coord(j + 1) % lattice) % obj % &
offset(t % filters(i) % offset, &
offset(distribcell_index, &
p % coord(j + 1) % lattice_x, &
p % coord(j + 1) % lattice_y, &
p % coord(j + 1) % lattice_z)

View file

@ -55,7 +55,6 @@ module tally_header
type TallyFilter
integer :: type = NONE
integer :: n_bins = 0
integer :: offset = 0 ! Only used for distribcell filters
integer, allocatable :: int_bins(:)
real(8), allocatable :: real_bins(:) ! Only used for energy filters
end type TallyFilter