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

Merge remote-tracking branch 'upstream/develop' into new-scatt-mat

Browse source
This commit is contained in:
Will Boyd 2017-03-07 13:53:33 -05:00
commit b3687b48ee
31 changed files with 1047 additions and 1196 deletions
Download patch file Download diff file Expand all files Collapse all files

6
CMakeLists.txt
View file

@ -225,7 +225,11 @@ endif()
#===============================================================================
add_library(pugixml src/pugixml/pugixml_c.cpp src/pugixml/pugixml.cpp)
set_property(TARGET pugixml PROPERTY CXX_STANDARD 11)
if(CMAKE_VERSION VERSION_LESS 3.1)
target_compile_options(pugixml PRIVATE -std=c++11)
else()
set_property(TARGET pugixml PROPERTY CXX_STANDARD 11)
endif()
add_library(pugixml_fortran src/pugixml/pugixml_f.F90)
target_link_libraries(pugixml_fortran pugixml)

13
docs/source/io_formats/summary.rst
View file

@ -38,9 +38,6 @@ The current version of the summary file format is 5.0.
is an array if the cell uses distributed materials, otherwise it is
a scalar.
- **temperature** (*double[]*) -- Temperature of the cell in Kelvin.
- **offset** (*int[]*) -- Offsets used for distribcell tally
filter. This dataset is present only if fill_type is set to
'universe'.
- **translation** (*double[3]*) -- Translation applied to the fill
universe. This dataset is present only if fill_type is set to
'universe'.
@ -50,15 +47,6 @@ The current version of the summary file format is 5.0.
- **lattice** (*int*) -- Unique ID of the lattice which fills the
cell. Only present if fill_type is set to 'lattice'.
- **region** (*char[]*) -- Region specification for the cell.
- **distribcell_index** (*int*) -- Index of this cell in distribcell
arrays. Only present if this cell is listed in a distribcell filter
or if it uses distributed materials.
- **paths** (*char[][]*) -- The paths traversed through the CSG tree
to reach each distribcell instance. This consists of the integer
IDs for each universe, cell and lattice delimited by '->'. Each
lattice cell is specified by its (x,y) or (x,y,z) indices. Only
present if this cell is listed in a distribcell filter or if it
uses distributed materials.
**/geometry/surfaces/surface <uid>/**
@ -87,7 +75,6 @@ The current version of the summary file format is 5.0.
- **pitch** (*double[]*) -- Pitch of the lattice in centimeters.
- **outer** (*int*) -- Outer universe assigned to lattice cells
outside the defined range.
- **offsets** (*int[]*) -- Offsets used for distribcell tally filter.
- **universes** (*int[][][]*) -- Three-dimensional array of universes
assigned to each cell of the lattice.
- **dimension** (*int[]*) -- The number of lattice cells in each

910
docs/source/pythonapi/examples/pandas-dataframes.ipynb
View file

File diff suppressed because one or more lines are too long

75
openmc/cell.py
View file

@ -79,13 +79,12 @@ class Cell(object):
translation : Iterable of float
If the cell is filled with a universe, this array specifies a vector
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
distribcell_paths : list of str
The paths traversed through the CSG tree to reach each distribcell
instance
paths : list of str
The paths traversed through the CSG tree to reach each cell
instance. This property is initialized by calling the
:meth:`Geometry.determine_paths` method.
num_instances : int
The number of instances of this cell throughout the geometry.
volume : float
Volume of the cell in cm^3. This can either be set manually or
calculated in a stochastic volume calculation and added via the
@ -103,9 +102,7 @@ class Cell(object):
self._rotation_matrix = None
self._temperature = None
self._translation = None
self._offsets = None
self._distribcell_index = None
self._distribcell_paths = None
self._paths = []
self._volume = None
self._atoms = None
@ -162,8 +159,6 @@ class Cell(object):
string += '\t{0: <15}=\t{1}\n'.format('Temperature',
self.temperature)
string += '{: <16}=\t{}\n'.format('\tTranslation', self.translation)
string += '{: <16}=\t{}\n'.format('\tOffset', self.offsets)
string += '{: <16}=\t{}\n'.format('\tDistribcell index', self.distribcell_index)
return string
@ -212,22 +207,21 @@ class Cell(object):
def translation(self):
return self._translation
@property
def offsets(self):
return self._offsets
@property
def distribcell_index(self):
return self._distribcell_index
@property
def distribcell_paths(self):
return self._distribcell_paths
@property
def volume(self):
return self._volume
@property
def paths(self):
if not self._paths:
raise ValueError('Cell instance paths have not been determined. '
'Call the Geometry.determine_paths() method.')
return self._paths
@property
def num_instances(self):
return len(self.paths)
@id.setter
def id(self, cell_id):
if cell_id is None:
@ -316,11 +310,6 @@ class Cell(object):
else:
self._temperature = temperature
@offsets.setter
def offsets(self, offsets):
cv.check_type('cell offsets', offsets, Iterable)
self._offsets = offsets
@region.setter
def region(self, region):
if region is not None:
@ -333,17 +322,6 @@ class Cell(object):
cv.check_type('cell volume', volume, Real)
self._volume = volume
@distribcell_index.setter
def distribcell_index(self, ind):
cv.check_type('distribcell index', ind, Integral)
self._distribcell_index = ind
@distribcell_paths.setter
def distribcell_paths(self, distribcell_paths):
cv.check_iterable_type('distribcell_paths', distribcell_paths,
string_types)
self._distribcell_paths = distribcell_paths
def add_surface(self, surface, halfspace):
"""Add a half-space to the list of half-spaces whose intersection defines the
cell.
@ -406,23 +384,6 @@ class Cell(object):
else:
raise ValueError('No volume information found for this cell.')
def get_cell_instance(self, path, distribcell_index):
# If the Cell is filled by a Material
if self.fill_type in ('material', 'distribmat', 'void'):
offset = 0
# If the Cell is filled by a Universe
elif self.fill_type == 'universe':
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_cell_instance(path, distribcell_index)
return offset
def get_nuclides(self):
"""Returns all nuclides in the cell

154
openmc/filter.py
View file

@ -395,7 +395,7 @@ class Filter(object):
Keyword arguments
-----------------
distribcell_paths : bool
paths : bool
Only used for DistirbcellFilter. If True (default), expand
distribcell indices into multi-index columns describing the path
to that distribcell through the CSG tree. NOTE: This option assumes
@ -1041,6 +1041,42 @@ class EnergyoutFilter(EnergyFilter):
"""
def _path_to_levels(path):
"""Convert distribcell path to list of levels
Parameters
----------
path : str
Distribcell path
Returns
-------
list
List of levels in path
"""
# Split path into universes/cells/lattices
path_items = path.split('->')
# Pair together universe and cell information from the same level
idx = [i for i, item in enumerate(path_items) if item.startswith('u')]
for i in reversed(idx):
univ_id = int(path_items.pop(i)[1:])
cell_id = int(path_items.pop(i)[1:])
path_items.insert(i, ('universe', univ_id, cell_id))
# Reformat lattice into tuple
idx = [i for i, item in enumerate(path_items) if isinstance(item, str)]
for i in idx:
item = path_items.pop(i)[1:-1]
lat_id, lat_xyz = item.split('(')
lat_id = int(lat_id)
lat_xyz = tuple(int(x) for x in lat_xyz.split(','))
path_items.insert(i, ('lattice', lat_id, lat_xyz))
return path_items
class DistribcellFilter(Filter):
"""Bins tally event locations on instances of repeated cells.
@ -1059,14 +1095,14 @@ class DistribcellFilter(Filter):
stride : Integral
The number of filter, nuclide and score bins within each of this
filter's bins.
distribcell_paths : list of str
paths : list of str
The paths traversed through the CSG tree to reach each distribcell
instance (for 'distribcell' filters only)
"""
def __init__(self, bins):
self._distribcell_paths = None
self._paths = None
super(DistribcellFilter, self).__init__(bins)
@classmethod
@ -1079,20 +1115,16 @@ class DistribcellFilter(Filter):
out = cls(group['bins'].value)
out.num_bins = group['n_bins'].value
if 'paths' in group:
out.distribcell_paths = [str(path.decode()) for path in
group['paths'].value]
return out
@property
def distribcell_paths(self):
return self._distribcell_paths
def paths(self):
return self._paths
@distribcell_paths.setter
def distribcell_paths(self, distribcell_paths):
cv.check_iterable_type('distribcell_paths', distribcell_paths, str)
self._distribcell_paths = distribcell_paths
@paths.setter
def paths(self, paths):
cv.check_iterable_type('paths', paths, str)
self._paths = paths
def check_bins(self, bins):
if not len(bins) == 1:
@ -1127,7 +1159,7 @@ class DistribcellFilter(Filter):
Keyword arguments
-----------------
distribcell_paths : bool
paths : bool
If True (default), expand distribcell indices into multi-index
columns describing the path to that distribcell through the CSG
tree. NOTE: This option assumes that all distribcell paths are of
@ -1165,47 +1197,35 @@ class DistribcellFilter(Filter):
level_df = None
distribcell_paths = kwargs.setdefault('distribcell_paths', True)
paths = kwargs.setdefault('paths', True)
# Create Pandas Multi-index columns for each level in CSG tree
if distribcell_paths:
if paths:
# Distribcell paths require linked metadata from the Summary
if self.distribcell_paths is None:
if self.paths is None:
msg = 'Unable to construct distribcell paths since ' \
'the Summary is not linked to the StatePoint'
raise ValueError(msg)
# Make copy of array of distribcell paths to use in
# Pandas Multi-index column construction
distribcell_paths = copy.deepcopy(self.distribcell_paths)
num_offsets = len(distribcell_paths)
num_offsets = len(self.paths)
paths = [_path_to_levels(p) for p in self.paths]
# Loop over CSG levels in the distribcell paths
level_counter = 0
levels_remain = True
while levels_remain:
num_levels = len(paths[0])
for i_level in range(num_levels):
# Use level key as first index in Pandas Multi-index column
level_counter += 1
level_key = 'level {}'.format(level_counter)
# Use the first distribcell path to determine if level
# is a universe/cell or lattice level
first_path = distribcell_paths[0]
next_index = first_path.index('-')
level = first_path[:next_index]
# Trim universe/lattice info from path
first_path = first_path[next_index+2:]
level_key = 'level {}'.format(i_level + 1)
# Create a dictionary for this level for Pandas Multi-index
level_dict = OrderedDict()
# This level is a lattice (e.g., ID(x,y,z))
if '(' in level:
level_type = 'lattice'
# Use the first distribcell path to determine if level
# is a universe/cell or lattice level
path = paths[0]
if path[i_level][0] == 'lattice':
# Initialize prefix Multi-index keys
lat_id_key = (level_key, 'lat', 'id')
lat_x_key = (level_key, 'lat', 'x')
@ -1217,12 +1237,10 @@ class DistribcellFilter(Filter):
level_dict[lat_id_key] = np.empty(num_offsets)
level_dict[lat_x_key] = np.empty(num_offsets)
level_dict[lat_y_key] = np.empty(num_offsets)
level_dict[lat_z_key] = np.empty(num_offsets)
if len(path[i_level][2]) == 3:
level_dict[lat_z_key] = np.empty(num_offsets)
# This level is a universe / cell (e.g., ID->ID)
else:
level_type = 'universe'
# Initialize prefix Multi-index keys
univ_key = (level_key, 'univ', 'id')
cell_key = (level_key, 'cell', 'id')
@ -1232,52 +1250,22 @@ class DistribcellFilter(Filter):
level_dict[univ_key] = np.empty(num_offsets)
level_dict[cell_key] = np.empty(num_offsets)
# Determine any levels remain in path
if '-' not in first_path:
levels_remain = False
# Populate Multi-index arrays with all distribcell paths
for i, path in enumerate(distribcell_paths):
if level_type == 'lattice':
# Extract lattice ID, indices from path
next_index = path.index('-')
lat_id_indices = path[:next_index]
# Trim lattice info from distribcell path
distribcell_paths[i] = path[next_index+2:]
# Extract the lattice cell indices from the path
i1 = lat_id_indices.index('(')
i2 = lat_id_indices.index(')')
i3 = lat_id_indices[i1+1:i2]
for i, path in enumerate(paths):
level = path[i_level]
if level[0] == 'lattice':
# Assign entry to Lattice Multi-index column
level_dict[lat_id_key][i] = path[:i1]
level_dict[lat_x_key][i] = int(i3.split(',')[0]) - 1
level_dict[lat_y_key][i] = int(i3.split(',')[1]) - 1
level_dict[lat_z_key][i] = int(i3.split(',')[2]) - 1
level_dict[lat_id_key][i] = level[1]
level_dict[lat_x_key][i] = level[2][0]
level_dict[lat_y_key][i] = level[2][1]
if len(level[2]) == 3:
level_dict[lat_z_key][i] = level[2][2]
else:
# Extract universe ID from path
next_index = path.index('-')
universe_id = int(path[:next_index])
# Trim universe info from distribcell path
path = path[next_index+2:]
# Extract cell ID from path
if '-' in path:
next_index = path.index('-')
cell_id = int(path[:next_index])
distribcell_paths[i] = path[next_index+2:]
else:
cell_id = int(path)
distribcell_paths[i] = ''
# Assign entry to Universe, Cell Multi-index columns
level_dict[univ_key][i] = universe_id
level_dict[cell_key][i] = cell_id
level_dict[univ_key][i] = level[1]
level_dict[cell_key][i] = level[2]
# Tile the Multi-index columns
for level_key, level_bins in level_dict.items():
@ -1569,7 +1557,7 @@ class AzimuthalFilter(RealFilter):
'increasing'.format(bins, type(self))
raise ValueError(msg)
def get_pandas_dataframe(self, data_size, distribcell_paths=True):
def get_pandas_dataframe(self, data_size, paths=True):
"""Builds a Pandas DataFrame for the Filter's bins.
This method constructs a Pandas DataFrame object for the filter with

94
openmc/geometry.py
View file

@ -1,6 +1,8 @@
from collections import OrderedDict
from collections import OrderedDict, Iterable
from xml.etree import ElementTree as ET
from six import string_types
import openmc
from openmc.clean_xml import sort_xml_elements, clean_xml_indentation
from openmc.checkvalue import check_type
@ -42,12 +44,6 @@ class Geometry(object):
@root_universe.setter
def root_universe(self, root_universe):
check_type('root universe', root_universe, openmc.Universe)
if root_universe.id != 0:
msg = 'Unable to add root Universe "{0}" to Geometry since ' \
'it has ID="{1}" instead of ' \
'ID=0'.format(root_universe, root_universe.id)
raise ValueError(msg)
self._root_universe = root_universe
def add_volume_information(self, volume_calc):
@ -110,53 +106,52 @@ class Geometry(object):
"""
return self.root_universe.find(point)
def get_cell_instance(self, path):
"""Return the instance number for the final cell in a geometry path.
def get_instances(self, paths):
"""Return the instance number(s) for a cell/material in a geometry path.
The instance is an index into tally distribcell filter arrays.
The instance numbers are used as indices into distributed
material/temperature arrays and tally distribcell filter arrays.
Parameters
----------
path : list
A list of IDs that form the path to the target. It should begin with
0 for the base universe, and should cover every universe, cell, and
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)
paths : str or iterable of str
The path traversed through the CSG tree to reach a cell or material
instance. For example, 'u0->c10->l20(2,2,1)->u5->c5' would indicate
the cell instance whose first level is universe 0 and cell 10,
second level is lattice 20 position (2,2,1), and third level is
universe 5 and cell 5.
Returns
-------
instance : int
Index in tally results array for distribcell filters
int or list of int
Instance number(s) for the given path(s)
"""
# Make sure we are working with an iterable
return_list = (isinstance(paths, Iterable) and
not isinstance(paths, string_types))
path_list = paths if return_list else [paths]
# Extract the cell id from the path
last_index = path.rfind('>')
cell_id = int(path[last_index+1:])
indices = []
for p in path_list:
# Extract the cell id from the path
last_index = p.rfind('>')
last_path = p[last_index+1:]
uid = int(last_path[1:])
# Find the distribcell index of the cell.
cells = self.get_all_cells()
for cell in cells:
if cell.id == cell_id:
distribcell_index = cell.distribcell_index
break
else:
raise RuntimeError('Could not find cell {} specified in a \
distribcell filter'.format(cell_id))
# Get corresponding cell/material
if last_path[0] == 'c':
obj = self.get_all_cells()[uid]
elif last_path[0] == 'm':
obj = self.get_all_materials()[uid]
# Return memoize'd offset if possible
if (path, distribcell_index) in self._offsets:
offset = self._offsets[(path, distribcell_index)]
# Determine index in paths array
try:
indices.append(obj.paths.index(p))
except ValueError:
indices.append(None)
# Begin recursive call to compute offset starting with the base Universe
else:
offset = self._root_universe.get_cell_instance(path,
distribcell_index)
self._offsets[(path, distribcell_index)] = offset
# Return the final offset
return offset
return indices if return_list else indices[0]
def get_all_cells(self):
"""Return all cells in the geometry.
@ -455,3 +450,20 @@ class Geometry(object):
lattices = list(lattices)
lattices.sort(key=lambda x: x.id)
return lattices
def determine_paths(self):
"""Determine paths through CSG tree for cells and materials.
This method recursively traverses the CSG tree to determine each unique
path that reaches every cell and material. The paths are stored in the
:attr:`Cell.paths` and :attr:`Material.paths` attributes.
"""
# (Re-)initialize all cell instances to 0
for cell in self.get_all_cells().values():
cell._paths = []
for material in self.get_all_materials().values():
material._paths = []
# Recursively traverse the CSG tree to count all cell instances
self.root_universe._determine_paths()

231
openmc/lattice.py
View file

@ -132,11 +132,6 @@ class Lattice(object):
name = group['name'].value.decode() if 'name' in group else ''
lattice_type = group['type'].value.decode()
if 'offsets' in group:
offsets = group['offsets'][...]
else:
offsets = None
if lattice_type == 'rectangular':
dimension = group['dimension'][...]
lower_left = group['lower_left'][...]
@ -169,10 +164,6 @@ class Lattice(object):
# Set the universes for the lattice
lattice.universes = uarray
# Set the distribcell offsets for the lattice
if offsets is not None:
lattice.offsets = offsets
elif lattice_type == 'hexagonal':
n_rings = group['n_rings'].value
n_axial = group['n_axial'].value
@ -256,9 +247,6 @@ class Lattice(object):
# Lattice is 2D; extract the only axial level
lattice.universes = uarray[0]
if offsets is not None:
lattice.offsets = offsets
return lattice
def get_unique_universes(self):
@ -378,6 +366,54 @@ class Lattice(object):
return all_universes
def get_universe(self, idx):
"""Return universe corresponding to a lattice element index
Parameters
----------
idx : Iterable of int
Lattice element indices. For a rectangular lattice, the indices are
given in the :math:`(x,y)` or :math:`(x,y,z)` coordinate system. For
hexagonal lattices, they are given in the :math:`x,\alpha` or
:math:`x,\alpha,z` coordinate systems.
Returns
-------
openmc.Universe
Universe with given indices
"""
idx_u = self.get_universe_index(idx)
if self.ndim == 2:
return self.universes[idx_u[0]][idx_u[1]]
else:
return self.universes[idx_u[0]][idx_u[1]][idx_u[2]]
def find(self, point):
"""Find cells/universes/lattices which contain a given point
Parameters
----------
point : 3-tuple of float
Cartesian coordinates of the point
Returns
-------
list
Sequence of universes, cells, and lattices which are traversed to
find the given point
"""
idx, p = self.find_element(point)
if self.is_valid_index(idx):
u = self.get_universe(idx)
else:
if self.outer is not None:
u = self.outer
else:
return []
return [(self, idx)] + u.find(p)
class RectLattice(Lattice):
"""A lattice consisting of rectangular prisms.
@ -443,7 +479,6 @@ class RectLattice(Lattice):
# Initialize Lattice class attributes
self._lower_left = None
self._offsets = None
def __eq__(self, other):
if not isinstance(other, RectLattice):
@ -492,19 +527,6 @@ class RectLattice(Lattice):
string = string.rstrip('\n')
if self._offsets is not None:
string += '{0: <16}\n'.format('\tOffsets')
# Lattice cell offsets
for i, offset in enumerate(np.ravel(self._offsets)):
string += '{0} '.format(offset)
# Add a newline character when we reach end of row of cells
if (i+1) % self.shape[0] == 0:
string += '\n'
string = string.rstrip('\n')
return string
@property
@ -516,6 +538,22 @@ class RectLattice(Lattice):
return list(np.broadcast(*np.ogrid[
:self.shape[2], :self.shape[1], :self.shape[0]]))
@property
def _natural_indices(self):
"""Iterate over all possible (x,y) or (x,y,z) lattice element indices.
This property is used when constructing distributed cell and material
paths. Most importantly, the iteration order matches that used on the
Fortran side.
"""
if self.ndim == 2:
nx, ny = self.shape
return np.broadcast(*np.ogrid[:nx, :ny])
else:
nx, ny, nz = self.shape
return np.broadcast(*np.ogrid[:nx, :ny, :nz])
@property
def lower_left(self):
return self._lower_left
@ -524,10 +562,6 @@ class RectLattice(Lattice):
def ndim(self):
return len(self.pitch)
@property
def offsets(self):
return self._offsets
@property
def shape(self):
return self._universes.shape[::-1]
@ -538,11 +572,6 @@ class RectLattice(Lattice):
cv.check_length('lattice lower left corner', lower_left, 2, 3)
self._lower_left = lower_left
@offsets.setter
def offsets(self, offsets):
cv.check_type('lattice offsets', offsets, Iterable)
self._offsets = offsets
@Lattice.pitch.setter
def pitch(self, pitch):
cv.check_type('lattice pitch', pitch, Iterable, Real)
@ -557,34 +586,6 @@ class RectLattice(Lattice):
min_depth=2, max_depth=3)
self._universes = np.asarray(universes)
def get_cell_instance(self, path, distribcell_index):
# Extract the lattice element from the path
next_index = path.index('-')
lat_id_indices = path[:next_index]
path = path[next_index+2:]
# Extract the lattice cell indices from the path
i1 = lat_id_indices.index('(')
i2 = lat_id_indices.index(')')
i = lat_id_indices[i1+1:i2]
lat_x = int(i.split(',')[0]) - 1
lat_y = int(i.split(',')[1]) - 1
lat_z = int(i.split(',')[2]) - 1
# For 2D Lattices
if self.ndim == 2:
offset = self._offsets[lat_z, lat_y, lat_x, distribcell_index-1]
offset += self._universes[lat_x][lat_y].get_cell_instance(
path, distribcell_index)
# For 3D Lattices
else:
offset = self._offsets[lat_z, lat_y, lat_x, distribcell_index-1]
offset += self._universes[lat_z][lat_y][lat_x].get_cell_instance(
path, distribcell_index)
return offset
def find_element(self, point):
"""Determine index of lattice element and local coordinates for a point
@ -681,32 +682,6 @@ class RectLattice(Lattice):
0 <= idx[1] < self.shape[1] and
0 <= idx[2] < self.shape[2])
def find(self, point):
"""Find cells/universes/lattices which contain a given point
Parameters
----------
point : 3-tuple of float
Cartesian coordinatesof the point
Returns
-------
list
Sequence of universes, cells, and lattices which are traversed to
find the given point
"""
idx, p = self.find_element(point)
if self.is_valid_index(idx):
idx_u = self.get_universe_index(idx)
u = self.universes[idx_u]
else:
if self.outer is not None:
u = self.outer
else:
return []
return [(self, idx)] + u.find(p)
def create_xml_subelement(self, xml_element):
# Determine if XML element already contains subelement for this Lattice
@ -922,6 +897,35 @@ class HexLattice(Lattice):
for r in range(self.num_rings)
for i in range(max(6*(self.num_rings - 1 - r), 1))]
@property
def _natural_indices(self):
"""Iterate over all possible (x,alpha) or (x,alpha,z) lattice element
indices.
This property is used when constructing distributed cell and material
paths. Most importantly, the iteration order matches that used on the
Fortran side.
"""
r = self.num_rings
if self.num_axial is None:
for a in range(-r + 1, r):
for x in range(-r + 1, r):
idx = (x, a)
if self.is_valid_index(idx):
yield idx
else:
for z in range(self.num_axial):
for a in range(-r + 1, r):
for x in range(-r + 1, r):
idx = (x, a, z)
if self.is_valid_index(idx):
yield idx
@property
def ndim(self):
return 2 if isinstance(self.universes[0][0], openmc.Universe) else 3
@center.setter
def center(self, center):
cv.check_type('lattice center', center, Iterable, Real)
@ -948,28 +952,15 @@ class HexLattice(Lattice):
# The Universes within each sub-list are ordered from the "top" in a
# clockwise fashion.
# Check to see if the given universes look like a 2D or a 3D array.
if isinstance(self._universes[0][0], openmc.Universe):
n_dims = 2
elif isinstance(self._universes[0][0][0], openmc.Universe):
n_dims = 3
else:
msg = 'HexLattice ID={0:d} does not appear to be either 2D or ' \
'3D. Make sure set_universes was given a two-deep or ' \
'three-deep iterable of universes.'.format(self._id)
raise RuntimeError(msg)
# Set the number of axial positions.
if n_dims == 3:
if self.ndim == 3:
self._num_axial = len(self._universes)
else:
self._num_axial = None
# Set the number of rings and make sure this number is consistent for
# all axial positions.
if n_dims == 3:
if self.ndim == 3:
self._num_rings = len(self._universes[0])
for rings in self._universes:
if len(rings) != self._num_rings:
@ -981,7 +972,7 @@ class HexLattice(Lattice):
self._num_rings = len(self._universes)
# Make sure there are the correct number of elements in each ring.
if n_dims == 3:
if self.ndim == 3:
for axial_slice in self._universes:
# Check the center ring.
if len(axial_slice[-1]) != 1:
@ -1153,36 +1144,6 @@ class HexLattice(Lattice):
else:
return g < self.num_rings and 0 <= idx[2] < self.num_axial
def find(self, point):
"""Find cells/universes/lattices which contain a given point
Parameters
----------
point : 3-tuple of float
Cartesian coordinatesof the point
Returns
-------
list
Sequence of universes, cells, and lattices which are traversed to
find the given point
"""
idx, p = self.find_element(point)
if self.is_valid_index(idx):
idx_u = self.get_universe_index(idx)
if self.num_axial is None:
u = self.universes[idx_u[0]][idx_u[1]]
else:
u = self.universes[idx_u[0]][idx_u[1]][idx_u[2]]
else:
if self.outer is not None:
u = self.outer
else:
return []
return [(self, idx)] + u.find(p)
def create_xml_subelement(self, xml_element):
# Determine if XML element already contains subelement for this Lattice
path = './hex_lattice[@id=\'{0}\']'.format(self._id)

23
openmc/material.py
View file

@ -75,6 +75,12 @@ class Material(object):
Volume of the material in cm^3. This can either be set manually or
calculated in a stochastic volume calculation and added via the
:meth:`Material.add_volume_information` method.
paths : list of str
The paths traversed through the CSG tree to reach each material
instance. This property is initialized by calling the
:meth:`Geometry.determine_paths` method.
num_instances : int
The number of instances of this material throughout the geometry.
"""
@ -86,6 +92,7 @@ class Material(object):
self._density = None
self._density_units = ''
self._depletable = False
self._paths = []
self._volume = None
self._atoms = {}
@ -195,6 +202,17 @@ class Material(object):
def depletable(self):
return self._depletable
@property
def paths(self):
if not self._paths:
raise ValueError('Material instance paths have not been determined. '
'Call the Geometry.determine_paths() method.')
return self._paths
@property
def num_instances(self):
return len(self.paths)
@property
def elements(self):
return self._elements
@ -275,6 +293,11 @@ class Material(object):
cv.check_type('material volume', volume, Real)
self._volume = volume
@num_instances.setter
def num_instances(self, num_instances):
cv.check_type('num_instances', num_instances, Integral)
self._num_instances = num_instances
@classmethod
def from_hdf5(cls, group):
"""Create material from HDF5 group

13
openmc/mgxs/mdgxs.py
View file

@ -774,7 +774,7 @@ class MDGXS(MGXS):
modified.write('\n\\end{document}')
def get_pandas_dataframe(self, groups='all', nuclides='all',
xs_type='macro', distribcell_paths=True,
xs_type='macro', paths=True,
delayed_groups='all'):
"""Build a Pandas DataFrame for the MDGXS data.
@ -795,7 +795,7 @@ class MDGXS(MGXS):
xs_type: {'macro', 'micro'}
Return macro or micro cross section in units of cm^-1 or barns.
Defaults to 'macro'.
distribcell_paths : bool, optional
paths : bool, optional
Construct columns for distribcell tally filters (default is True).
The geometric information in the Summary object is embedded into
a Multi-index column with a geometric "path" to each distribcell
@ -830,8 +830,7 @@ class MDGXS(MGXS):
# Use tally summation to sum across all nuclides
xs_tally = self.xs_tally.summation(nuclides=self.get_nuclides())
df = xs_tally.get_pandas_dataframe(
distribcell_paths=distribcell_paths)
df = xs_tally.get_pandas_dataframe(paths=paths)
# Remove nuclide column since it is homogeneous and redundant
if self.domain_type == 'mesh':
@ -842,13 +841,11 @@ class MDGXS(MGXS):
# If the user requested a specific set of nuclides
elif self.by_nuclide and nuclides != 'all':
xs_tally = self.xs_tally.get_slice(nuclides=nuclides)
df = xs_tally.get_pandas_dataframe(
distribcell_paths=distribcell_paths)
df = xs_tally.get_pandas_dataframe(paths=paths)
# If the user requested all nuclides, keep nuclide column in dataframe
else:
df = self.xs_tally.get_pandas_dataframe(
distribcell_paths=distribcell_paths)
df = self.xs_tally.get_pandas_dataframe(paths=paths)
# Remove the score column since it is homogeneous and redundant
if self.domain_type == 'mesh':

25
openmc/mgxs/mgxs.py
View file

@ -1845,7 +1845,7 @@ class MGXS(object):
modified.write('\n\\end{document}')
def get_pandas_dataframe(self, groups='all', nuclides='all',
xs_type='macro', distribcell_paths=True):
xs_type='macro', paths=True):
"""Build a Pandas DataFrame for the MGXS data.
This method leverages :meth:`openmc.Tally.get_pandas_dataframe`, but
@ -1865,7 +1865,7 @@ class MGXS(object):
xs_type: {'macro', 'micro'}
Return macro or micro cross section in units of cm^-1 or barns.
Defaults to 'macro'.
distribcell_paths : bool, optional
paths : bool, optional
Construct columns for distribcell tally filters (default is True).
The geometric information in the Summary object is embedded into
a Multi-index column with a geometric "path" to each distribcell
@ -1895,8 +1895,7 @@ class MGXS(object):
# Use tally summation to sum across all nuclides
xs_tally = self.xs_tally.summation(nuclides=self.get_nuclides())
df = xs_tally.get_pandas_dataframe(
distribcell_paths=distribcell_paths)
df = xs_tally.get_pandas_dataframe(paths=paths)
# Remove nuclide column since it is homogeneous and redundant
if self.domain_type == 'mesh':
@ -1907,13 +1906,11 @@ class MGXS(object):
# If the user requested a specific set of nuclides
elif self.by_nuclide and nuclides != 'all':
xs_tally = self.xs_tally.get_slice(nuclides=nuclides)
df = xs_tally.get_pandas_dataframe(
distribcell_paths=distribcell_paths)
df = xs_tally.get_pandas_dataframe(paths=paths)
# If the user requested all nuclides, keep nuclide column in dataframe
else:
df = self.xs_tally.get_pandas_dataframe(
distribcell_paths=distribcell_paths)
df = self.xs_tally.get_pandas_dataframe(paths=paths)
# Remove the score column since it is homogeneous and redundant
if self.domain_type == 'mesh':
@ -4388,7 +4385,7 @@ class ScatterMatrixXS(MatrixMGXS):
return xs
def get_pandas_dataframe(self, groups='all', nuclides='all', moment='all',
xs_type='macro', distribcell_paths=True):
xs_type='macro', paths=True):
"""Build a Pandas DataFrame for the MGXS data.
This method leverages :meth:`openmc.Tally.get_pandas_dataframe`, but
@ -4412,7 +4409,7 @@ class ScatterMatrixXS(MatrixMGXS):
xs_type: {'macro', 'micro'}
Return macro or micro cross section in units of cm^-1 or barns.
Defaults to 'macro'.
distribcell_paths : bool, optional
paths : bool, optional
Construct columns for distribcell tally filters (default is True).
The geometric information in the Summary object is embedded into a
Multi-index column with a geometric "path" to each distribcell
@ -4432,7 +4429,7 @@ class ScatterMatrixXS(MatrixMGXS):
"""
df = super(ScatterMatrixXS, self).get_pandas_dataframe(
groups, nuclides, xs_type, distribcell_paths)
groups, nuclides, xs_type, paths)
if self.scatter_format == 'legendre':
# Add a moment column to dataframe
@ -5614,7 +5611,7 @@ class Chi(MGXS):
return xs
def get_pandas_dataframe(self, groups='all', nuclides='all',
xs_type='macro', distribcell_paths=False):
xs_type='macro', paths=False):
"""Build a Pandas DataFrame for the MGXS data.
This method leverages :meth:`openmc.Tally.get_pandas_dataframe`, but
@ -5634,7 +5631,7 @@ class Chi(MGXS):
xs_type: {'macro', 'micro'}
Return macro or micro cross section in units of cm^-1 or barns.
Defaults to 'macro'.
distribcell_paths : bool, optional
paths : bool, optional
Construct columns for distribcell tally filters (default is True).
The geometric information in the Summary object is embedded into
a Multi-index column with a geometric "path" to each distribcell
@ -5655,7 +5652,7 @@ class Chi(MGXS):
# Build the dataframe using the parent class method
df = super(Chi, self).get_pandas_dataframe(
groups, nuclides, xs_type, distribcell_paths=distribcell_paths)
groups, nuclides, xs_type, paths=paths)
# If user requested micro cross sections, multiply by the atom
# densities to cancel out division made by the parent class method

4
openmc/statepoint.py
View file

@ -648,6 +648,8 @@ class StatePoint(object):
if isinstance(tally_filter, (openmc.DistribcellFilter)):
cell_id = tally_filter.bins[0]
cell = cells[cell_id]
tally_filter.distribcell_paths = cell.distribcell_paths
if not cell._paths:
summary.geometry.determine_paths()
tally_filter.paths = cell.paths
self._summary = summary

27
openmc/summary.py
View file

@ -121,10 +121,6 @@ class Summary(object):
cell = openmc.Cell(cell_id=cell_id, name=name)
if fill_type == 'universe':
if 'offset' in group:
offset = group['offset'][...]
cell.offsets = offset
if 'translation' in group:
translation = group['translation'][...]
translation = np.asarray(translation, dtype=np.float64)
@ -145,14 +141,6 @@ class Summary(object):
if region:
cell.region = Region.from_expression(region, surfaces)
# Get the distribcell data
if 'distribcell_index' in group:
ind = group['distribcell_index'].value
cell.distribcell_index = ind
paths = group['paths'][...]
paths = [str(path.decode()) for path in paths]
cell.distribcell_paths = paths
# Add the Cell to the global dictionary of all Cells
cells[cell.id] = cell
@ -175,6 +163,11 @@ class Summary(object):
def _finalize_geometry(self, cells, cell_fills, universes, lattices):
materials = {m.id: m for m in self.materials}
# Keep track of universes that are used as fills. That way, we can
# determine which universe is NOT used as a fill (and hence is the root
# universe)
fill_univ_ids = set()
# Iterate over all Cells and add fill Materials, Universes and Lattices
for cell_id, (fill_type, fill_id) in cell_fills.items():
# Retrieve the object corresponding to the fill type and ID
@ -186,14 +179,20 @@ class Summary(object):
fill = materials[fill_id] if fill_id > 0 else None
elif fill_type == 'universe':
fill = universes[fill_id]
fill_univ_ids.add(fill_id)
else:
fill = lattices[fill_id]
for idx in fill._natural_indices:
univ = fill.get_universe(idx)
fill_univ_ids.add(univ.id)
# Set the fill for the Cell
cells[cell_id].fill = fill
# Set the root universe for the Geometry
self.geometry.root_universe = universes[0]
# Determine root universe for geometry
non_fill = set(universes.keys()) - fill_univ_ids
self.geometry.root_universe = universes[non_fill.pop()]
def add_volume_information(self, volume_calc):
"""Add volume information to the geometry within the summary file

82
openmc/surface.py
View file

@ -1882,7 +1882,8 @@ class Halfspace(Region):
else str(self.surface.id)
def get_rectangular_prism(width, height, axis='z', origin=(0., 0.)):
def get_rectangular_prism(width, height, axis='z', origin=(0., 0.),
boundary_type='transmission'):
"""Get an infinite rectangular prism from four planar surfaces.
Parameters
@ -1900,6 +1901,9 @@ def get_rectangular_prism(width, height, axis='z', origin=(0., 0.)):
Origin of the prism. The two floats correspond to (y,z), (x,z) or
(x,y) for prisms parallel to the x, y or z axis, respectively.
Defaults to (0., 0.).
boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
Boundary condition that defines the behavior for particles hitting the
surfaces comprising the rectangular prism (default is 'transmission').
Returns
-------
@ -1914,28 +1918,41 @@ def get_rectangular_prism(width, height, axis='z', origin=(0., 0.)):
check_type('origin', origin, Iterable, Real)
if axis == 'x':
min_y = YPlane(name='minimum y', y0=-width/2.+origin[0])
max_y = YPlane(name='maximum y', y0=+width/2.+origin[0])
min_z = ZPlane(name='minimum z', z0=-height/2.+origin[1])
max_z = ZPlane(name='maximum z', z0=+height/2.+origin[1])
min_y = YPlane(name='minimum y', y0=-width/2.+origin[0],
boundary_type=boundary_type)
max_y = YPlane(name='maximum y', y0=+width/2.+origin[0],
boundary_type=boundary_type)
min_z = ZPlane(name='minimum z', z0=-height/2.+origin[1],
boundary_type=boundary_type)
max_z = ZPlane(name='maximum z', z0=+height/2.+origin[1],
boundary_type=boundary_type)
prism = +min_y & -max_y & +min_z & -max_z
elif axis == 'y':
min_x = XPlane(name='minimum x', x0=-width/2.+origin[0])
max_x = XPlane(name='maximum x', x0=+width/2.+origin[0])
min_z = ZPlane(name='minimum z', z0=-height/2.+origin[1])
max_z = ZPlane(name='maximum z', z0=+height/2.+origin[1])
min_x = XPlane(name='minimum x', x0=-width/2.+origin[0],
boundary_type=boundary_type)
max_x = XPlane(name='maximum x', x0=+width/2.+origin[0],
boundary_type=boundary_type)
min_z = ZPlane(name='minimum z', z0=-height/2.+origin[1],
boundary_type=boundary_type)
max_z = ZPlane(name='maximum z', z0=+height/2.+origin[1],
boundary_type=boundary_type)
prism = +min_x & -max_x & +min_z & -max_z
else:
min_x = XPlane(name='minimum x', x0=-width/2.+origin[0])
max_x = XPlane(name='maximum x', x0=+width/2.+origin[0])
min_y = YPlane(name='minimum y', y0=-height/2.+origin[1])
max_y = YPlane(name='maximum y', y0=+height/2.+origin[1])
min_x = XPlane(name='minimum x', x0=-width/2.+origin[0],
boundary_type=boundary_type)
max_x = XPlane(name='maximum x', x0=+width/2.+origin[0],
boundary_type=boundary_type)
min_y = YPlane(name='minimum y', y0=-height/2.+origin[1],
boundary_type=boundary_type)
max_y = YPlane(name='maximum y', y0=+height/2.+origin[1],
boundary_type=boundary_type)
prism = +min_x & -max_x & +min_y & -max_y
return prism
def get_hexagonal_prism(edge_length=1., orientation='y'):
def get_hexagonal_prism(edge_length=1., orientation='y',
boundary_type='transmission'):
"""Create a hexagon region from six surface planes.
Parameters
@ -1946,6 +1963,9 @@ def get_hexagonal_prism(edge_length=1., orientation='y'):
An 'x' orientation means that two sides of the hexagon are parallel to
the x-axis and a 'y' orientation means that two sides of the hexagon are
parallel to the y-axis.
boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
Boundary condition that defines the behavior for particles hitting the
surfaces comprising the hexagonal prism (default is 'transmission').
Returns
-------
@ -1960,10 +1980,18 @@ def get_hexagonal_prism(edge_length=1., orientation='y'):
right = XPlane(x0=sqrt(3.)/2.*l)
left = XPlane(x0=-sqrt(3.)/2.*l)
c = sqrt(3.)/3.
upper_right = Plane(A=c, B=1., D=l) # y = -x/sqrt(3) + a
upper_left = Plane(A=-c, B=1., D=l) # y = x/sqrt(3) + a
lower_right = Plane(A=-c, B=1., D=-l) # y = x/sqrt(3) - a
lower_left = Plane(A=c, B=1., D=-l) # y = -x/sqrt(3) - a
# y = -x/sqrt(3) + a
upper_right = Plane(A=c, B=1., D=l, boundary_type=boundary_type)
# y = x/sqrt(3) + a
upper_left = Plane(A=-c, B=1., D=l, boundary_type=boundary_type)
# y = x/sqrt(3) - a
lower_right = Plane(A=-c, B=1., D=-l, boundary_type=boundary_type)
# y = -x/sqrt(3) - a
lower_left = Plane(A=c, B=1., D=-l, boundary_type=boundary_type)
return Intersection(-right, +left, -upper_right, -upper_left,
+lower_right, +lower_left)
@ -1971,9 +1999,17 @@ def get_hexagonal_prism(edge_length=1., orientation='y'):
top = YPlane(y0=sqrt(3.)/2.*l)
bottom = YPlane(y0=-sqrt(3.)/2.*l)
c = sqrt(3.)
upper_right = Plane(A=c, B=1., D=c*l) # y = -sqrt(3)*(x - a)
lower_right = Plane(A=-c, B=1., D=-c*l) # y = sqrt(3)*(x + a)
lower_left = Plane(A=c, B=1., D=-c*l) # y = -sqrt(3)*(x + a)
upper_left = Plane(A=-c, B=1., D=c*l) # y = sqrt(3)*(x + a)
# y = -sqrt(3)*(x - a)
upper_right = Plane(A=c, B=1., D=c*l, boundary_type=boundary_type)
# y = sqrt(3)*(x + a)
lower_right = Plane(A=-c, B=1., D=-c*l, boundary_type=boundary_type)
# y = -sqrt(3)*(x + a)
lower_left = Plane(A=c, B=1., D=-c*l, boundary_type=boundary_type)
# y = sqrt(3)*(x + a)
upper_left = Plane(A=-c, B=1., D=c*l, boundary_type=boundary_type)
return Intersection(-top, +bottom, -upper_right, +lower_right,
+lower_left, -upper_left)
+lower_left, -upper_left)

7
openmc/tallies.py
View file

@ -1517,8 +1517,7 @@ class Tally(object):
return data
def get_pandas_dataframe(self, filters=True, nuclides=True, scores=True,
derivative=True, distribcell_paths=True,
float_format='{:.2e}'):
derivative=True, paths=True, float_format='{:.2e}'):
"""Build a Pandas DataFrame for the Tally data.
This method constructs a Pandas DataFrame object for the Tally data
@ -1538,7 +1537,7 @@ class Tally(object):
Include columns with score bin information (default is True).
derivative : bool
Include columns with differential tally info (default is True).
distribcell_paths : bool, optional
paths : bool, optional
Construct columns for distribcell tally filters (default is True).
The geometric information in the Summary object is embedded into a
Multi-index column with a geometric "path" to each distribcell
@ -1581,7 +1580,7 @@ class Tally(object):
# Append each Filter's DataFrame to the overall DataFrame
for self_filter in self.filters:
filter_df = self_filter.get_pandas_dataframe(
data_size, distribcell_paths=distribcell_paths)
data_size, paths=paths)
df = pd.concat([df, filter_df], axis=1)
# Include DataFrame column for nuclides if user requested it

63
openmc/universe.py
View file

@ -60,10 +60,6 @@ class Universe(object):
# Values - Cells
self._cells = OrderedDict()
# Keys - Cell IDs
# Values - Offsets
self._cell_offsets = OrderedDict()
if cells is not None:
self.add_cells(cells)
@ -392,27 +388,6 @@ class Universe(object):
self._cells.clear()
def get_cell_instance(self, path, distribcell_index):
# Pop off the root Universe ID from the path
next_index = path.index('-')
path = path[next_index+2:]
# Extract the Cell ID from the path
if '-' in path:
next_index = path.index('-')
cell_id = int(path[:next_index])
path = path[next_index+2:]
else:
cell_id = int(path)
path = ''
# Make a recursive call to the Cell within this Universe
offset = self.cells[cell_id].get_cell_instance(path, distribcell_index)
# Return the offset computed at all nested Universe levels
return offset
def get_nuclides(self):
"""Returns all nuclides in the universe
@ -519,3 +494,41 @@ class Universe(object):
# Append the Universe ID to the subelement and add to Element
cell_element.set("universe", str(self._id))
xml_element.append(cell_element)
def _determine_paths(self, path=''):
"""Count the number of instances for each cell in the universe, and
record the count in the :attr:`Cell.num_instances` properties."""
univ_path = path + 'u{}'.format(self.id)
for cell in self.cells.values():
cell_path = '{}->c{}'.format(univ_path, cell.id)
# If universe-filled, recursively count cells in filling universe
if cell.fill_type == 'universe':
cell.fill._determine_paths(cell_path + '->')
# If lattice-filled, recursively call for all universes in lattice
elif cell.fill_type == 'lattice':
latt = cell.fill
# Count instances in each universe in the lattice
for index in latt._natural_indices:
latt_path = '{}->l{}({})->'.format(
cell_path, latt.id, ",".join(str(x) for x in index))
univ = latt.get_universe(index)
univ._determine_paths(latt_path)
else:
if cell.fill_type == 'material':
mat = cell.fill
elif cell.fill_type == 'distribmat':
mat = cell.fill[len(cell._paths)]
else:
mat = None
if mat is not None:
mat._paths.append('{}->m{}'.format(cell_path, mat.id))
# Append current path
cell._paths.append(cell_path)

143
src/geometry.F90
View file

@ -156,7 +156,7 @@ contains
do j = 1, n_coord
p % n_coord = j
univ => universes(p % coord(j) % universe)
n = univ % n_cells
n = size(univ % cells)
! loop through each cell on this level
do i = 1, n
@ -216,7 +216,7 @@ contains
else
use_search_cells = .false.
univ => universes(p % coord(j) % universe)
n = univ % n_cells
n = size(univ % cells)
end if
CELL_LOOP: do i = 1, n
@ -1087,9 +1087,9 @@ contains
! routine is called once upon initialization.
!===============================================================================
subroutine calc_offsets(goal, map, univ, counts, found)
subroutine calc_offsets(univ_id, map, univ, counts, found)
integer, intent(in) :: goal ! target universe ID
integer, intent(in) :: univ_id ! target universe ID
integer, intent(in) :: map ! map index in vector of maps
type(Universe), intent(in) :: univ ! universe searching in
integer, intent(inout) :: counts(:,:) ! target count
@ -1104,7 +1104,7 @@ contains
type(Universe), pointer :: next_univ ! next universe to cycle through
class(Lattice), pointer :: lat ! pointer to current lattice
n = univ % n_cells
n = size(univ % cells)
offset = 0
do i = 1, n
@ -1126,7 +1126,7 @@ contains
! Count contents of this cell
next_univ => universes(c % fill)
offset = offset + count_target(next_univ, counts, found, goal, map)
offset = offset + count_target(next_univ, counts, found, univ_id, map)
! Move into the next universe
next_univ => universes(c % fill)
@ -1150,7 +1150,7 @@ contains
lat % offset(map, j, k, m) = offset
next_univ => universes(lat % universes(j, k, m))
offset = offset + &
count_target(next_univ, counts, found, goal, map)
count_target(next_univ, counts, found, univ_id, map)
end do
end do
end do
@ -1171,7 +1171,7 @@ contains
lat % offset(map, j, k, m) = offset
next_univ => universes(lat % universes(j, k, m))
offset = offset + &
count_target(next_univ, counts, found, goal, map)
count_target(next_univ, counts, found, univ_id, map)
end if
end do
end do
@ -1188,12 +1188,12 @@ contains
! universe ID beginning with the universe given.
!===============================================================================
recursive function count_target(univ, counts, found, goal, map) result(count)
recursive function count_target(univ, counts, found, univ_id, map) result(count)
type(Universe), intent(in) :: univ ! universe to search through
integer, intent(inout) :: counts(:,:) ! target count
logical, intent(inout) :: found(:,:) ! target found
integer, intent(in) :: goal ! target universe ID
integer, intent(in) :: univ_id ! target universe ID
integer, intent(in) :: map ! current map
integer :: i ! index over cells
@ -1213,7 +1213,7 @@ contains
! If this is the target, it can't contain itself.
! Count = 1, then quit
if (univ % id == goal) then
if (univ % id == univ_id) then
count = 1
counts(universe_dict % get_key(univ % id), map) = 1
found(universe_dict % get_key(univ % id), map) = .true.
@ -1221,7 +1221,7 @@ contains
end if
count = 0
n = univ % n_cells
n = size(univ % cells)
do i = 1, n
@ -1241,12 +1241,12 @@ contains
next_univ => universes(c % fill)
! Found target - stop since target cannot contain itself
if (next_univ % id == goal) then
if (next_univ % id == univ_id) then
count = count + 1
return
end if
count = count + count_target(next_univ, counts, found, goal, map)
count = count + count_target(next_univ, counts, found, univ_id, map)
c => cells(cell_index)
! ====================================================================
@ -1267,13 +1267,13 @@ contains
next_univ => universes(lat % universes(j, k, m))
! Found target - stop since target cannot contain itself
if (next_univ % id == goal) then
if (next_univ % id == univ_id) then
count = count + 1
cycle
end if
count = count + &
count_target(next_univ, counts, found, goal, map)
count_target(next_univ, counts, found, univ_id, map)
end do
end do
@ -1295,13 +1295,13 @@ contains
next_univ => universes(lat % universes(j, k, m))
! Found target - stop since target cannot contain itself
if (next_univ % id == goal) then
if (next_univ % id == univ_id) then
count = count + 1
cycle
end if
count = count + &
count_target(next_univ, counts, found, goal, map)
count_target(next_univ, counts, found, univ_id, map)
end if
end do
end do
@ -1326,81 +1326,64 @@ contains
type(Universe), intent(in) :: univ ! universe to search through
integer :: i ! index over cells
integer :: j, k, m ! indices in lattice
integer :: n ! number of cells to search
integer :: cell_index ! index in cells array
type(Cell), pointer :: c ! pointer to current cell
type(Universe), pointer :: next_univ ! next universe to loop through
class(Lattice), pointer :: lat ! pointer to current lattice
integer :: i ! index over cells
integer :: j, k, m ! indices in lattice
integer :: n ! number of cells to search
n = univ % n_cells
n = size(univ % cells)
do i = 1, n
associate (c => cells(univ % cells(i)))
c % instances = c % instances + 1
cell_index = univ % cells(i)
! ====================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
if (c % type == FILL_UNIVERSE) then
! get pointer to cell
c => cells(cell_index)
c % instances = c % instances + 1
call count_instance(universes(c % fill))
! ====================================================================
! AT LOWEST UNIVERSE, TERMINATE SEARCH
if (c % type == FILL_MATERIAL) then
! ====================================================================
! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
elseif (c % type == FILL_LATTICE) then
! ====================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
elseif (c % type == FILL_UNIVERSE) then
! Set current lattice
associate (lat => lattices(c % fill) % obj)
next_univ => universes(c % fill)
select type (lat)
type is (RectLattice)
call count_instance(next_univ)
c => cells(cell_index)
! ====================================================================
! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
elseif (c % type == FILL_LATTICE) then
! Set current lattice
lat => lattices(c % fill) % obj
select type (lat)
type is (RectLattice)
! Loop over lattice coordinates
do j = 1, lat % n_cells(1)
do k = 1, lat % n_cells(2)
do m = 1, lat % n_cells(3)
next_univ => universes(lat % universes(j, k, m))
call count_instance(next_univ)
! Loop over lattice coordinates
do j = 1, lat % n_cells(1)
do k = 1, lat % n_cells(2)
do m = 1, lat % n_cells(3)
call count_instance(universes(lat % universes(j, k, m)))
end do
end do
end do
end do
end do
type is (HexLattice)
type is (HexLattice)
! Loop over lattice coordinates
do m = 1, lat % n_axial
do k = 1, 2*lat % n_rings - 1
do j = 1, 2*lat % n_rings - 1
! This array location is never used
if (j + k < lat % n_rings + 1) then
cycle
! This array location is never used
else if (j + k > 3*lat % n_rings - 1) then
cycle
else
next_univ => universes(lat % universes(j, k, m))
call count_instance(next_univ)
end if
! Loop over lattice coordinates
do m = 1, lat % n_axial
do k = 1, 2*lat % n_rings - 1
do j = 1, 2*lat % n_rings - 1
! This array location is never used
if (j + k < lat % n_rings + 1) then
cycle
! This array location is never used
else if (j + k > 3*lat % n_rings - 1) then
cycle
else
call count_instance(universes(lat % universes(j, k, m)))
end if
end do
end do
end do
end do
end do
end select
end if
end select
end associate
end if
end associate
end do
end subroutine count_instance
@ -1423,7 +1406,7 @@ contains
class(Lattice), pointer :: lat ! pointer to current lattice
levels_below = 0
do i = 1, univ % n_cells
do i = 1, size(univ % cells)
c => cells(univ % cells(i))
! ====================================================================

5
src/geometry_header.F90
View file

@ -17,7 +17,6 @@ module geometry_header
type Universe
integer :: id ! Unique ID
integer :: type ! Type
integer :: n_cells ! # of cells within
integer, allocatable :: cells(:) ! List of cells within
real(8) :: x0 ! Translation in x-coordinate
real(8) :: y0 ! Translation in y-coordinate
@ -153,8 +152,8 @@ module geometry_header
real(8), allocatable :: rotation_matrix(:,:)
end type Cell
! array index of universe 0
integer :: BASE_UNIVERSE
! array index of the root universe
integer :: root_universe = -1
contains

88
src/initialize.F90
View file

@ -16,11 +16,11 @@ module initialize
use geometry, only: neighbor_lists, count_instance, calc_offsets, &
maximum_levels
use geometry_header, only: Cell, Universe, Lattice, RectLattice, HexLattice,&
&BASE_UNIVERSE
root_universe
use global
use hdf5_interface, only: file_open, read_attribute, file_close, &
hdf5_bank_t, hdf5_integer8_t
use input_xml, only: read_input_xml, cells_in_univ_dict, read_plots_xml
use input_xml, only: read_input_xml, read_plots_xml
use material_header, only: Material
use message_passing
use mgxs_data, only: read_mgxs, create_macro_xs
@ -81,9 +81,6 @@ contains
! XML files because we need the PRNG to be initialized first
if (run_mode == MODE_PLOTTING) call read_plots_xml()
! Set up universe structures
call prepare_universes()
! Use dictionaries to redefine index pointers
call adjust_indices()
@ -97,7 +94,7 @@ contains
! Check to make sure there are not too many nested coordinate levels in the
! geometry since the coordinate list is statically allocated for performance
! reasons
if (maximum_levels(universes(BASE_UNIVERSE)) > MAX_COORD) then
if (maximum_levels(universes(root_universe)) > MAX_COORD) then
call fatal_error("Too many nested coordinate levels in the geometry. &
&Try increasing the maximum number of coordinate levels by &
&providing the CMake -Dmaxcoord= option.")
@ -416,79 +413,6 @@ contains
end subroutine read_command_line
!===============================================================================
! PREPARE_UNIVERSES allocates the universes array and determines the cells array
! for each universe.
!===============================================================================
subroutine prepare_universes()
integer :: i ! index in cells array
integer :: i_univ ! index in universes array
integer :: n_cells_in_univ ! number of cells in a universe
integer, allocatable :: index_cell_in_univ(:) ! the index in the univ%cells
! array for each universe
type(ElemKeyValueII), pointer :: pair_list => null()
type(ElemKeyValueII), pointer :: current => null()
type(ElemKeyValueII), pointer :: next => null()
type(Universe), pointer :: univ => null()
type(Cell), pointer :: c => null()
allocate(universes(n_universes))
! We also need to allocate the cell count lists for each universe. The logic
! for this is a little more convoluted. In universe_dict, the (key,value)
! pairs are the id of the universe and the index in the array. In
! cells_in_univ_dict, it's the id of the universe and the number of cells.
pair_list => universe_dict%keys()
current => pair_list
do while (associated(current))
! Find index of universe in universes array
i_univ = current%value
univ => universes(i_univ)
univ%id = current%key
! Check for lowest level universe
if (univ%id == 0) BASE_UNIVERSE = i_univ
! Find cell count for this universe
n_cells_in_univ = cells_in_univ_dict%get_key(univ%id)
! Allocate cell list for universe
allocate(univ%cells(n_cells_in_univ))
univ%n_cells = n_cells_in_univ
! Move to next universe
next => current%next
deallocate(current)
current => next
end do
! Also allocate a list for keeping track of where cells have been assigned
! in each universe
allocate(index_cell_in_univ(n_universes))
index_cell_in_univ = 0
do i = 1, n_cells
c => cells(i)
! Get pointer to corresponding universe
i_univ = universe_dict%get_key(c%universe)
univ => universes(i_univ)
! Increment the index for the cells array within the Universe object and
! then store the index of the Cell object in that array
index_cell_in_univ(i_univ) = index_cell_in_univ(i_univ) + 1
univ%cells(index_cell_in_univ(i_univ)) = i
end do
! Clear dictionary
call cells_in_univ_dict%clear()
end subroutine prepare_universes
!===============================================================================
! ADJUST_INDICES changes the values for 'surfaces' for each cell and the
! material index assigned to each to the indices in the surfaces and material
@ -794,7 +718,7 @@ contains
if (.not. distribcell_active) return
! Count the number of instances of each cell.
call count_instance(universes(BASE_UNIVERSE))
call count_instance(universes(root_universe))
! Set the number of bins in all distribcell filters.
do i = 1, n_tallies
@ -881,7 +805,7 @@ contains
! Compute the number of unique universes containing these distribcells
! to determine the number of offset tables to allocate
do i = 1, n_universes
do j = 1, universes(i) % n_cells
do j = 1, size(universes(i) % cells)
if (cell_list % contains(universes(i) % cells(j))) then
n_maps = n_maps + 1
end if
@ -904,7 +828,7 @@ contains
! unique distribcell array index.
k = 1
do i = 1, n_universes
do j = 1, universes(i) % n_cells
do j = 1, size(universes(i) % cells)
if (cell_list % contains(universes(i) % cells(j))) then
cells(universes(i) % cells(j)) % distribcell_index = k
univ_list(k) = universes(i) % id

88
src/input_xml.F90
View file

@ -10,7 +10,7 @@ module input_xml
use energy_grid, only: grid_method, n_log_bins
use error, only: fatal_error, warning
use geometry_header, only: Cell, Lattice, RectLattice, HexLattice, &
get_temperatures
get_temperatures, root_universe
use global
use hdf5_interface
use list_header, only: ListChar, ListInt, ListReal
@ -35,9 +35,6 @@ module input_xml
implicit none
save
type(DictIntInt) :: cells_in_univ_dict ! Used to count how many cells each
! universe contains
contains
!===============================================================================
@ -1061,7 +1058,7 @@ contains
integer :: i, j, k, m, i_x, i_a, input_index
integer :: n, n_mats, n_x, n_y, n_z, n_rings, n_rlats, n_hlats
integer :: universe_num
integer :: univ_id
integer :: n_cells_in_univ
integer :: coeffs_reqd
integer :: i_xmin, i_xmax, i_ymin, i_ymax, i_zmin, i_zmax
@ -1089,6 +1086,11 @@ contains
type(XMLNode), allocatable :: node_hlat_list(:)
type(VectorInt) :: tokens
type(VectorInt) :: rpn
type(VectorInt) :: fill_univ_ids ! List of fill universe IDs
type(VectorInt) :: univ_ids ! List of all universe IDs
type(DictIntInt) :: cells_in_univ_dict ! Used to count how many cells each
! universe contains
! Display output message
call write_message("Reading geometry XML file...", 5)
@ -1153,10 +1155,12 @@ contains
if (check_for_node(node_cell, "universe")) then
call get_node_value(node_cell, "universe", c % universe)
else
c % universe = NONE
c % universe = 0
end if
if (check_for_node(node_cell, "fill")) then
call get_node_value(node_cell, "fill", c % fill)
if (find(fill_univ_ids, c % fill) == -1) &
call fill_univ_ids % push_back(c % fill)
else
c % fill = NONE
end if
@ -1352,15 +1356,16 @@ contains
! For cells, we also need to check if there's a new universe --
! also for every cell add 1 to the count of cells for the
! specified universe
universe_num = c % universe
if (.not. cells_in_univ_dict % has_key(universe_num)) then
univ_id = c % universe
if (.not. cells_in_univ_dict % has_key(univ_id)) then
n_universes = n_universes + 1
n_cells_in_univ = 1
call universe_dict % add_key(universe_num, n_universes)
call universe_dict % add_key(univ_id, n_universes)
call univ_ids % push_back(univ_id)
else
n_cells_in_univ = 1 + cells_in_univ_dict % get_key(universe_num)
n_cells_in_univ = 1 + cells_in_univ_dict % get_key(univ_id)
end if
call cells_in_univ_dict % add_key(universe_num, n_cells_in_univ)
call cells_in_univ_dict % add_key(univ_id, n_cells_in_univ)
end do
@ -1769,7 +1774,9 @@ contains
do k = 0, n_y - 1
do j = 1, n_x
lat % universes(j, n_y - k, m) = &
&temp_int_array(j + n_x*k + n_x*n_y*(m-1))
temp_int_array(j + n_x*k + n_x*n_y*(m-1))
if (find(fill_univ_ids, lat % universes(j, n_y - k, m)) == -1) &
call fill_univ_ids % push_back(lat % universes(j, n_y - k, m))
end do
end do
end do
@ -1779,6 +1786,8 @@ contains
lat % outer = NO_OUTER_UNIVERSE
if (check_for_node(node_lat, "outer")) then
call get_node_value(node_lat, "outer", lat % outer)
if (find(fill_univ_ids, lat % outer) == -1) &
call fill_univ_ids % push_back(lat % outer)
end if
! Check for 'outside' nodes which are no longer supported.
@ -1895,7 +1904,9 @@ contains
do j = 1, k
! Place universe in array.
lat % universes(i_x + n_rings, i_a + n_rings, m) = &
&temp_int_array(input_index)
temp_int_array(input_index)
if (find(fill_univ_ids, temp_int_array(input_index)) == -1) &
call fill_univ_ids % push_back(temp_int_array(input_index))
! Walk index to closest non-adjacent right neighbor.
i_x = i_x + 2
i_a = i_a - 1
@ -1920,7 +1931,9 @@ contains
do j = 1, n_rings - mod(k-1, 2)
! Place universe in array.
lat % universes(i_x + n_rings, i_a + n_rings, m) = &
&temp_int_array(input_index)
temp_int_array(input_index)
if (find(fill_univ_ids, temp_int_array(input_index)) == -1) &
call fill_univ_ids % push_back(temp_int_array(input_index))
! Walk index to closest non-adjacent right neighbor.
i_x = i_x + 2
i_a = i_a - 1
@ -1940,7 +1953,9 @@ contains
do j = 1, n_rings - k
! Place universe in array.
lat % universes(i_x + n_rings, i_a + n_rings, m) = &
&temp_int_array(input_index)
temp_int_array(input_index)
if (find(fill_univ_ids, temp_int_array(input_index)) == -1) &
call fill_univ_ids % push_back(temp_int_array(input_index))
! Walk index to closest non-adjacent right neighbor.
i_x = i_x + 2
i_a = i_a - 1
@ -1958,6 +1973,8 @@ contains
lat % outer = NO_OUTER_UNIVERSE
if (check_for_node(node_lat, "outer")) then
call get_node_value(node_lat, "outer", lat % outer)
if (find(fill_univ_ids, lat % outer) == -1) &
call fill_univ_ids % push_back(lat % outer)
end if
! Check for 'outside' nodes which are no longer supported.
@ -1974,6 +1991,47 @@ contains
end select
end do HEX_LATTICES
! ==========================================================================
! SETUP UNIVERSES
! Allocate universes, universe cell arrays, and assign base universe
allocate(universes(n_universes))
do i = 1, n_universes
associate (u => universes(i))
u % id = univ_ids % data(i)
! Allocate cell list
n_cells_in_univ = cells_in_univ_dict % get_key(u % id)
allocate(u % cells(n_cells_in_univ))
u % cells(:) = 0
! Check whether universe is a fill universe
if (find(fill_univ_ids, u % id) == -1) then
if (root_universe > 0) then
call fatal_error("Two or more universes are not used as fill &
&universes, so it is not possible to distinguish which one &
&is the root universe.")
else
root_universe = i
end if
end if
end associate
end do
do i = 1, n_cells
! Get index in universes array
j = universe_dict % get_key(cells(i) % universe)
! Set the first zero entry in the universe cells array to the index in the
! global cells array
associate (u => universes(j))
u % cells(find(u % cells, 0)) = i
end associate
end do
! Clear dictionary
call cells_in_univ_dict%clear()
! Close geometry XML file
call doc % clear()

2
src/output.F90
View file

@ -6,7 +6,7 @@ module output
use endf, only: reaction_name
use error, only: fatal_error, warning
use geometry_header, only: Cell, Universe, Lattice, RectLattice, &
HexLattice, BASE_UNIVERSE
HexLattice
use global
use math, only: t_percentile
use mesh_header, only: RegularMesh

4
src/particle_header.F90
View file

@ -4,7 +4,7 @@ module particle_header
use constants, only: NEUTRON, ONE, NONE, ZERO, MAX_SECONDARY, &
MAX_DELAYED_GROUPS, ERROR_REAL
use error, only: fatal_error
use geometry_header, only: BASE_UNIVERSE
use geometry_header, only: root_universe
implicit none
@ -144,7 +144,7 @@ contains
this % g = 1
! Set up base level coordinates
this % coord(1) % universe = BASE_UNIVERSE
this % coord(1) % universe = root_universe
this % n_coord = 1
end subroutine initialize_particle

1
src/particle_restart.F90
View file

@ -4,7 +4,6 @@ module particle_restart
use bank_header, only: Bank
use constants
use geometry_header, only: BASE_UNIVERSE
use global
use hdf5_interface, only: file_open, file_close, read_dataset
use output, only: write_message, print_particle

6
src/plot.F90
View file

@ -3,7 +3,7 @@ module plot
use constants
use error, only: fatal_error
use geometry, only: find_cell, check_cell_overlap
use geometry_header, only: Cell, BASE_UNIVERSE
use geometry_header, only: Cell, root_universe
use global
use hdf5_interface
use mesh, only: get_mesh_indices
@ -163,7 +163,7 @@ contains
call p % initialize()
p % coord(1) % xyz = xyz
p % coord(1) % uvw = [ HALF, HALF, HALF ]
p % coord(1) % universe = BASE_UNIVERSE
p % coord(1) % universe = root_universe
!$omp parallel do firstprivate(p) private(x, rgb, id) reduction(+ : data)
do y = 1, height
@ -367,7 +367,7 @@ contains
call p % initialize()
p % coord(1) % xyz = ll
p % coord(1) % uvw = [ HALF, HALF, HALF ]
p % coord(1) % universe = BASE_UNIVERSE
p % coord(1) % universe = root_universe
! Open binary plot file for writing
file_id = file_create(pl%path_plot)

1
src/source.F90
View file

@ -12,7 +12,6 @@ module source
use distribution_multivariate, only: SpatialBox
use error, only: fatal_error
use geometry, only: find_cell
use geometry_header, only: BASE_UNIVERSE
use global
use hdf5_interface, only: file_create, file_open, file_close, read_dataset
use message_passing, only: rank

55
src/summary.F90
View file

@ -4,7 +4,7 @@ module summary
use constants
use endf, only: reaction_name
use geometry_header, only: BASE_UNIVERSE, Cell, Universe, Lattice, &
use geometry_header, only: root_universe, Cell, Universe, Lattice, &
RectLattice, HexLattice
use global
use hdf5_interface
@ -123,11 +123,8 @@ contains
integer(HID_T) :: lattices_group, lattice_group
real(8), allocatable :: coeffs(:)
character(:), allocatable :: region_spec
character(MAX_LINE_LEN), allocatable :: paths(:)
character(MAX_LINE_LEN) :: path
type(Cell), pointer :: c
class(Surface), pointer :: s
type(Universe), pointer :: u
class(Lattice), pointer :: lat
! Use H5LT interface to write number of geometry objects
@ -227,22 +224,6 @@ contains
end do
call write_dataset(cell_group, "region", adjustl(region_spec))
! Write distribcell data
if (c % distribcell_index /= NONE) then
call write_dataset(cell_group, "distribcell_index", &
c % distribcell_index)
allocate(paths(c % instances))
do k = 1, c % instances
path = ''
offset = 1
call find_offset(i, universes(BASE_UNIVERSE), k, offset, path)
paths(k) = path
end do
call write_dataset(cell_group, "paths", paths)
deallocate(paths)
end if
call close_group(cell_group)
end do CELL_LOOP
@ -354,21 +335,22 @@ contains
! Write information on each universe
UNIVERSE_LOOP: do i = 1, n_universes
u => universes(i)
univ_group = create_group(universes_group, "universe " // &
trim(to_str(u%id)))
associate (u => universes(i))
univ_group = create_group(universes_group, "universe " // &
trim(to_str(u%id)))
! Write list of cells in this universe
if (u % n_cells > 0) then
allocate(cell_ids(u % n_cells))
do j = 1, u % n_cells
cell_ids(j) = cells(u % cells(j)) % id
end do
call write_dataset(univ_group, "cells", cell_ids)
deallocate(cell_ids)
end if
! Write list of cells in this universe
if (size(u % cells) > 0) then
allocate(cell_ids(size(u % cells)))
do j = 1, size(u % cells)
cell_ids(j) = cells(u % cells(j)) % id
end do
call write_dataset(univ_group, "cells", cell_ids)
deallocate(cell_ids)
end if
call close_group(univ_group)
call close_group(univ_group)
end associate
end do UNIVERSE_LOOP
call close_group(universes_group)
@ -393,13 +375,6 @@ contains
call write_dataset(lattice_group, "outer", lat % outer)
end if
! Write distribcell offsets if present
if (allocated(lat%offset)) then
if (size(lat%offset) > 0) then
call write_dataset(lattice_group, "offsets", lat%offset)
end if
end if
select type (lat)
type is (RectLattice)
! Write lattice type.

108
src/tally_filter.F90
View file

@ -3,7 +3,7 @@ module tally_filter
use algorithm, only: binary_search
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 geometry_header, only: root_universe, RectLattice, HexLattice
use global
use hdf5_interface
use mesh_header, only: RegularMesh
@ -791,7 +791,7 @@ contains
integer :: offset
type(Universe), pointer :: univ
univ => universes(BASE_UNIVERSE)
univ => universes(root_universe)
offset = 0
label = ''
call find_offset(this % cell, univ, bin-1, offset, label)
@ -1370,11 +1370,11 @@ contains
! the target cell with the given offset
!===============================================================================
recursive subroutine find_offset(goal, univ, final, offset, path)
recursive subroutine find_offset(i_cell, univ, target_offset, offset, path)
integer, intent(in) :: goal ! The target cell index
integer, intent(in) :: i_cell ! The target cell index
type(Universe), intent(in) :: univ ! Universe to begin search
integer, intent(in) :: final ! Target offset
integer, intent(in) :: target_offset ! Target offset
integer, intent(inout) :: offset ! Current offset
character(*), intent(inout) :: path ! Path to offset
@ -1387,6 +1387,7 @@ contains
integer :: cell_index ! Index in cells array
integer :: lat_offset ! Offset from lattice
integer :: temp_offset ! Looped sum of offsets
integer :: i_univ ! index in universes array
logical :: this_cell = .false. ! Advance in this cell?
logical :: later_cell = .false. ! Fill cells after this one?
type(Cell), pointer :: c ! Pointer to current cell
@ -1394,24 +1395,25 @@ contains
class(Lattice), pointer :: lat ! Pointer to current lattice
! Get the distribcell index for this cell
map = cells(goal) % distribcell_index
map = cells(i_cell) % distribcell_index
n = univ % n_cells
n = size(univ % cells)
! Write to the geometry stack
if (univ%id == 0) then
path = trim(path) // to_str(univ%id)
i_univ = universe_dict % get_key(univ % id)
if (i_univ == root_universe) then
path = trim(path) // "u" // to_str(univ%id)
else
path = trim(path) // "->" // to_str(univ%id)
path = trim(path) // "->u" // to_str(univ%id)
end if
! Look through all cells in this universe
do i = 1, n
! 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
if (univ % cells(i) == i_cell .and. offset == target_offset) then
c => cells(univ % cells(i))
path = trim(path) // "->" // to_str(c % id)
path = trim(path) // "->c" // to_str(c % id)
return
end if
end do
@ -1421,8 +1423,7 @@ contains
later_cell = .false.
cell_index = univ % cells(i)
c => cells(cell_index)
c => cells(univ % cells(i))
this_cell = .false.
@ -1433,13 +1434,10 @@ contains
do j = i+1, n
cell_index = univ % cells(j)
c => cells(cell_index)
c => cells(univ % cells(j))
! Skip normal cells which do not have offsets
if (c % type == FILL_MATERIAL) then
cycle
end if
if (c % type == FILL_MATERIAL) cycle
! Break loop once we've found the next cell with an offset
exit
@ -1463,7 +1461,7 @@ contains
! If the final offset is in the range of offset - temp_offset+offset
! then the goal is in this cell
if (final < temp_offset + offset) then
if (target_offset < temp_offset + offset) then
this_cell = .true.
end if
end if
@ -1483,7 +1481,7 @@ contains
cell_index = univ % cells(i)
c => cells(cell_index)
path = trim(path) // "->" // to_str(c%id)
path = trim(path) // "->c" // to_str(c%id)
! ====================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
@ -1493,7 +1491,7 @@ contains
offset = c % offset(map) + offset
next_univ => universes(c % fill)
call find_offset(goal, next_univ, final, offset, path)
call find_offset(i_cell, next_univ, target_offset, offset, path)
return
! ====================================================================
@ -1510,7 +1508,7 @@ contains
type is (RectLattice)
! Write to the geometry stack
path = trim(path) // "->" // to_str(lat%id)
path = trim(path) // "->l" // to_str(lat%id)
n_x = lat % n_cells(1)
n_y = lat % n_cells(2)
@ -1524,16 +1522,21 @@ contains
do l = 1, n_y
do m = 1, n_z
if (final >= lat % offset(map, k, l, m) + offset) then
if (target_offset >= lat % offset(map, k, l, m) + offset) then
if (k == n_x .and. l == n_y .and. m == n_z) then
! This is last lattice cell, so target must be here
lat_offset = lat % offset(map, k, l, m)
offset = offset + lat_offset
next_univ => universes(lat % universes(k, l, m))
path = trim(path) // "(" // trim(to_str(k)) // &
"," // trim(to_str(l)) // "," // &
trim(to_str(m)) // ")"
call find_offset(goal, next_univ, final, offset, path)
if (lat % is_3d) then
path = trim(path) // "(" // trim(to_str(k-1)) // &
"," // trim(to_str(l-1)) // "," // &
trim(to_str(m-1)) // ")"
else
path = trim(path) // "(" // trim(to_str(k-1)) // &
"," // trim(to_str(l-1)) // ")"
end if
call find_offset(i_cell, next_univ, target_offset, offset, path)
return
else
old_m = m
@ -1546,10 +1549,15 @@ contains
lat_offset = lat % offset(map, old_k, old_l, old_m)
offset = offset + lat_offset
next_univ => universes(lat % universes(old_k, old_l, old_m))
path = trim(path) // "(" // trim(to_str(old_k)) // &
"," // trim(to_str(old_l)) // "," // &
trim(to_str(old_m)) // ")"
call find_offset(goal, next_univ, final, offset, path)
if (lat % is_3d) then
path = trim(path) // "(" // trim(to_str(old_k-1)) // &
"," // trim(to_str(old_l-1)) // "," // &
trim(to_str(old_m-1)) // ")"
else
path = trim(path) // "(" // trim(to_str(old_k-1)) // &
"," // trim(to_str(old_l-1)) // ")"
end if
call find_offset(i_cell, next_univ, target_offset, offset, path)
return
end if
@ -1562,7 +1570,7 @@ contains
type is (HexLattice)
! Write to the geometry stack
path = trim(path) // "->" // to_str(lat%id)
path = trim(path) // "->l" // to_str(lat%id)
n_z = lat % n_axial
n_y = 2 * lat % n_rings - 1
@ -1584,17 +1592,23 @@ contains
cycle
end if
if (final >= lat % offset(map, k, l, m) + offset) then
if (target_offset >= lat % offset(map, k, l, m) + offset) then
if (k == lat % n_rings .and. l == n_y .and. m == n_z) then
! This is last lattice cell, so target must be here
lat_offset = lat % offset(map, k, l, m)
offset = offset + lat_offset
next_univ => universes(lat % universes(k, l, m))
path = trim(path) // "(" // &
trim(to_str(k - lat % n_rings)) // "," // &
trim(to_str(l - lat % n_rings)) // "," // &
trim(to_str(m)) // ")"
call find_offset(goal, next_univ, final, offset, path)
if (lat % is_3d) then
path = trim(path) // "(" // &
trim(to_str(k - lat % n_rings)) // "," // &
trim(to_str(l - lat % n_rings)) // "," // &
trim(to_str(m - 1)) // ")"
else
path = trim(path) // "(" // &
trim(to_str(k - lat % n_rings)) // "," // &
trim(to_str(l - lat % n_rings)) // ")"
end if
call find_offset(i_cell, next_univ, target_offset, offset, path)
return
else
old_m = m
@ -1607,11 +1621,17 @@ contains
lat_offset = lat % offset(map, old_k, old_l, old_m)
offset = offset + lat_offset
next_univ => universes(lat % universes(old_k, old_l, old_m))
path = trim(path) // "(" // &
trim(to_str(old_k - lat % n_rings)) // "," // &
trim(to_str(old_l - lat % n_rings)) // "," // &
trim(to_str(old_m)) // ")"
call find_offset(goal, next_univ, final, offset, path)
if (lat % is_3d) then
path = trim(path) // "(" // &
trim(to_str(old_k - lat % n_rings)) // "," // &
trim(to_str(old_l - lat % n_rings)) // "," // &
trim(to_str(old_m - 1)) // ")"
else
path = trim(path) // "(" // &
trim(to_str(old_k - lat % n_rings)) // "," // &
trim(to_str(old_l - lat % n_rings)) // ")"
end if
call find_offset(i_cell, next_univ, target_offset, offset, path)
return
end if

1
src/tracking.F90
View file

@ -5,7 +5,6 @@ module tracking
use error, only: fatal_error, warning
use geometry, only: find_cell, distance_to_boundary, cross_surface, &
cross_lattice, check_cell_overlap
use geometry_header, only: Universe, BASE_UNIVERSE
use global
use output, only: write_message
use message_passing

2
tests/test_asymmetric_lattice/results_true.dat
View file

@ -1 +1 @@
3a43588f986f577151b488f64d8e7ef8e568bf9b3df923d2b8267d662d176360c99c4f3788d95676861fdb0f2ef8eb7569db37b8c6f0a849de05ae96a5844e5b
f388c07481a7333bde8044d2e8805954ec41cc3e73d00f1b92b6acff03dad06e62a74e8b7595994e877f69a051ff34d8e498de69ac0efec6e396434007c10b0c

8
tests/test_asymmetric_lattice/test_asymmetric_lattice.py
View file

@ -90,14 +90,10 @@ class AsymmetricLatticeTestHarness(PyAPITestHarness):
# Extract fuel assembly lattices from the summary
cells = sp.summary.geometry.get_all_cells()
core = cells[1]
fuel = cells[80]
fuel = fuel.fill
core = core.fill
fuel_cell = cells[27]
# Append a string of lattice distribcell offsets to the string
outstr += ', '.join(map(str, fuel.offsets.flatten())) + '\n'
outstr += ', '.join(map(str, core.offsets.flatten())) + '\n'
outstr += '\n'.join(fuel_cell.paths) + '\n'
# Hash the results if necessary
if hash_output:

2
tests/test_distribmat/results_true.dat
View file

@ -7,5 +7,3 @@ Cell
Region = -10000
Rotation = None
Translation = None
Offset = None
Distribcell index= 1

2
tests/test_multipole/results_true.dat
View file

@ -8,5 +8,3 @@ Cell
Rotation = None
Temperature = [ 500. 0. 700. 800.]
Translation = None
Offset = None
Distribcell index= 1