diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index fb58efbab6..3325e0b157 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -922,13 +922,13 @@ class MGXS(object): # Override the domain object that loaded from an OpenMC summary file # NOTE: This is necessary for micro cross-sections which require # the isotopic number densities as computed by OpenMC - geom = statepoint.summary.geometry + su = statepoint.summary if self.domain_type in ('cell', 'distribcell'): - self.domain = geom.get_all_cells()[self.domain.id] + self.domain = su._fast_cells[self.domain.id] elif self.domain_type == 'universe': - self.domain = geom.get_all_universes()[self.domain.id] + self.domain = su._fast_universes[self.domain.id] elif self.domain_type == 'material': - self.domain = geom.get_all_materials()[self.domain.id] + self.domain = su._fast_materials[self.domain.id] elif self.domain_type == 'mesh': self.domain = statepoint.meshes[self.domain.id] else: diff --git a/openmc/openmoc_compatible.py b/openmc/openmoc_compatible.py index fcedc945a9..98012fae96 100644 --- a/openmc/openmoc_compatible.py +++ b/openmc/openmoc_compatible.py @@ -227,14 +227,14 @@ def get_openmc_surface(openmoc_surface): cv.check_type('openmoc_surface', openmoc_surface, openmoc.Surface) - surface_id = openmoc_surface.id + surface_id = openmoc_surface.getId() # If this Surface was already created, use it if surface_id in OPENMC_SURFACES: return OPENMC_SURFACES[surface_id] # Create an OpenMC Surface to represent this OpenMOC Surface - name = openmoc_surface.name + name = openmoc_surface.getName() # Correct for OpenMC's syntax for Surfaces dividing Cells boundary = openmoc_surface.getBoundaryType() @@ -248,6 +248,7 @@ def get_openmc_surface(openmoc_surface): boundary = 'transmission' if openmoc_surface.getSurfaceType() == openmoc.PLANE: + openmoc_surface = openmoc.castSurfaceToPlane(openmoc_surface) A = openmoc_surface.getA() B = openmoc_surface.getB() C = openmoc_surface.getC() @@ -255,21 +256,25 @@ def get_openmc_surface(openmoc_surface): openmc_surface = openmc.Plane(surface_id, boundary, A, B, C, D, name) elif openmoc_surface.getSurfaceType() == openmoc.XPLANE: + openmoc_surface = openmoc.castSurfaceToXPlane(openmoc_surface) x0 = openmoc_surface.getX() openmc_surface = openmc.XPlane(surface_id, boundary, x0, name) elif openmoc_surface.getSurfaceType() == openmoc.YPLANE: + openmoc_surface = openmoc.castSurfaceToYPlane(openmoc_surface) y0 = openmoc_surface.getY() openmc_surface = openmc.YPlane(surface_id, boundary, y0, name) elif openmoc_surface.getSurfaceType() == openmoc.ZPLANE: + openmoc_surface = openmoc.castSurfaceToZPlane(openmoc_surface) z0 = openmoc_surface.getZ() openmc_surface = openmc.ZPlane(surface_id, boundary, z0, name) elif openmoc_surface.getSurfaceType() == openmoc.ZCYLINDER: + openmoc_surface = openmoc.castSurfaceToZCylinder(openmoc_surface) x0 = openmoc_surface.getX0() y0 = openmoc_surface.getY0() - R = openmoc_surface.getR() + R = openmoc_surface.getRadius() openmc_surface = openmc.ZCylinder(surface_id, boundary, x0, y0, R, name) # Add the OpenMC Surface to the global collection of all OpenMC Surfaces @@ -324,26 +329,9 @@ def get_openmoc_cell(openmc_cell): translation = np.asarray(openmc_cell.translation, dtype=np.float64) openmoc_cell.setTranslation(translation) - # Add surfaces to OpenMOC cell from OpenMC cell region. Right now this only - # works if the region is a single half-space or an intersection of - # half-spaces, i.e., no complex cells. - region = openmc_cell.region - if region is not None: - if isinstance(region, openmc.Halfspace): - surface = region.surface - halfspace = -1 if region.side == '-' else 1 - openmoc_cell.addSurface(halfspace, get_openmoc_surface(surface)) - elif isinstance(region, openmc.Intersection): - for node in region.nodes: - if not isinstance(node, openmc.Halfspace): - raise NotImplementedError("Complex cells not yet " - "supported in OpenMOC.") - surface = node.surface - halfspace = -1 if node.side == '-' else 1 - openmoc_cell.addSurface(halfspace, get_openmoc_surface(surface)) - else: - raise NotImplementedError("Complex cells not yet supported " - "in OpenMOC.") + # Convert OpenMC's cell region to an equivalent OpenMOC region + if openmc_cell.region is not None: + openmoc_cell.setRegion(get_openmoc_region(openmc_cell.region)) # Add the OpenMC Cell to the global collection of all OpenMC Cells OPENMC_CELLS[cell_id] = openmc_cell @@ -354,6 +342,85 @@ def get_openmoc_cell(openmc_cell): return openmoc_cell +def get_openmoc_region(openmc_region): + """Return an OpenMOC region corresponding to an OpenMC region. + + Parameters + ---------- + openmc_region : openmc.Region + OpenMC region + + Returns + ------- + openmoc_region : openmoc.Region + Equivalent OpenMOC region + + """ + + cv.check_type('openmc_region', openmc_region, openmc.Region) + + # Recursively instantiate a region of the appropriate type + if isinstance(openmc_region, openmc.Halfspace): + surface = openmc_region.surface + halfspace = -1 if openmc_region.side == '-' else 1 + openmoc_region = \ + openmoc.Halfspace(halfspace, get_openmoc_surface(surface)) + elif isinstance(openmc_region, openmc.Intersection): + openmoc_region = openmoc.Intersection() + for openmc_node in openmc_region.nodes: + openmoc_region.addNode(get_openmoc_region(openmc_node)) + elif isinstance(openmc_region, openmc.Union): + openmoc_region = openmoc.Union() + for openmc_node in openmc_region.nodes: + openmoc_region.addNode(get_openmoc_region(openmc_node)) + elif isinstance(openmc_region, openmc.Complement): + openmoc_region = openmoc.Complement() + openmoc_region.addNode(get_openmoc_region(openmc_region.node)) + + return openmoc_region + + +def get_openmc_region(openmoc_region): + """Return an OpenMC region corresponding to an OpenMOC region. + + Parameters + ---------- + openmoc_region : openmoc.Region + OpenMOC region + + Returns + ------- + openmc_region : openmc.Region + Equivalent OpenMC region + + """ + + cv.check_type('openmoc_region', openmoc_region, openmoc.Region) + + # Recursively instantiate a region of the appropriate type + if openmoc_region.getRegionType() == openmoc.HALFSPACE: + openmoc_region = openmoc.castRegionToHalfspace(openmoc_region) + surface = get_openmc_surface(openmoc_region.getSurface()) + side = '-' if openmoc_region.getHalfspace() == -1 else '+' + openmc_region = openmc.Halfspace(surface, side) + elif openmoc_region.getRegionType() == openmoc.INTERSECTION: + openmc_region = openmc.Intersection() + for openmoc_node in openmoc_region.getNodes(): + openmc_node = get_openmc_region(openmoc_node) + openmc_region.nodes.append(openmc_node) + elif openmoc_region.getRegionType() == openmoc.UNION: + openmc_region = openmc.Union() + for openmoc_node in openmoc_region.getNodes(): + openmc_node = get_openmc_region(openmoc_node) + openmc_region.nodes.append(openmc_node) + elif openmoc_region.getRegionType() == openmoc.COMPLEMENT: + openmoc_nodes = openmoc_region.getNodes() + openmc_node = get_openmc_region(openmoc_nodes[0]) + openmc_region = openmc.Complement(openmc_node) + + return openmc_region + + def get_openmc_cell(openmoc_cell): """Return an OpenMC cell corresponding to an OpenMOC cell. @@ -386,24 +453,23 @@ def get_openmc_cell(openmoc_cell): openmc_cell.fill = get_openmc_material(fill) elif (openmoc_cell.getType() == openmoc.FILL): fill = openmoc_cell.getFillUniverse() - if isinstance(fill, openmoc.Lattice): + if fill.getType() == openmoc.LATTICE: + fill = openmoc.castUniverseToLattice(fill) openmc_cell.fill = get_openmc_lattice(fill) else: openmc_cell.fill = get_openmc_universe(fill) if openmoc_cell.isRotated(): - rotation = openmoc_cell.getRotation(3) + rotation = openmoc_cell.retrieveRotation(3) openmc_cell.rotation = rotation if openmoc_cell.isTranslated(): - translation = openmoc_cell.getTranslation(3) + translation = openmoc_cell.retrieveTranslation(3) openmc_cell.translation = translation - regions = [] - for surf_id, surf_halfspace in openmoc_cell.getSurfaces().values(): - halfspace = surf_halfspace._halfspace - surface = get_openmc_surface(surf_halfspace._surface) - regions.append(-surface if halfspace == -1 else +surface) - openmc_cell.region = openmc.Intersection(*regions) + # Convert OpenMC's cell region to an equivalent OpenMOC region + openmoc_region = openmoc_cell.getRegion() + if openmoc_region is not None: + openmc_cell.region = get_openmc_region(openmoc_region) # Add the OpenMC Cell to the global collection of all OpenMC Cells OPENMC_CELLS[cell_id] = openmc_cell @@ -485,7 +551,7 @@ def get_openmc_universe(openmoc_universe): openmc_universe = openmc.Universe(universe_id, name) # Convert all OpenMOC Cells in this Universe to OpenMC Cells - for openmoc_cell in openmoc_universe.getCells(): + for openmoc_cell in openmoc_universe.getCells().values(): openmc_cell = get_openmc_cell(openmoc_cell) openmc_universe.add_cell(openmc_cell) @@ -542,7 +608,8 @@ def get_openmoc_lattice(openmc_lattice): # Convert 2D lower left to 3D for OpenMOC if len(lower_left) == 2: - new_lower_left = np.ones(3, dtype=np.float64) * np.finfo(np.float64).min + new_lower_left = np.ones(3, dtype=np.float64) + new_lower_left *= np.finfo(np.float64).min / 2. new_lower_left[:2] = lower_left lower_left = new_lower_left @@ -566,7 +633,7 @@ def get_openmoc_lattice(openmc_lattice): for y in range(dimension[1]): for x in range(dimension[0]): universe_id = universes[z][y][x].id - universe_array[z][dimension[1]-y-1][x] = unique_universes[universe_id] + universe_array[z][y][x] = unique_universes[universe_id] openmoc_lattice = openmoc.Lattice(lattice_id, name) openmoc_lattice.setWidth(pitch[0], pitch[1], pitch[2]) @@ -610,15 +677,20 @@ def get_openmc_lattice(openmoc_lattice): return OPENMC_LATTICES[lattice_id] name = openmoc_lattice.getName() - dimension = [1, openmoc_lattice.getNumY(), openmoc_lattice.getNumX()] - width = [1, openmoc_lattice.getWidthY(), openmoc_lattice.getWidthX()] + dimension = [openmoc_lattice.getNumX(), + openmoc_lattice.getNumY(), + openmoc_lattice.getNumZ()] + width = [openmoc_lattice.getWidthX(), + openmoc_lattice.getWidthY(), + openmoc_lattice.getWidthZ()] offset = openmoc_lattice.getOffset() + offset = [offset.getX(), offset.getY(), offset.getZ()] lower_left = np.array(offset, dtype=np.float64) + \ ((np.array(width, dtype=np.float64) * np.array(dimension, dtype=np.float64))) / -2.0 # Initialize an empty array for the OpenMOC nested Universes in this Lattice - universe_array = np.ndarray(tuple(np.array(dimension)[::-1]), + universe_array = np.ndarray(tuple(np.array(dimension)), dtype=openmoc.Universe) # Create OpenMOC Universes for each unique nested Universe in this Lattice @@ -628,12 +700,24 @@ def get_openmc_lattice(openmoc_lattice): unique_universes[universe_id] = get_openmc_universe(universe) # Build the nested Universe array - for z in range(dimension[2]): + for x in range(dimension[0]): for y in range(dimension[1]): - for x in range(dimension[0]): - universe = openmoc_lattice.getUniverse(x, y) + for z in range(dimension[2]): + universe = openmoc_lattice.getUniverse(x, y, z) universe_id = universe.getId() - universe_array[z][dimension[1]-y-1][x] = unique_universes[universe_id] + universe_array[x][y][z] = \ + unique_universes[universe_id] + + universe_array = np.swapaxes(universe_array, 0, 1) + universe_array = universe_array[::-1, :, :] + + # Convert axially infinite 3D OpenMOC lattice to a 2D OpenMC lattice + if width[2] == np.finfo(np.float64).max: + dimension = dimension[:2] + width = width[:2] + offset = offset[:2] + lower_left = lower_left[:2] + universe_array = np.squeeze(universe_array, 2) openmc_lattice = openmc.RectLattice(lattice_id=lattice_id, name=name) openmc_lattice.pitch = width @@ -729,8 +813,8 @@ def get_openmc_geometry(openmoc_geometry): OPENMC_LATTICES.clear() OPENMOC_LATTICES.clear() - openmoc_root_universe = \ - openmoc_geometry.getRootUniverse(openmoc_root_universe) + openmoc_root_universe = openmoc_geometry.getRootUniverse() + openmc_root_universe = get_openmc_universe(openmoc_root_universe) openmc_geometry = openmc.Geometry() openmc_geometry.root_universe = openmc_root_universe diff --git a/openmc/summary.py b/openmc/summary.py index 82275007c0..61596cc38e 100644 --- a/openmc/summary.py +++ b/openmc/summary.py @@ -41,6 +41,13 @@ class Summary(object): cv.check_filetype_version(self._f, 'summary', _VERSION_SUMMARY) self._geometry = openmc.Geometry() + + self._fast_materials = {} + self._fast_surfaces = {} + self._fast_cells = {} + self._fast_universes = {} + self._fast_lattices = {} + self._materials = openmc.Materials() self._nuclides = {} @@ -62,7 +69,7 @@ class Summary(object): @property def nuclides(self): return self._nuclides - + @property def version(self): return tuple(self._f.attrs['openmc_version']) @@ -76,11 +83,11 @@ class Summary(object): def _read_geometry(self): # Read in and initialize the Materials and Geometry self._read_materials() - surfaces = self._read_surfaces() - cells, cell_fills = self._read_cells(surfaces) - universes = self._read_universes(cells) - lattices = self._read_lattices(universes) - self._finalize_geometry(cells, cell_fills, universes, lattices) + self._read_surfaces() + cell_fills = self._read_cells() + self._read_universes() + self._read_lattices() + self._finalize_geometry(cell_fills) def _read_materials(self): for group in self._f['materials'].values(): @@ -89,16 +96,15 @@ class Summary(object): # Add the material to the Materials collection self.materials.append(material) + # Store in the dictionary of materials for fast queries + self._fast_materials[material.id] = material + def _read_surfaces(self): - surfaces = {} for group in self._f['geometry/surfaces'].values(): surface = openmc.Surface.from_hdf5(group) - surfaces[surface.id] = surface + self._fast_surfaces[surface.id] = surface - return surfaces - - def _read_cells(self, surfaces): - cells = {} + def _read_cells(self): # Initialize dictionary for each Cell's fill cell_fills = {} @@ -139,29 +145,24 @@ class Summary(object): # Generate Region object given infix expression if region: - cell.region = Region.from_expression(region, surfaces) + cell.region = Region.from_expression(region, self._fast_surfaces) # Add the Cell to the global dictionary of all Cells - cells[cell.id] = cell + self._fast_cells[cell.id] = cell - return cells, cell_fills + return cell_fills - def _read_universes(self, cells): - universes = {} + def _read_universes(self): for group in self._f['geometry/universes'].values(): - universe = openmc.Universe.from_hdf5(group, cells) - universes[universe.id] = universe - return universes + universe = openmc.Universe.from_hdf5(group, self._fast_cells) + self._fast_universes[universe.id] = universe - def _read_lattices(self, universes): - lattices = {} + def _read_lattices(self): for group in self._f['geometry/lattices'].values(): - lattice = openmc.Lattice.from_hdf5(group, universes) - lattices[lattice.id] = lattice - return lattices + lattice = openmc.Lattice.from_hdf5(group, self._fast_universes) + self._fast_lattices[lattice.id] = lattice - def _finalize_geometry(self, cells, cell_fills, universes, lattices): - materials = {m.id: m for m in self.materials} + def _finalize_geometry(self, cell_fills): # 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 @@ -173,15 +174,15 @@ class Summary(object): # Retrieve the object corresponding to the fill type and ID if fill_type == 'material': if isinstance(fill_id, Iterable): - fill = [materials[mat] if mat > 0 else None + fill = [self._fast_materials[mat] if mat > 0 else None for mat in fill_id] else: - fill = materials[fill_id] if fill_id > 0 else None + fill = self._fast_materials[fill_id] if fill_id > 0 else None elif fill_type == 'universe': - fill = universes[fill_id] + fill = self._fast_universes[fill_id] fill_univ_ids.add(fill_id) else: - fill = lattices[fill_id] + fill = self._fast_lattices[fill_id] for idx in fill._natural_indices: univ = fill.get_universe(idx) fill_univ_ids.add(univ.id) @@ -189,12 +190,12 @@ class Summary(object): fill_univ_ids.add(fill.outer.id) # Set the fill for the Cell - cells[cell_id].fill = fill + self._fast_cells[cell_id].fill = fill # Determine root universe for geometry - non_fill = set(universes.keys()) - fill_univ_ids + non_fill = set(self._fast_universes.keys()) - fill_univ_ids - self.geometry.root_universe = universes[non_fill.pop()] + self.geometry.root_universe = self._fast_universes[non_fill.pop()] def add_volume_information(self, volume_calc): """Add volume information to the geometry within the summary file