crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-29 06:35:48 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge branch 'develop' into sparse-tallies

Browse source
This commit is contained in:
wbinventor@gmail.com 2016-01-04 19:23:39 -05:00
commit f62b057cfb
20 changed files with 186 additions and 169 deletions
Download patch file Download diff file Expand all files Collapse all files

4
docs/source/usersguide/output/statepoint.rst
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.

9
docs/source/usersguide/output/summary.rst
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.
@ -306,6 +310,11 @@ The current revision of the summary file format is 1.
than the number of user-specified scores since each score might have
multiple scoring bins, e.g., scatter-PN.
**/tallies/tally <uid>/moment_orders** (*char[][]*)
Tallying moment orders for Legendre and spherical harmonic tally expansions
(*e.g.*, 'P2', 'Y1,2', etc.).
**/tallies/tally <uid>/score_bins** (*char[][]*)
Scoring bins for the tally.

16
openmc/filter.py
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

31
openmc/geometry.py
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

26
openmc/opencg_compatible.py
View file

@ -725,11 +725,11 @@ def get_openmc_cell(opencg_cell):
else:
openmc_cell.fill = get_openmc_material(fill)
if opencg_cell.rotation:
if opencg_cell.rotation is not None:
rotation = np.asarray(opencg_cell.rotation, dtype=np.float64)
openmc_cell.rotation = rotation
if opencg_cell.translation:
if opencg_cell.translation is not None:
translation = np.asarray(opencg_cell.translation, dtype=np.float64)
openmc_cell.translation = translation
@ -881,16 +881,28 @@ def get_opencg_lattice(openmc_lattice):
universes = openmc_lattice.universes
outer = openmc_lattice.outer
# Convert 2D dimension to 3D for OpenCG
if len(dimension) == 2:
new_dimension = np.ones(3, dtype=np.int)
new_dimension[:2] = dimension
dimension = new_dimension
# Convert 2D pitch to 3D for OpenCG
if len(pitch) == 2:
new_pitch = np.ones(3, dtype=np.float64) * np.inf
new_pitch = np.ones(3, dtype=np.float64) * np.finfo(np.float64).max
new_pitch[:2] = pitch
pitch = new_pitch
# Convert 2D lower left to 3D for OpenCG
if len(lower_left) == 2:
new_lower_left = np.ones(3, dtype=np.float64)
new_lower_left = np.ones(3, dtype=np.float64) * np.finfo(np.float64).min
new_lower_left[:2] = lower_left
lower_left = new_lower_left
# Convert 2D universes array to 3D for OpenCG
if len(universes.shape) == 2:
universes.shape = (1,) + universes.shape
# Initialize an empty array for the OpenCG nested Universes in this Lattice
universe_array = np.ndarray(tuple(np.array(dimension)[::-1]),
dtype=opencg.Universe)
@ -905,7 +917,7 @@ def get_opencg_lattice(openmc_lattice):
for z in range(dimension[2]):
for y in range(dimension[1]):
for x in range(dimension[0]):
universe_id = universes[x][dimension[1]-y-1][z].id
universe_id = universes[z][y][x].id
universe_array[z][y][x] = unique_universes[universe_id]
opencg_lattice = opencg.Lattice(lattice_id, name)
@ -963,7 +975,7 @@ def get_openmc_lattice(opencg_lattice):
outer = opencg_lattice.outside
# Initialize an empty array for the OpenMC nested Universes in this Lattice
universe_array = np.ndarray(tuple(np.array(dimension)),
universe_array = np.ndarray(tuple(np.array(dimension)[::-1]),
dtype=openmc.Universe)
# Create OpenMC Universes for each unique nested Universe in this Lattice
@ -977,7 +989,7 @@ def get_openmc_lattice(opencg_lattice):
for y in range(dimension[1]):
for x in range(dimension[0]):
universe_id = universes[z][y][x].id
universe_array[x][y][z] = unique_universes[universe_id]
universe_array[z][y][x] = unique_universes[universe_id]
# Reverse y-dimension in array to match ordering in OpenCG
universe_array = universe_array[:, ::-1, :]

6
openmc/statepoint.py
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':
@ -404,7 +400,7 @@ class StatePoint(object):
for j in range(i+1, n_filters):
filter.stride *= tally.filters[j].num_bins
# Read scattering moment order strings (e.g., P3, Y-1,2, etc.)
# Read scattering moment order strings (e.g., P3, Y1,2, etc.)
moments = self._f['{0}{1}/moment_orders'.format(
base, tally_key)].value

35
openmc/summary.py
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
@ -330,11 +335,8 @@ class Summary(object):
self._f['geometry/lattices'][key]['lower_left'][...]
pitch = self._f['geometry/lattices'][key]['pitch'][...]
outer = self._f['geometry/lattices'][key]['outer'].value
universe_ids = \
self._f['geometry/lattices'][key]['universes'][...]
universe_ids = np.swapaxes(universe_ids, 0, 1)
universe_ids = np.swapaxes(universe_ids, 1, 2)
self._f['geometry/lattices'][key]['universes'][...]
# Create the Lattice
lattice = openmc.RectLattice(lattice_id=lattice_id, name=name)
@ -350,22 +352,22 @@ class Summary(object):
universes = \
np.ndarray(tuple(universe_ids.shape), dtype=openmc.Universe)
for x in range(universe_ids.shape[0]):
for z in range(universe_ids.shape[0]):
for y in range(universe_ids.shape[1]):
for z in range(universe_ids.shape[2]):
universes[x, y, z] = \
self.get_universe_by_id(universe_ids[x, y, z])
for x in range(universe_ids.shape[2]):
universes[z, y, x] = \
self.get_universe_by_id(universe_ids[z, y, x])
# Transpose, reverse y-dimension for appropriate ordering
shape = universes.shape
universes = np.transpose(universes, (1, 0, 2))
universes.shape = shape
universes = universes[:, ::-1, :]
# Use 2D NumPy array to store lattice universes for 2D lattices
if len(dimension) == 2:
universes = np.squeeze(universes)
universes = np.atleast_2d(universes)
# Set the universes for the lattice
lattice.universes = universes
if offsets is not None:
offsets = np.swapaxes(offsets, 0, 1)
offsets = np.swapaxes(offsets, 1, 2)
offsets = np.swapaxes(offsets, 0, 2)
lattice.offsets = offsets
# Add the Lattice to the global dictionary of all Lattices
@ -521,7 +523,7 @@ class Summary(object):
# Create Tally object and assign basic properties
tally = openmc.Tally(tally_id, tally_name)
# Read scattering moment order strings (e.g., P3, Y-1,2, etc.)
# Read scattering moment order strings (e.g., P3, Y1,2, etc.)
moments = self._f['{0}/moment_orders'.format(subbase)].value
# Read score metadata
@ -532,7 +534,6 @@ class Summary(object):
# If this is a moment, use generic moment order
pattern = r'-n$|-pn$|-yn$'
score = re.sub(pattern, '-' + moments[j].decode(), score)
tally.add_score(score)
# Read filter metadata

40
openmc/tallies.py
View file

@ -67,7 +67,9 @@ class Tally(object):
triggers : list of openmc.trigger.Trigger
List of tally triggers
num_scores : Integral
Total number of user-specified scores
Total number of scores, accounting for the fact that a single
user-specified score, e.g. scatter-P3 or flux-Y2,2, might have multiple
bins
num_filter_bins : Integral
Total number of filter bins accounting for all filters
num_bins : Integral
@ -1938,7 +1940,7 @@ class Tally(object):
data = self.get_values(
filters=filters, filter_bins=filter_bins, value='mean')
indices = self.get_filter_indices(filters, filter_bins)
self._mean[indices, :, :] = data
self.mean[indices, :, :] = data
# Adjust the std_dev data array to relect the new filter order
if self.std_dev is not None:
@ -1947,7 +1949,7 @@ class Tally(object):
data = self.get_values(
filters=filters, filter_bins=filter_bins, value='std_dev')
indices = self.get_filter_indices(filters, filter_bins)
self._std_dev[indices, :, :] = data
self.std_dev[indices, :, :] = data
def _swap_nuclides(self, nuclide1, nuclide2):
"""Reverse the ordering of two nuclides in this tally
@ -2004,17 +2006,17 @@ class Tally(object):
# Adjust the mean data array to relect the new nuclide order
if self.mean is not None:
nuclide1_mean = self._mean[:, nuclide1_index, :].copy()
nuclide2_mean = self._mean[:, nuclide2_index, :].copy()
self._mean[:, nuclide2_index, :] = nuclide1_mean
self._mean[:, nuclide1_index, :] = nuclide2_mean
nuclide1_mean = self.mean[:, nuclide1_index, :].copy()
nuclide2_mean = self.mean[:, nuclide2_index, :].copy()
self.mean[:, nuclide2_index, :] = nuclide1_mean
self.mean[:, nuclide1_index, :] = nuclide2_mean
# Adjust the std_dev data array to relect the new nuclide order
if self.std_dev is not None:
nuclide1_std_dev = self._std_dev[:, nuclide1_index, :].copy()
nuclide2_std_dev = self._std_dev[:, nuclide2_index, :].copy()
self._std_dev[:, nuclide2_index, :] = nuclide1_std_dev
self._std_dev[:, nuclide1_index, :] = nuclide2_std_dev
nuclide1_std_dev = self.std_dev[:, nuclide1_index, :].copy()
nuclide2_std_dev = self.std_dev[:, nuclide2_index, :].copy()
self.std_dev[:, nuclide2_index, :] = nuclide1_std_dev
self.std_dev[:, nuclide1_index, :] = nuclide2_std_dev
def _swap_scores(self, score1, score2):
"""Reverse the ordering of two scores in this tally
@ -2076,17 +2078,17 @@ class Tally(object):
# Adjust the mean data array to relect the new nuclide order
if self.mean is not None:
score1_mean = self._mean[:, :, score1_index].copy()
score2_mean = self._mean[:, :, score2_index].copy()
self._mean[:, :, score2_index] = score1_mean
self._mean[:, :, score1_index] = score2_mean
score1_mean = self.mean[:, :, score1_index].copy()
score2_mean = self.mean[:, :, score2_index].copy()
self.mean[:, :, score2_index] = score1_mean
self.mean[:, :, score1_index] = score2_mean
# Adjust the std_dev data array to relect the new nuclide order
if self.std_dev is not None:
score1_std_dev = self._std_dev[:, :, score1_index].copy()
score2_std_dev = self._std_dev[:, :, score2_index].copy()
self._std_dev[:, :, score2_index] = score1_std_dev
self._std_dev[:, :, score1_index] = score2_std_dev
score1_std_dev = self.std_dev[:, :, score1_index].copy()
score2_std_dev = self.std_dev[:, :, score2_index].copy()
self.std_dev[:, :, score2_index] = score1_std_dev
self.std_dev[:, :, score1_index] = score2_std_dev
def __add__(self, other):
"""Adds this tally to another tally or scalar value.

45
openmc/universe.py
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
@ -807,7 +821,7 @@ class Lattice(object):
def universes(self, universes):
cv.check_iterable_type('lattice universes', universes, Universe,
min_depth=2, max_depth=3)
self._universes = universes
self._universes = np.asarray(universes)
def get_unique_universes(self):
"""Determine all unique universes in the lattice
@ -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
@ -1119,7 +1134,7 @@ class RectLattice(Lattice):
for z in range(self._dimension[2]):
for y in range(self._dimension[1]):
for x in range(self._dimension[0]):
universe = self._universes[x][y][z]
universe = self._universes[z][y][x]
# Append Universe ID to the Lattice XML subelement
universe_ids += '{0} '.format(universe._id)
@ -1137,7 +1152,7 @@ class RectLattice(Lattice):
else:
for y in range(self._dimension[1]):
for x in range(self._dimension[0]):
universe = self._universes[x][y]
universe = self._universes[y][x]
# Append Universe ID to Lattice XML subelement
universe_ids += '{0} '.format(universe._id)

4
src/ace_header.F90
View file

@ -48,7 +48,7 @@ module ace_header
integer :: MT ! ENDF MT value
real(8) :: Q_value ! Reaction Q value
integer :: multiplicity ! Number of secondary particles released
type(Tab1), allocatable :: multiplicity_E ! Energy-dependent neutron yield
type(Tab1), pointer :: multiplicity_E => null() ! Energy-dependent neutron yield
integer :: threshold ! Energy grid index of threshold
logical :: scatter_in_cm ! scattering system in center-of-mass?
logical :: multiplicity_with_E = .false. ! Flag to indicate E-dependent multiplicity
@ -308,6 +308,8 @@ module ace_header
class(Reaction), intent(inout) :: this ! The Reaction object to clear
if (associated(this % multiplicity_E)) deallocate(this % multiplicity_E)
if (associated(this % edist)) then
call this % edist % clear()
deallocate(this % edist)

2
src/constants.F90
View file

@ -14,7 +14,7 @@ module constants
integer, parameter :: REVISION_STATEPOINT = 14
integer, parameter :: REVISION_PARTICLE_RESTART = 1
integer, parameter :: REVISION_TRACK = 1
integer, parameter :: REVISION_SUMMARY = 1
integer, parameter :: REVISION_SUMMARY = 2
! ============================================================================
! ADJUSTABLE PARAMETERS

9
src/geometry_header.F90
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(:)

2
src/hdf5_interface.F90
View file

@ -1898,7 +1898,7 @@ contains
logical :: mpio
integer :: hdf5_err
integer :: driver
integer(HID_T) :: driver
integer(HID_T) :: file_id
integer(HID_T) :: fapl_id

60
src/initialize.F90
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

3
src/input_xml.F90
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)

