diff --git a/docs/source/usersguide/output/statepoint.rst b/docs/source/usersguide/output/statepoint.rst index d3c1729af6..15bc79f739 100644 --- a/docs/source/usersguide/output/statepoint.rst +++ b/docs/source/usersguide/output/statepoint.rst @@ -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 /filter /offset** (*int*) - - Filter offset (used for distribcell filter). - **/tallies/tally /filter /n_bins** (*int*) Number of bins for the j-th filter. diff --git a/docs/source/usersguide/output/summary.rst b/docs/source/usersguide/output/summary.rst index f87f60c4a0..66a5cad154 100644 --- a/docs/source/usersguide/output/summary.rst +++ b/docs/source/usersguide/output/summary.rst @@ -121,6 +121,10 @@ The current revision of the summary file format is 1. Region specification for the cell. +**/geometry/cells/cell /distribcell_index** (*int*) + + Index of this cell in distribcell filter arrays. + **/geometry/surfaces/surface /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 /moment_orders** (*char[][]*) + + Tallying moment orders for Legendre and spherical harmonic tally expansions + (*e.g.*, 'P2', 'Y1,2', etc.). + **/tallies/tally /score_bins** (*char[][]*) Scoring bins for the tally. diff --git a/openmc/filter.py b/openmc/filter.py index 04935b8edc..d61fb648f6 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -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 diff --git a/openmc/geometry.py b/openmc/geometry.py index e848e0cddf..fc8e19f071 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -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 diff --git a/openmc/opencg_compatible.py b/openmc/opencg_compatible.py index afa57c78d9..61485b5099 100644 --- a/openmc/opencg_compatible.py +++ b/openmc/opencg_compatible.py @@ -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, :] diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 845937f594..2b1e57dfa4 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -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 diff --git a/openmc/summary.py b/openmc/summary.py index 45742b0c8a..46fec01d50 100644 --- a/openmc/summary.py +++ b/openmc/summary.py @@ -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 diff --git a/openmc/tallies.py b/openmc/tallies.py index 040cbb278d..8764e2e17c 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -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. diff --git a/openmc/universe.py b/openmc/universe.py index 98367c381b..b997e8845c 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -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) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 985371ff18..d7c298979c 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -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) diff --git a/src/constants.F90 b/src/constants.F90 index aef2d61776..8eabdaf604 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -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 diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 6ca13740ba..3dad3ba395 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -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(:) diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index fc7a462e6a..01e50d983a 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -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 diff --git a/src/initialize.F90 b/src/initialize.F90 index 9d63eb4e93..30654be9d2 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -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 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index abd6e29092..07d739b889 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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) diff --git a/src/output.F90 b/src/output.F90 index c31b20b70b..66e6831371 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -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 diff --git a/src/state_point.F90 b/src/state_point.F90 index f9f4b5b757..046c857fd2 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -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. & diff --git a/src/summary.F90 b/src/summary.F90 index 217c454b81..8708d0b56f 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -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 diff --git a/src/tally.F90 b/src/tally.F90 index 09e30a2426..7210efc5c8 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -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) diff --git a/src/tally_header.F90 b/src/tally_header.F90 index 18b9219522..dcbba3d89e 100644 --- a/src/tally_header.F90 +++ b/src/tally_header.F90 @@ -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