Merge pull request #6 from smharper/python-api

Hex lattice for Python API
This commit is contained in:
Will Boyd 2015-03-16 21:28:43 -04:00
commit 842f1c8ebd
2 changed files with 466 additions and 108 deletions

View file

@ -0,0 +1,158 @@
import openmc
###############################################################################
# Simulation Input File Parameters
###############################################################################
# OpenMC simulation parameters
batches = 20
inactive = 10
particles = 10000
###############################################################################
# Exporting to OpenMC materials.xml File
###############################################################################
# Instantiate some Nuclides
h1 = openmc.Nuclide('H-1')
o16 = openmc.Nuclide('O-16')
u235 = openmc.Nuclide('U-235')
fe56 = openmc.Nuclide('Fe-56')
# Instantiate some Materials and register the appropriate Nuclides
fuel = openmc.Material(material_id=1, name='fuel')
fuel.set_density('g/cc', 4.5)
fuel.add_nuclide(u235, 1.)
moderator = openmc.Material(material_id=2, name='moderator')
moderator.set_density('g/cc', 1.0)
moderator.add_nuclide(h1, 2.)
moderator.add_nuclide(o16, 1.)
moderator.add_s_alpha_beta('HH2O', '71t')
iron = openmc.Material(material_id=3, name='iron')
iron.set_density('g/cc', 7.9)
iron.add_nuclide(fe56, 1.)
# Instantiate a MaterialsFile, register all Materials, and export to XML
materials_file = openmc.MaterialsFile()
materials_file.set_default_xs('71c')
materials_file.add_materials([moderator, fuel, iron])
materials_file.export_to_xml()
###############################################################################
# Exporting to OpenMC geometry.xml File
###############################################################################
# Instantiate Surfaces
left = openmc.XPlane(surface_id=1, x0=-3, name='left')
right = openmc.XPlane(surface_id=2, x0=3, name='right')
bottom = openmc.YPlane(surface_id=3, y0=-4, name='bottom')
top = openmc.YPlane(surface_id=4, y0=4, name='top')
fuel_surf = openmc.ZCylinder(surface_id=5, x0=0, y0=0, R=0.4)
left.set_boundary_type('vacuum')
right.set_boundary_type('vacuum')
top.set_boundary_type('vacuum')
bottom.set_boundary_type('vacuum')
# Instantiate Cells
cell1 = openmc.Cell(cell_id=1, name='Cell 1')
cell2 = openmc.Cell(cell_id=101, name='cell 2')
cell3 = openmc.Cell(cell_id=102, name='cell 3')
cell4 = openmc.Cell(cell_id=500, name='cell 4')
cell5 = openmc.Cell(cell_id=600, name='cell 5')
cell6 = openmc.Cell(cell_id=601, name='cell 6')
# Register Surfaces with Cells
cell1.add_surface(left, halfspace=+1)
cell1.add_surface(right, halfspace=-1)
cell1.add_surface(bottom, halfspace=+1)
cell1.add_surface(top, halfspace=-1)
cell2.add_surface(fuel_surf, halfspace=-1)
cell3.add_surface(fuel_surf, halfspace=+1)
cell5.add_surface(fuel_surf, halfspace=-1)
cell6.add_surface(fuel_surf, halfspace=+1)
# Register Materials with Cells
cell2.set_fill(fuel)
cell3.set_fill(moderator)
cell4.set_fill(moderator)
cell5.set_fill(iron)
cell6.set_fill(moderator)
# Instantiate Universe
univ1 = openmc.Universe(universe_id=1)
univ2 = openmc.Universe(universe_id=3)
univ3 = openmc.Universe(universe_id=4)
root = openmc.Universe(universe_id=0, name='root universe')
# Register Cells with Universe
univ1.add_cells([cell2, cell3])
univ2.add_cells([cell4])
univ3.add_cells([cell5, cell6])
root.add_cell(cell1)
# Instantiate a Lattice
lattice = openmc.HexLattice(lattice_id=5)
lattice.set_center([0., 0., 0.])
lattice.set_pitch([1., 2.])
lattice.set_universes([
[ [univ2] + [univ3]*11, [univ2] + [univ3]*5, [univ3] ],
[ [univ2] + [univ1]*11, [univ2] + [univ1]*5, [univ1] ],
[ [univ2] + [univ3]*11, [univ2] + [univ3]*5, [univ3] ]])
lattice.set_outer(univ2)
# Fill Cell with the Lattice
cell1.set_fill(lattice)
# Instantiate a Geometry and register the root Universe
geometry = openmc.Geometry()
geometry.set_root_universe(root)
# Instantiate a GeometryFile, register Geometry, and export to XML
geometry_file = openmc.GeometryFile()
geometry_file.set_geometry(geometry)
geometry_file.export_to_xml()
###############################################################################
# Exporting to OpenMC settings.xml File
###############################################################################
# Instantiate a SettingsFile, set all runtime parameters, and export to XML
settings_file = openmc.SettingsFile()
settings_file.set_batches(batches)
settings_file.set_inactive(inactive)
settings_file.set_particles(particles)
settings_file.set_source_space('box', [-1, -1, -1, 1, 1, 1])
settings_file.export_to_xml()
###############################################################################
# Exporting to OpenMC plots.xml File
###############################################################################
plot_xy = openmc.Plot(plot_id=1)
plot_xy.set_filename('plot_xy')
plot_xy.set_origin([0, 0, 0])
plot_xy.set_width([6, 6])
plot_xy.set_pixels([400, 400])
plot_xy.set_color('mat')
plot_yz = openmc.Plot(plot_id=2)
plot_yz.set_filename('plot_yz')
plot_yz.set_basis('yz')
plot_yz.set_origin([0, 0, 0])
plot_yz.set_width([8, 8])
plot_yz.set_pixels([400, 400])
plot_yz.set_color('mat')
# Instantiate a PlotsFile, add Plot, and export to XML
plot_file = openmc.PlotsFile()
plot_file.add_plot(plot_xy)
plot_file.add_plot(plot_yz)
plot_file.export_to_xml()