37
src/output.F90
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

1
src/state_point.F90
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. &

17
src/summary.F90
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
@ -355,16 +357,22 @@ contains
call write_dataset(lattice_group, "type", "rectangular")
! Write lattice dimensions, lower left corner, and pitch
call write_dataset(lattice_group, "dimension", lat%n_cells)
call write_dataset(lattice_group, "lower_left", lat%lower_left)
if (lat % is_3d) then
call write_dataset(lattice_group, "dimension", lat % n_cells)
call write_dataset(lattice_group, "lower_left", lat % lower_left)
else
call write_dataset(lattice_group, "dimension", lat % n_cells(1:2))
call write_dataset(lattice_group, "lower_left", lat % lower_left)
end if
! Write lattice universes.
allocate(lattice_universes(lat%n_cells(1), lat%n_cells(2), &
&lat%n_cells(3)))
do j = 1, lat%n_cells(1)
do k = 1, lat%n_cells(2)
do k = 0, lat%n_cells(2) - 1
do m = 1, lat%n_cells(3)
lattice_universes(j,k,m) = universes(lat%universes(j,k,m))%id
lattice_universes(j, k+1, m) = &
universes(lat%universes(j, lat%n_cells(2) - k, m))%id
end do
end do
end do
@ -542,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

7
src/tally.F90
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)

1
src/tally_header.F90
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