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

Merge pull request #839 from wbinventor/openmoc-complex-regions

Browse source 
OpenMOC compatibility module now supports complex regions
This commit is contained in:
Paul Romano 2017-05-05 06:22:44 -05:00 • committed by GitHub
commit 14e4f119f3
3 changed files with 167 additions and 82 deletions
Download patch file Download diff file Expand all files Collapse all files

8
openmc/mgxs/mgxs.py
View file

@ -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:

172
openmc/openmoc_compatible.py
View file

@ -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

69
openmc/summary.py
View file

@ -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