View file

@ -610,36 +610,6 @@ class Lattice(object):
self._name = name
def set_pitch(self, pitch):
if not isinstance(pitch, (tuple, list, np.ndarray)):
msg = 'Unable to set Lattice ID={0} pitch to {1} since ' \
'it is not a Python tuple/list or NumPy ' \
'array'.format(self._id, pitch)
raise ValueError(msg)
elif len(pitch) != 2 and len(pitch) != 3:
msg = 'Unable to set Lattice ID={0} pitch to {1} since it does ' \
'not contain 2 or 3 coordinates'.format(self._id, pitch)
raise ValueError(msg)
for dim in pitch:
if not is_integer(dim) and not is_float(dim):
msg = 'Unable to set Lattice ID={0} pitch to {1} since ' \
'it is not an an integer or floating point ' \
'value'.format(self._id, dim)
raise ValueError(msg)
elif dim < 0:
msg = 'Unable to set Lattice ID={0} pitch to {1} since it ' \
'is a negative value'.format(self._id, dim)
raise ValueError(msg)
self._pitch = pitch
def set_outer(self, outer):
if not isinstance(outer, Universe):
@ -716,37 +686,6 @@ class Lattice(object):
return all_universes
def create_xml_subelement(self, xml_element):
# Determine if XML element already contains subelement for this Lattice
path = './lattice[@id=\'{0}\']'.format(self._id)
test = xml_element.find(path)
# If the element does contain the Lattice subelement, then return
if not test is None:
return
lattice_subelement = ET.Element("lattice")
lattice_subelement.set("id", str(self._id))
# Export the Lattice cell pitch
if len(self._pitch) == 3:
pitch = ET.SubElement(lattice_subelement, "pitch")
pitch.text = '{0} {1} {2}'.format(self._pitch[0], \
self._pitch[1], \
self._pitch[2])
else:
pitch = ET.SubElement(lattice_subelement, "pitch")
pitch.text = '{0} {1}'.format(self._pitch[0], \
self._pitch[1])
# Export the Lattice outer Universe (if specified)
if self._outer is not None:
outer = ET.SubElement(lattice_subelement, "outer")
outer.text = '{0}'.format(self._outer._id)
return lattice_subelement
################################################################################
@ -834,6 +773,36 @@ class RectLattice(Lattice):
self._offsets = offsets
def set_pitch(self, pitch):
if not isinstance(pitch, (tuple, list, np.ndarray)):
msg = 'Unable to set Lattice ID={0} pitch to {1} since ' \
'it is not a Python tuple/list or NumPy ' \
'array'.format(self._id, pitch)
raise ValueError(msg)
elif len(pitch) != 2 and len(pitch) != 3:
msg = 'Unable to set Lattice ID={0} pitch to {1} since it does ' \
'not contain 2 or 3 coordinates'.format(self._id, pitch)
raise ValueError(msg)
for dim in pitch:
if not is_integer(dim) and not is_float(dim):
msg = 'Unable to set Lattice ID={0} pitch to {1} since ' \
'it is not an an integer or floating point ' \
'value'.format(self._id, dim)
raise ValueError(msg)
elif dim < 0:
msg = 'Unable to set Lattice ID={0} pitch to {1} since it ' \
'is a negative value'.format(self._id, dim)
raise ValueError(msg)
self._pitch = pitch
def get_offset(self, path, filter_offset):
# Get the current element and remove it from the list
@ -911,8 +880,24 @@ class RectLattice(Lattice):
if not test is None:
return
lattice_subelement = \
super(RectLattice, self).create_xml_subelement(xml_element)
lattice_subelement = ET.Element("lattice")
lattice_subelement.set("id", str(self._id))
# Export the Lattice cell pitch
if len(self._pitch) == 3:
pitch = ET.SubElement(lattice_subelement, "pitch")
pitch.text = '{0} {1} {2}'.format(self._pitch[0], \
self._pitch[1], \
self._pitch[2])
else:
pitch = ET.SubElement(lattice_subelement, "pitch")
pitch.text = '{0} {1}'.format(self._pitch[0], \
self._pitch[1])
# Export the Lattice outer Universe (if specified)
if self._outer is not None:
outer = ET.SubElement(lattice_subelement, "outer")
outer.text = '{0}'.format(self._outer._id)
# Export Lattice cell dimensions
if len(self._dimension) == 3:
@ -1050,6 +1035,36 @@ class HexLattice(Lattice):
self._center = center
def set_pitch(self, pitch):
if not isinstance(pitch, (tuple, list, np.ndarray)):
msg = 'Unable to set Lattice ID={0} pitch to {1} since ' \
'it is not a Python tuple/list or NumPy ' \
'array'.format(self._id, pitch)
raise ValueError(msg)
elif len(pitch) != 1 and len(pitch) != 2:
msg = 'Unable to set Lattice ID={0} pitch to {1} since it does ' \
'not contain 2 or 3 coordinates'.format(self._id, pitch)
raise ValueError(msg)
for dim in pitch:
if not is_integer(dim) and not is_float(dim):
msg = 'Unable to set Lattice ID={0} pitch to {1} since ' \
'it is not an an integer or floating point ' \
'value'.format(self._id, dim)
raise ValueError(msg)
elif dim < 0:
msg = 'Unable to set Lattice ID={0} pitch to {1} since it ' \
'is a negative value'.format(self._id, dim)
raise ValueError(msg)
self._pitch = pitch
def set_universes(self, universes):
super(HexLattice, self).set_universes(universes)
@ -1058,7 +1073,78 @@ class HexLattice(Lattice):
# The sub-lists are ordered from outermost ring to innermost ring.
# The Universes within each sub-list are ordered from the "top" in a
# clockwise fashion.
self.set_num_rings(universes.shape[0])
# Check to see if the given universes look like a 2D or a 3D array.
if isinstance(self._universes[0][0], Universe):
n_dims = 2
elif isinstance(self._universes[0][0][0], 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:
self.set_num_axial(self._universes.shape[0])
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:
self.set_num_rings(len(self._universes[0]))
for rings in self._universes:
if len(rings) != self._num_rings:
msg = 'HexLattice ID={0:d} has an inconsistent number of ' \
'rings per axial positon'.format(self._id)
raise ValueError(msg)
else:
self.set_num_rings(self._universes.shape[0])
# Make sure there are the correct number of elements in each ring.
if n_dims == 3:
for axial_slice in self._universes:
# Check the center ring.
if len(axial_slice[-1]) != 1:
msg = 'HexLattice ID={0:d} has the wrong number of ' \
'elements in the innermost ring. Only 1 element is ' \
'allowed in the innermost ring.'.format(self._id)
raise ValueError(msg)
# Check the outer rings.
for r in range(self._num_rings-1):
if len(axial_slice[r]) != 6*(self._num_rings - 1 - r):
msg = 'HexLattice ID={0:d} has the wrong number of ' \
'elements in ring number {1:d} (counting from the '\
'outermost ring). This ring should have {2:d} ' \
'elements.'.format(self._id, r,
6*(self._num_rings - 1 - r))
raise ValueError(msg)
else:
axial_slice = self._universes
# Check the center ring.
if len(axial_slice[-1]) != 1:
msg = 'HexLattice ID={0:d} has the wrong number of ' \
'elements in the innermost ring. Only 1 element is ' \
'allowed in the innermost ring.'.format(self._id)
raise ValueError(msg)
# Check the outer rings.
for r in range(self._num_rings-1):
if len(axial_slice[r]) != 6*(self._num_rings - 1 - r):
msg = 'HexLattice ID={0:d} has the wrong number of ' \
'elements in ring number {1:d} (counting from the '\
'outermost ring). This ring should have {2:d} ' \
'elements.'.format(self._id, r,
6*(self._num_rings - 1 - r))
raise ValueError(msg)
def __repr__(self):
@ -1081,16 +1167,12 @@ class HexLattice(Lattice):
string += '{0: <16}\n'.format('\tUniverses')
# FIXME: This loop must be revised for hexagonal lattice ordering
# Lattice nested Universe IDs - column major for Fortran
for i, universe in enumerate(np.ravel(self._universes)):
string += '{0} '.format(universe._id)
if self._num_axial is not None:
slices = [self._repr_axial_slice(x) for x in self._universes]
string += '\n'.join(slices)
# Add a newline character every time we reach end of row of cells
if (i+1) % self._dimension[-1] == 0:
string += '\n'
string = string.rstrip('\n')
else:
string += self._repr_axial_slice(self._universes)
return string
@ -1098,18 +1180,34 @@ class HexLattice(Lattice):
def create_xml_subelement(self, xml_element):
# Determine if XML element already contains subelement for this Lattice
path = './lattice[@id=\'{0}\']'.format(self._id)
path = './hex_lattice[@id=\'{0}\']'.format(self._id)
test = xml_element.find(path)
# If the element does contain the Lattice subelement, then return
if not test is None:
return
lattice_subelement = \
super(HexLattice, self).create_xml_subelement(xml_element)
lattice_subelement = ET.Element("hex_lattice")
lattice_subelement.set("id", str(self._id))
# Export the Lattice cell pitch
if len(self._pitch) == 2:
pitch = ET.SubElement(lattice_subelement, "pitch")
pitch.text = '{0} {1}'.format(self._pitch[0], \
self._pitch[1])
else:
pitch = ET.SubElement(lattice_subelement, "pitch")
pitch.text = '{0}'.format(self._pitch[0])
# Export the Lattice outer Universe (if specified)
if self._outer is not None:
outer = ET.SubElement(lattice_subelement, "outer")
outer.text = '{0}'.format(self._outer._id)
lattice_subelement.set("n_rings", str(self._num_rings))
lattice_subelement.set("n_axial", str(self._num_axial))
if self._num_axial is not None:
lattice_subelement.set("n_axial", str(self._num_axial))
# Export Lattice cell center
if len(self._center) == 3:
@ -1122,54 +1220,156 @@ class HexLattice(Lattice):
dimension.text = '{0} {1}'.format(self._center[0], \
self._center[1])
# FIXME: The following must be revised for hexagonal lattice ordering
# Export the Lattice nested Universe IDs - column major for Fortran
universe_ids = '\n'
# Export the Lattice nested Universe IDs.
# 3D Lattices
if len(self._dimension) == 3:
for z in range(self._dimension[2]):
for y in range(self._dimension[1]):
for x in range(self._dimension[0]):
if self._num_axial is not None:
slices = []
for z in range(self._num_axial):
# Initialize the center universe.
universe = self._universes[z][-1][0]
universe.create_xml_subelement(xml_element)
universe = self._universes[x][y][z]
# Append Universe ID to the Lattice XML subelement
universe_ids += '{0} '.format(universe._id)
# Create XML subelement for this Universe
# Initialize the remaining universes.
for r in range(self._num_rings-1):
for theta in range(6*(self._num_rings - 1 - r)):
universe = self._universes[z][r][theta]
universe.create_xml_subelement(xml_element)
# Add newline character when we reach end of row of cells
universe_ids += '\n'
# Get a string representation of the universe IDs.
slices.append(self._repr_axial_slice(self._universes[z]))
# Add newline character when we reach end of row of cells
universe_ids += '\n'
# Collapse the list of axial slices into a single string.
universe_ids = '\n'.join(slices)
# 2D Lattices
else:
for y in range(self._dimension[1]):
for x in range(self._dimension[0]):
# Initialize the center universe.
universe = self._universes[-1][0]
universe.create_xml_subelement(xml_element)
universe = self._universes[x][y]
# Append Universe ID to Lattice XML subelement
universe_ids += '{0} '.format(universe._id)
# Create XML subelement for this Universe
# Initialize the remaining universes.
for r in range(self._num_rings-1):
for theta in range(2*(self._num_rings - r)):
universe = self._universes[r][theta]
universe.create_xml_subelement(xml_element)
# Add newline character when we reach end of row of cells
universe_ids += '\n'
# Remove trailing newline character from Universe IDs string
universe_ids = universe_ids.rstrip('\n')
# Get a string representation of the universe IDs.
universe_ids = self._repr_axial_slice(self._universes)
universes = ET.SubElement(lattice_subelement, "universes")
universes.text = universe_ids
universes.text = '\n' + universe_ids
if len(self._name) > 0:
xml_element.append(ET.Comment(self._name))
# Append the XML subelement for this Lattice to the XML element
xml_element.append(lattice_subelement)
def _repr_axial_slice(self, universes):
"""Return string representation for the given 2D group of universes.
The 'universes' argument should be a list of lists of universes where
each sub-list represents a single ring. The first list should be the
outer ring.
"""
# Find the largest universe ID and count the number of digits so we can
# properly pad the output string later.
largest_id = max([max([univ._id for univ in ring])
for ring in universes])
n_digits = len(str(largest_id))
pad = ' '*n_digits
id_form = '{: ^' + str(n_digits) + 'd}'
# Initialize the list for each row.
rows = [ [] for i in range(1 + 4 * (self._num_rings-1)) ]
middle = 2 * (self._num_rings - 1)
# Start with the degenerate first ring.
universe = universes[-1][0]
rows[middle] = [id_form.format(universe._id)]
# Add universes one ring at a time.
for r in range(1, self._num_rings):
# r_prime increments down while r increments up.
r_prime = self._num_rings - 1 - r
theta = 0
y = middle + 2*r
# Climb down the top-right.
for i in range(r):
# Add the universe.
universe = universes[r_prime][theta]
rows[y].append(id_form.format(universe._id))
# Translate the indices.
y -= 1
theta += 1
# Climb down the right.
for i in range(r):
# Add the universe.
universe = universes[r_prime][theta]
rows[y].append(id_form.format(universe._id))
# Translate the indices.
y -= 2
theta += 1
# Climb down the bottom-right.
for i in range(r):
# Add the universe.
universe = universes[r_prime][theta]
rows[y].append(id_form.format(universe._id))
# Translate the indices.
y -= 1
theta += 1
# Climb up the bottom-left.
for i in range(r):
# Add the universe.
universe = universes[r_prime][theta]
rows[y].insert(0, id_form.format(universe._id))
# Translate the indices.
y += 1
theta += 1
# Climb up the left.
for i in range(r):
# Add the universe.
universe = universes[r_prime][theta]
rows[y].insert(0, id_form.format(universe._id))
# Translate the indices.
y += 2
theta += 1
# Climb up the top-left.
for i in range(r):
# Add the universe.
universe = universes[r_prime][theta]
rows[y].insert(0, id_form.format(universe._id))
# Translate the indices.
y += 1
theta += 1
# Flip the rows and join each row into a single string.
rows = [pad.join(x) for x in rows[::-1]]
# Pad the beginning of the rows so they line up properly.
for y in range(self._num_rings - 1):
rows[y] = (self._num_rings - 1 - y)*pad + rows[y]
rows[-1 - y] = (self._num_rings - 1 - y)*pad + rows[-1 - y]
for y in range(self._num_rings % 2, self._num_rings, 2):
rows[middle + y] = pad + rows[middle + y]
if y != 0: rows[middle - y] = pad + rows[middle - y]
# Join the rows together and return the string.
universe_ids = '\n'.join(rows)
return universe_ids