mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 06:05:58 -04:00
Add geometry plotting capability and restructure lattice attributes
This commit is contained in:
Paul Romano
parent
168f1269fd
commit
e40a369693
14 changed files with 527 additions and 100 deletions
|
|
@ -260,7 +260,6 @@
|
|||
"source": [
|
||||
"# Create fuel assembly Lattice\n",
|
||||
"assembly = openmc.RectLattice(name='1.6% Fuel Assembly')\n",
|
||||
"assembly.dimension = (17, 17)\n",
|
||||
"assembly.pitch = (1.26, 1.26)\n",
|
||||
"assembly.lower_left = [-1.26 * 17. / 2.0] * 2"
|
||||
]
|
||||
|
|
@ -1597,21 +1596,21 @@
|
|||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python 2",
|
||||
"language": "python2",
|
||||
"name": "python2"
|
||||
"display_name": "Python 3",
|
||||
"language": "python",
|
||||
"name": "python3"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
"name": "ipython",
|
||||
"version": 2
|
||||
"version": 3
|
||||
},
|
||||
"file_extension": ".py",
|
||||
"mimetype": "text/x-python",
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython2",
|
||||
"version": "2.7.11"
|
||||
"pygments_lexer": "ipython3",
|
||||
"version": "3.5.1"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
|
|
|
|||
|
|
@ -244,7 +244,6 @@
|
|||
"source": [
|
||||
"# Create fuel assembly Lattice\n",
|
||||
"assembly = openmc.RectLattice(name='1.6% Fuel Assembly')\n",
|
||||
"assembly.dimension = (17, 17)\n",
|
||||
"assembly.pitch = (1.26, 1.26)\n",
|
||||
"assembly.lower_left = [-1.26 * 17. / 2.0] * 2"
|
||||
]
|
||||
|
|
|
|||
|
|
@ -199,7 +199,6 @@
|
|||
"source": [
|
||||
"# Create fuel assembly Lattice\n",
|
||||
"assembly = openmc.RectLattice(name='1.6% Fuel - 0BA')\n",
|
||||
"assembly.dimension = (17, 17)\n",
|
||||
"assembly.pitch = (1.26, 1.26)\n",
|
||||
"assembly.lower_left = [-1.26 * 17. / 2.0] * 2\n",
|
||||
"assembly.universes = [[pin_cell_universe] * 17] * 17"
|
||||
|
|
@ -2194,21 +2193,21 @@
|
|||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python 2",
|
||||
"display_name": "Python 3",
|
||||
"language": "python",
|
||||
"name": "python2"
|
||||
"name": "python3"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
"name": "ipython",
|
||||
"version": 2
|
||||
"version": 3
|
||||
},
|
||||
"file_extension": ".py",
|
||||
"mimetype": "text/x-python",
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython2",
|
||||
"version": "2.7.6"
|
||||
"pygments_lexer": "ipython3",
|
||||
"version": "3.5.1"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
|
|
|
|||
|
|
@ -98,14 +98,12 @@ univ4.add_cell(cell2)
|
|||
|
||||
# Instantiate nested Lattices
|
||||
lattice1 = openmc.RectLattice(lattice_id=4, name='4x4 assembly')
|
||||
lattice1.dimension = [2, 2]
|
||||
lattice1.lower_left = [-1., -1.]
|
||||
lattice1.pitch = [1., 1.]
|
||||
lattice1.universes = [[univ1, univ2],
|
||||
[univ2, univ3]]
|
||||
|
||||
lattice2 = openmc.RectLattice(lattice_id=6, name='4x4 core')
|
||||
lattice2.dimension = [2, 2]
|
||||
lattice2.lower_left = [-2., -2.]
|
||||
lattice2.pitch = [2., 2.]
|
||||
lattice2.universes = [[univ4, univ4],
|
||||
|
|
|
|||
|
|
@ -94,7 +94,6 @@ root.add_cell(cell1)
|
|||
|
||||
# Instantiate a Lattice
|
||||
lattice = openmc.RectLattice(lattice_id=5)
|
||||
lattice.dimension = [4, 4]
|
||||
lattice.lower_left = [-2., -2.]
|
||||
lattice.pitch = [1., 1.]
|
||||
lattice.universes = [[univ1, univ2, univ1, univ2],
|
||||
|
|
|
|||
|
|
@ -97,7 +97,10 @@ class Cell(object):
|
|||
self.region = region
|
||||
|
||||
def __contains__(self, point):
|
||||
return point in self.region
|
||||
if self.region is None:
|
||||
return True
|
||||
else:
|
||||
return point in self.region
|
||||
|
||||
def __eq__(self, other):
|
||||
if not isinstance(other, Cell):
|
||||
|
|
|
|||
|
|
@ -62,6 +62,23 @@ class Geometry(object):
|
|||
tree.write("geometry.xml", xml_declaration=True, encoding='utf-8',
|
||||
method="xml")
|
||||
|
||||
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
|
||||
|
||||
"""
|
||||
return self.root_universe.find(point)
|
||||
|
||||
def get_cell_instance(self, path):
|
||||
"""Return the instance number for the final cell in a geometry path.
|
||||
|
||||
|
|
|
|||
|
|
@ -1,8 +1,12 @@
|
|||
from __future__ import division
|
||||
|
||||
import abc
|
||||
from collections import OrderedDict, Iterable
|
||||
from math import sqrt, floor
|
||||
from numbers import Real, Integral
|
||||
from xml.etree import ElementTree as ET
|
||||
import sys
|
||||
import warnings
|
||||
|
||||
import numpy as np
|
||||
|
||||
|
|
@ -113,12 +117,6 @@ class Lattice(object):
|
|||
cv.check_type('outer universe', outer, openmc.Universe)
|
||||
self._outer = outer
|
||||
|
||||
@universes.setter
|
||||
def universes(self, universes):
|
||||
cv.check_iterable_type('lattice universes', universes, openmc.Universe,
|
||||
min_depth=2, max_depth=3)
|
||||
self._universes = np.asarray(universes)
|
||||
|
||||
def get_unique_universes(self):
|
||||
"""Determine all unique universes in the lattice
|
||||
|
||||
|
|
@ -239,6 +237,11 @@ class Lattice(object):
|
|||
class RectLattice(Lattice):
|
||||
"""A lattice consisting of rectangular prisms.
|
||||
|
||||
To completely define a rectangular lattice, the
|
||||
:attr:`RectLattice.lower_left` :attr:`RectLattice.pitch`,
|
||||
:attr:`RectLattice.outer`, and :attr:`RectLattice.universes` properties need
|
||||
to be set.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
lattice_id : int, optional
|
||||
|
|
@ -253,12 +256,6 @@ class RectLattice(Lattice):
|
|||
Unique identifier for the lattice
|
||||
name : str
|
||||
Name of the lattice
|
||||
dimension : Iterable of int
|
||||
An array of two or three integers representing the number of lattice
|
||||
cells in the x- and y- (and z-) directions, respectively.
|
||||
lower_left : Iterable of float
|
||||
The coordinates of the lower-left corner of the lattice. If the lattice
|
||||
is two-dimensional, only the x- and y-coordinates are specified.
|
||||
pitch : Iterable of float
|
||||
Pitch of the lattice in the x, y, and (if applicable) z directions in
|
||||
cm.
|
||||
|
|
@ -266,7 +263,25 @@ class RectLattice(Lattice):
|
|||
A universe to fill all space outside the lattice
|
||||
universes : Iterable of Iterable of openmc.Universe
|
||||
A two- or three-dimensional list/array of universes filling each element
|
||||
of the lattice
|
||||
of the lattice. The first dimension corresponds to the z-direction (if
|
||||
applicable), the second dimension corresponds to the y-direction, and
|
||||
the third dimension corresponds to the x-direction. Note that for the
|
||||
y-direction, a higher index corresponds to a lower physical
|
||||
y-value. Each z-slice in the array can be thought of as a top-down view
|
||||
of the lattice.
|
||||
lower_left : Iterable of float
|
||||
The Cartesian coordinates of the lower-left corner of the lattice. If
|
||||
the lattice is two-dimensional, only the x- and y-coordinates are
|
||||
specified.
|
||||
indices : list of tuple
|
||||
A list of all possible (z,y,x) or (y,x) lattice element indices. These
|
||||
indices correspond to indices in the :attr:`RectLattice.universes`
|
||||
property.
|
||||
ndim : int
|
||||
The number of dimensions of the lattice
|
||||
shape : Iterable of int
|
||||
An array of two or three integers representing the number of lattice
|
||||
cells in the x- and y- (and z-) directions, respectively.
|
||||
|
||||
"""
|
||||
|
||||
|
|
@ -274,7 +289,6 @@ class RectLattice(Lattice):
|
|||
super(RectLattice, self).__init__(lattice_id, name)
|
||||
|
||||
# Initialize Lattice class attributes
|
||||
self._dimension = None
|
||||
self._lower_left = None
|
||||
self._offsets = None
|
||||
|
||||
|
|
@ -283,7 +297,7 @@ class RectLattice(Lattice):
|
|||
return False
|
||||
elif not super(RectLattice, self).__eq__(other):
|
||||
return False
|
||||
elif self.dimension != other.dimension:
|
||||
elif self.shape != other.shape:
|
||||
return False
|
||||
elif self.lower_left != other.lower_left:
|
||||
return False
|
||||
|
|
@ -300,8 +314,8 @@ class RectLattice(Lattice):
|
|||
string = 'RectLattice\n'
|
||||
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
|
||||
string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name)
|
||||
string += '{0: <16}{1}{2}\n'.format('\tDimension', '=\t',
|
||||
self._dimension)
|
||||
string += '{0: <16}{1}{2}\n'.format('\tShape', '=\t',
|
||||
self.shape)
|
||||
string += '{0: <16}{1}{2}\n'.format('\tLower Left', '=\t',
|
||||
self._lower_left)
|
||||
string += '{0: <16}{1}{2}\n'.format('\tPitch', '=\t', self._pitch)
|
||||
|
|
@ -320,7 +334,7 @@ class RectLattice(Lattice):
|
|||
string += '{0} '.format(universe._id)
|
||||
|
||||
# Add a newline character every time we reach end of row of cells
|
||||
if (i+1) % self._dimension[-1] == 0:
|
||||
if (i+1) % self.shape[0] == 0:
|
||||
string += '\n'
|
||||
|
||||
string = string.rstrip('\n')
|
||||
|
|
@ -333,7 +347,7 @@ class RectLattice(Lattice):
|
|||
string += '{0} '.format(offset)
|
||||
|
||||
# Add a newline character when we reach end of row of cells
|
||||
if (i+1) % self._dimension[-1] == 0:
|
||||
if (i+1) % self.shape[0] == 0:
|
||||
string += '\n'
|
||||
|
||||
string = string.rstrip('\n')
|
||||
|
|
@ -341,24 +355,29 @@ class RectLattice(Lattice):
|
|||
return string
|
||||
|
||||
@property
|
||||
def dimension(self):
|
||||
return self._dimension
|
||||
def indices(self):
|
||||
if self.ndim == 2:
|
||||
return list(np.broadcast(*np.ogrid[
|
||||
:self.shape[1], :self.shape[0]]))
|
||||
else:
|
||||
return list(np.broadcast(*np.ogrid[
|
||||
:self.shape[2], :self.shape[1], :self.shape[0]]))
|
||||
|
||||
@property
|
||||
def lower_left(self):
|
||||
return self._lower_left
|
||||
|
||||
@property
|
||||
def ndim(self):
|
||||
return len(self.pitch)
|
||||
|
||||
@property
|
||||
def offsets(self):
|
||||
return self._offsets
|
||||
|
||||
@dimension.setter
|
||||
def dimension(self, dimension):
|
||||
cv.check_type('lattice dimension', dimension, Iterable, Integral)
|
||||
cv.check_length('lattice dimension', dimension, 2, 3)
|
||||
for dim in dimension:
|
||||
cv.check_greater_than('lattice dimension', dim, 0)
|
||||
self._dimension = dimension
|
||||
@property
|
||||
def shape(self):
|
||||
return self._universes.shape[::-1]
|
||||
|
||||
@lower_left.setter
|
||||
def lower_left(self, lower_left):
|
||||
|
|
@ -379,8 +398,13 @@ class RectLattice(Lattice):
|
|||
cv.check_greater_than('lattice pitch', dim, 0.0)
|
||||
self._pitch = pitch
|
||||
|
||||
def get_cell_instance(self, path, distribcell_index):
|
||||
@Lattice.universes.setter
|
||||
def universes(self, universes):
|
||||
cv.check_iterable_type('lattice universes', universes, openmc.Universe,
|
||||
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]
|
||||
|
|
@ -395,7 +419,7 @@ class RectLattice(Lattice):
|
|||
lat_z = int(i.split(',')[2]) - 1
|
||||
|
||||
# For 2D Lattices
|
||||
if len(self._dimension) == 2:
|
||||
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)
|
||||
|
|
@ -408,6 +432,128 @@ class RectLattice(Lattice):
|
|||
|
||||
return offset
|
||||
|
||||
def find_element(self, point):
|
||||
"""Determine index of lattice element and local coordinates for a point
|
||||
|
||||
Parameters
|
||||
----------
|
||||
point : Iterable of float
|
||||
Cartesian coordinates of point
|
||||
|
||||
Returns
|
||||
-------
|
||||
2- or 3-tuple of int
|
||||
A tuple of the corresponding (x,y,z) lattice element indices
|
||||
3-tuple of float
|
||||
Carestian coordinates of the point in the corresponding lattice
|
||||
element coordinate system
|
||||
|
||||
"""
|
||||
ix = floor((point[0] - self._lower_left[0])/self._pitch[0])
|
||||
iy = floor((point[1] - self._lower_left[1])/self._pitch[1])
|
||||
if self.ndim == 2:
|
||||
idx = (ix, iy)
|
||||
else:
|
||||
iz = floor((point[2] - self._lower_left[2])/self._pitch[2])
|
||||
idx = (ix, iy, iz)
|
||||
return idx, self.get_local_coordinates(point, idx)
|
||||
|
||||
def get_local_coordinates(self, point, idx):
|
||||
"""Determine local coordinates of a point within a lattice element
|
||||
|
||||
Parameters
|
||||
----------
|
||||
point : Iterable of float
|
||||
Cartesian coordinates of point
|
||||
idx : Iterable of int
|
||||
(x,y,z) indices of lattice element. If the lattice is 2D, the z
|
||||
index can be omitted.
|
||||
|
||||
Returns
|
||||
-------
|
||||
3-tuple of float
|
||||
Cartesian coordinates of point in the lattice element coordinate
|
||||
system
|
||||
|
||||
"""
|
||||
x = point[0] - (self._lower_left[0] + (idx[0] + 0.5)*self._pitch[0])
|
||||
y = point[1] - (self._lower_left[1] + (idx[1] + 0.5)*self._pitch[1])
|
||||
if self.ndim == 2:
|
||||
z = point[2]
|
||||
else:
|
||||
z = point[2] - (self._lower_left[2] + (idx[2] + 0.5)*self._pitch[2])
|
||||
return (x, y, z)
|
||||
|
||||
def get_universe_index(self, idx):
|
||||
"""Return index in the universes array corresponding to a lattice element index
|
||||
|
||||
Parameters
|
||||
----------
|
||||
idx : Iterable of int
|
||||
Lattice element indices in the :math:`(x,y,z)` coordinate system
|
||||
|
||||
Returns
|
||||
-------
|
||||
2- or 3-tuple of int
|
||||
Indices used when setting the :attr:`RectLattice.universes` property
|
||||
|
||||
"""
|
||||
max_y = self.shape[1] - 1
|
||||
if self.ndim == 2:
|
||||
x, y = idx
|
||||
return (max_y - y, x)
|
||||
else:
|
||||
x, y, z = idx
|
||||
return (z, max_y - y, x)
|
||||
|
||||
def is_valid_index(self, idx):
|
||||
"""Determine whether lattice element index is within defined range
|
||||
|
||||
Parameters
|
||||
----------
|
||||
idx : Iterable of int
|
||||
Lattice element indices in the :math:`(x,y,z)` coordinate system
|
||||
|
||||
Returns
|
||||
-------
|
||||
bool
|
||||
Whether index is valid
|
||||
|
||||
"""
|
||||
if self.ndim == 2:
|
||||
return (0 <= idx[0] < self.shape[0] and
|
||||
0 <= idx[1] < self.shape[1])
|
||||
else:
|
||||
return (0 <= idx[0] < self.shape[0] and
|
||||
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
|
||||
|
|
@ -436,7 +582,7 @@ class RectLattice(Lattice):
|
|||
|
||||
# Export Lattice cell dimensions
|
||||
dimension = ET.SubElement(lattice_subelement, "dimension")
|
||||
dimension.text = ' '.join(map(str, self._dimension))
|
||||
dimension.text = ' '.join(map(str, self.shape))
|
||||
|
||||
# Export Lattice lower left
|
||||
lower_left = ET.SubElement(lattice_subelement, "lower_left")
|
||||
|
|
@ -446,10 +592,10 @@ class RectLattice(Lattice):
|
|||
universe_ids = '\n'
|
||||
|
||||
# 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.ndim == 3:
|
||||
for z in range(self.shape[2]):
|
||||
for y in range(self.shape[1]):
|
||||
for x in range(self.shape[0]):
|
||||
universe = self._universes[z][y][x]
|
||||
|
||||
# Append Universe ID to the Lattice XML subelement
|
||||
|
|
@ -466,8 +612,8 @@ class RectLattice(Lattice):
|
|||
|
||||
# 2D Lattices
|
||||
else:
|
||||
for y in range(self._dimension[1]):
|
||||
for x in range(self._dimension[0]):
|
||||
for y in range(self.shape[1]):
|
||||
for x in range(self.shape[0]):
|
||||
universe = self._universes[y][x]
|
||||
|
||||
# Append Universe ID to Lattice XML subelement
|
||||
|
|
@ -492,6 +638,10 @@ class RectLattice(Lattice):
|
|||
class HexLattice(Lattice):
|
||||
"""A lattice consisting of hexagonal prisms.
|
||||
|
||||
To completely define a hexagonal lattice, the :attr:`HexLattice.center`,
|
||||
:attr:`HexLattice.pitch`, :attr:`HexLattice.universes`, and
|
||||
:attr:`HexLattice.outer` properties need to be set.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
lattice_id : int, optional
|
||||
|
|
@ -506,26 +656,31 @@ class HexLattice(Lattice):
|
|||
Unique identifier for the lattice
|
||||
name : str
|
||||
Name of the lattice
|
||||
num_rings : int
|
||||
Number of radial ring positions in the xy-plane
|
||||
num_axial : int
|
||||
Number of positions along the z-axis.
|
||||
center : Iterable of float
|
||||
Coordinates of the center of the lattice. If the lattice does not have
|
||||
axial sections then only the x- and y-coordinates are specified
|
||||
pitch : Iterable of float
|
||||
Pitch of the lattice in cm. The first item in the iterable specifies the
|
||||
pitch in the radial direction and, if the lattice is 3D, the second item
|
||||
in the iterable specifies the pitch in the axial direction.
|
||||
outer : openmc.Universe
|
||||
A universe to fill all space outside the lattice
|
||||
universes : Iterable of Iterable of openmc.Universe
|
||||
universes : Nested Iterable of openmc.Universe
|
||||
A two- or three-dimensional list/array of universes filling each element
|
||||
of the lattice. Each sub-list corresponds to one ring of universes and
|
||||
should be ordered from outermost ring to innermost ring. The universes
|
||||
within each sub-list are ordered from the "top" and proceed in a
|
||||
clockwise fashion. The :meth:`HexLattice.show_indices` method can be
|
||||
used to help figure out indices for this property.
|
||||
center : Iterable of float
|
||||
Coordinates of the center of the lattice. If the lattice does not have
|
||||
axial sections then only the x- and y-coordinates are specified
|
||||
indices : list of tuple
|
||||
A list of all possible (z,r,i) or (r,i) lattice element indices that are
|
||||
possible, where z is the axial index, r is in the ring index (starting
|
||||
from the outermost ring), and i is the index with a ring starting from
|
||||
the top and proceeding clockwise.
|
||||
num_rings : int
|
||||
Number of radial ring positions in the xy-plane
|
||||
num_axial : int
|
||||
Number of positions along the z-axis.
|
||||
|
||||
"""
|
||||
|
||||
|
|
@ -597,17 +752,15 @@ class HexLattice(Lattice):
|
|||
def center(self):
|
||||
return self._center
|
||||
|
||||
@num_rings.setter
|
||||
def num_rings(self, num_rings):
|
||||
cv.check_type('number of rings', num_rings, Integral)
|
||||
cv.check_greater_than('number of rings', num_rings, 0)
|
||||
self._num_rings = num_rings
|
||||
|
||||
@num_axial.setter
|
||||
def num_axial(self, num_axial):
|
||||
cv.check_type('number of axial', num_axial, Integral)
|
||||
cv.check_greater_than('number of axial', num_axial, 0)
|
||||
self._num_axial = num_axial
|
||||
@property
|
||||
def indices(self):
|
||||
if self.num_axial is None:
|
||||
return [(r, i) for r in range(self._num_rings)
|
||||
for i in range(max(6*(self._num_rings - 1 - r), 1))]
|
||||
else:
|
||||
return [(z, r, i) for z in range(self._num_axial)
|
||||
for r in range(self._num_rings)
|
||||
for i in range(max(6*(self._num_rings - 1 - r), 1))]
|
||||
|
||||
@center.setter
|
||||
def center(self, center):
|
||||
|
|
@ -625,8 +778,9 @@ class HexLattice(Lattice):
|
|||
|
||||
@Lattice.universes.setter
|
||||
def universes(self, universes):
|
||||
# Call Lattice.universes parent class setter property
|
||||
Lattice.universes.fset(self, universes)
|
||||
cv.check_iterable_type('lattice universes', universes, openmc.Universe,
|
||||
min_depth=2, max_depth=3)
|
||||
self._universes = universes
|
||||
|
||||
# NOTE: This routine assumes that the user creates a "ragged" list of
|
||||
# lists, where each sub-list corresponds to one ring of Universes.
|
||||
|
|
@ -649,14 +803,14 @@ class HexLattice(Lattice):
|
|||
|
||||
# Set the number of axial positions.
|
||||
if n_dims == 3:
|
||||
self.num_axial = len(self._universes)
|
||||
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:
|
||||
self.num_rings = len(self._universes[0])
|
||||
self._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 ' \
|
||||
|
|
@ -664,7 +818,7 @@ class HexLattice(Lattice):
|
|||
raise ValueError(msg)
|
||||
|
||||
else:
|
||||
self.num_rings = len(self._universes)
|
||||
self._num_rings = len(self._universes)
|
||||
|
||||
# Make sure there are the correct number of elements in each ring.
|
||||
if n_dims == 3:
|
||||
|
|
@ -705,6 +859,170 @@ class HexLattice(Lattice):
|
|||
6*(self._num_rings - 1 - r))
|
||||
raise ValueError(msg)
|
||||
|
||||
def find_element(self, point):
|
||||
"""Determine index of lattice element and local coordinates for a point
|
||||
|
||||
Parameters
|
||||
----------
|
||||
point : Iterable of float
|
||||
Cartesian coordinates of point
|
||||
|
||||
Returns
|
||||
-------
|
||||
3-tuple of int
|
||||
Indices of corresponding lattice element in (x,:math:`alpha`,z)
|
||||
bases
|
||||
numpy.ndarray
|
||||
Carestian coordinates of the point in the corresponding lattice
|
||||
element coordinate system
|
||||
|
||||
"""
|
||||
# Convert coordinates to skewed bases
|
||||
x = point[0] - self._center[0]
|
||||
y = point[1] - self._center[1]
|
||||
if self._num_axial is None:
|
||||
iz = 1
|
||||
else:
|
||||
z = point[2] - self._center[2]
|
||||
iz = floor(z/self._pitch[1] + 0.5*self._num_axial)
|
||||
alpha = y - x/sqrt(3.)
|
||||
ix = floor(x/(sqrt(0.75) * self._pitch[0]))
|
||||
ia = floor(alpha/self._pitch[0])
|
||||
|
||||
# Check four lattice elements to see which one is closest based on local
|
||||
# coordinates
|
||||
d_min = np.inf
|
||||
for idx in [(ix, ia, iz), (ix + 1, ia, iz), (ix, ia + 1, iz),
|
||||
(ix + 1, ia + 1, iz)]:
|
||||
p = self.get_local_coordinates(point, idx)
|
||||
d = p[0]**2 + p[1]**2
|
||||
if d < d_min:
|
||||
d_min = d
|
||||
idx_min = idx
|
||||
p_min = p
|
||||
|
||||
return idx_min, p_min
|
||||
|
||||
def get_local_coordinates(self, point, idx):
|
||||
"""Determine local coordinates of a point within a lattice element
|
||||
|
||||
Parameters
|
||||
----------
|
||||
point : Iterable of float
|
||||
Cartesian coordinates of point
|
||||
idx : Iterable of int
|
||||
Indices of lattice element in (x,:math:`alpha`,z) bases
|
||||
|
||||
Returns
|
||||
-------
|
||||
3-tuple of float
|
||||
Cartesian coordinates of point in the lattice element coordinate
|
||||
system
|
||||
|
||||
"""
|
||||
x = point[0] - (self._center[0] + sqrt(0.75)*self._pitch[0]*idx[0])
|
||||
y = point[1] - (self._center[1] + (0.5*idx[0] + idx[1])*self._pitch[0])
|
||||
if self._num_axial is None:
|
||||
z = point[2]
|
||||
else:
|
||||
z = point[2] - (self._center[2] + (idx[2] + 0.5 - 0.5*self._num_axial)*
|
||||
self._pitch[1])
|
||||
return (x, y, z)
|
||||
|
||||
def get_universe_index(self, idx):
|
||||
"""Return index in the universes array corresponding to a lattice element index
|
||||
|
||||
Parameters
|
||||
----------
|
||||
idx : Iterable of int
|
||||
Lattice element indices in the :math:`(x,\alpha,z)` coordinate
|
||||
system
|
||||
|
||||
Returns
|
||||
-------
|
||||
2- or 3-tuple of int
|
||||
Indices used when setting the :attr:`HexLattice.universes` property
|
||||
|
||||
"""
|
||||
|
||||
# First we determine which ring the index corresponds to.
|
||||
x = idx[0]
|
||||
a = idx[1]
|
||||
z = -a - x
|
||||
g = max(abs(x), abs(a), abs(z))
|
||||
|
||||
# Next we use a clever method to figure out where along the ring we are.
|
||||
i_ring = self._num_rings - 1 - g
|
||||
if x >= 0:
|
||||
if a >= 0:
|
||||
i_within = x
|
||||
else:
|
||||
i_within = 2*g + z
|
||||
else:
|
||||
if a <= 0:
|
||||
i_within = 3*g - x
|
||||
else:
|
||||
i_within = 5*g - z
|
||||
|
||||
if self.num_axial is None:
|
||||
return (i_ring, i_within)
|
||||
else:
|
||||
return (idx[2], i_ring, i_within)
|
||||
|
||||
def is_valid_index(self, idx):
|
||||
"""Determine whether lattice element index is within defined range
|
||||
|
||||
Parameters
|
||||
----------
|
||||
idx : Iterable of int
|
||||
Lattice element indices in the :math:`(x,\alpha,z)` coordinate
|
||||
system
|
||||
|
||||
Returns
|
||||
-------
|
||||
bool
|
||||
Whether index is valid
|
||||
|
||||
"""
|
||||
x = idx[0]
|
||||
y = idx[1]
|
||||
z = 0 - y - x
|
||||
g = max(abs(x), abs(y), abs(z))
|
||||
if self.num_axial is None:
|
||||
return g < self.num_rings
|
||||
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)
|
||||
|
|
@ -736,8 +1054,8 @@ class HexLattice(Lattice):
|
|||
lattice_subelement.set("n_axial", str(self._num_axial))
|
||||
|
||||
# Export Lattice cell center
|
||||
dimension = ET.SubElement(lattice_subelement, "center")
|
||||
dimension.text = ' '.join(map(str, self._center))
|
||||
center = ET.SubElement(lattice_subelement, "center")
|
||||
center.text = ' '.join(map(str, self._center))
|
||||
|
||||
# Export the Lattice nested Universe IDs.
|
||||
|
||||
|
|
|
|||
|
|
@ -861,8 +861,8 @@ def get_opencg_lattice(openmc_lattice):
|
|||
universes = new_universes
|
||||
|
||||
# Initialize an empty array for the OpenCG nested Universes in this Lattice
|
||||
universe_array = np.ndarray(tuple(np.array(dimension)[::-1]),
|
||||
dtype=opencg.Universe)
|
||||
universe_array = np.empty(tuple(np.array(dimension)[::-1]),
|
||||
dtype=opencg.Universe)
|
||||
|
||||
# Create OpenCG Universes for each unique nested Universe in this Lattice
|
||||
unique_universes = openmc_lattice.get_unique_universes()
|
||||
|
|
@ -929,8 +929,8 @@ 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)[::-1]),
|
||||
dtype=openmc.Universe)
|
||||
universe_array = np.empty(tuple(np.array(dimension)[::-1]),
|
||||
dtype=openmc.Universe)
|
||||
|
||||
# Create OpenMC Universes for each unique nested Universe in this Lattice
|
||||
unique_universes = opencg_lattice.get_unique_universes()
|
||||
|
|
@ -953,7 +953,6 @@ def get_openmc_lattice(opencg_lattice):
|
|||
np.array(dimension, dtype=np.float64))) / -2.0
|
||||
|
||||
openmc_lattice = openmc.RectLattice(lattice_id=lattice_id)
|
||||
openmc_lattice.dimension = dimension
|
||||
openmc_lattice.pitch = width
|
||||
openmc_lattice.universes = universe_array
|
||||
openmc_lattice.lower_left = lower_left
|
||||
|
|
|
|||
|
|
@ -358,7 +358,6 @@ class Summary(object):
|
|||
|
||||
# Create the Lattice
|
||||
lattice = openmc.RectLattice(lattice_id=lattice_id, name=name)
|
||||
lattice.dimension = tuple(dimension)
|
||||
lattice.lower_left = lower_left
|
||||
lattice.pitch = pitch
|
||||
|
||||
|
|
@ -368,7 +367,7 @@ class Summary(object):
|
|||
|
||||
# Build array of Universe pointers for the Lattice
|
||||
universes = \
|
||||
np.ndarray(tuple(universe_ids.shape), dtype=openmc.Universe)
|
||||
np.empty(tuple(universe_ids.shape), dtype=openmc.Universe)
|
||||
|
||||
for z in range(universe_ids.shape[0]):
|
||||
for y in range(universe_ids.shape[1]):
|
||||
|
|
@ -403,8 +402,6 @@ class Summary(object):
|
|||
|
||||
# Create the Lattice
|
||||
lattice = openmc.HexLattice(lattice_id=lattice_id, name=name)
|
||||
lattice.num_rings = n_rings
|
||||
lattice.num_axial = n_axial
|
||||
lattice.center = center
|
||||
lattice.pitch = pitch
|
||||
|
||||
|
|
@ -417,12 +414,12 @@ class Summary(object):
|
|||
# (x, alpha, z) to the Python API's format of a ragged nested
|
||||
# list of (z, ring, theta).
|
||||
universes = []
|
||||
for z in range(lattice.num_axial):
|
||||
for z in range(n_axial):
|
||||
# Add a list for this axial level.
|
||||
universes.append([])
|
||||
x = lattice.num_rings - 1
|
||||
a = 2*lattice.num_rings - 2
|
||||
for r in range(lattice.num_rings - 1, 0, -1):
|
||||
x = n_rings - 1
|
||||
a = 2*n_rings - 2
|
||||
for r in range(n_rings - 1, 0, -1):
|
||||
# Add a list for this ring.
|
||||
universes[-1].append([])
|
||||
|
||||
|
|
|
|||
|
|
@ -1,6 +1,7 @@
|
|||
from collections import OrderedDict, Iterable
|
||||
from numbers import Integral
|
||||
from xml.etree import ElementTree as ET
|
||||
import random
|
||||
import sys
|
||||
import warnings
|
||||
|
||||
|
|
@ -124,6 +125,110 @@ class Universe(object):
|
|||
else:
|
||||
self._name = ''
|
||||
|
||||
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
|
||||
|
||||
"""
|
||||
p = np.asarray(point)
|
||||
for cell in self._cells.values():
|
||||
if p in cell:
|
||||
if cell._type in ('normal', 'void'):
|
||||
return [self, cell]
|
||||
elif cell._type == 'fill':
|
||||
if cell.translation is not None:
|
||||
p -= cell.translation
|
||||
if cell.rotation is not None:
|
||||
p[:] = cell.rotation_matrix.dot(p)
|
||||
return [self, cell] + cell.fill.find(p)
|
||||
else:
|
||||
return [self, cell] + cell.fill.find(p)
|
||||
return []
|
||||
|
||||
def plot(self, center=(0., 0., 0.), width=(1., 1.), pixels=(200, 200),
|
||||
basis='xy', color_by='cell'):
|
||||
"""Display a slice plot of the universe.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
center : Iterable of float
|
||||
Coordinates at the center of the plot
|
||||
width : Iterable of float
|
||||
Width of the plot in each basis direction
|
||||
pixels : Iterable of int
|
||||
Number of pixels to use in each basis direction
|
||||
basis : {'xy', 'xz', 'yz'}
|
||||
The basis directions for the plot
|
||||
color_by : {'cell', 'material'}
|
||||
Indicate whether the plot should be colored by cell or by material
|
||||
|
||||
"""
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
if basis == 'xy':
|
||||
x_min = center[0] - 0.5*width[0]
|
||||
x_max = center[0] + 0.5*width[0]
|
||||
y_min = center[1] - 0.5*width[1]
|
||||
y_max = center[1] + 0.5*width[1]
|
||||
elif basis == 'yz':
|
||||
# The x-axis will correspond to physical y and the y-axis will correspond to physical z
|
||||
x_min = center[1] - 0.5*width[0]
|
||||
x_max = center[1] + 0.5*width[0]
|
||||
y_min = center[2] - 0.5*width[1]
|
||||
y_max = center[2] + 0.5*width[1]
|
||||
elif basis == 'xz':
|
||||
# The y-axis will correspond to physical z
|
||||
x_min = center[0] - 0.5*width[0]
|
||||
x_max = center[0] + 0.5*width[0]
|
||||
y_min = center[2] - 0.5*width[1]
|
||||
y_max = center[2] + 0.5*width[1]
|
||||
|
||||
# Determine locations to determine cells at
|
||||
x_coords = np.linspace(x_min, x_max, pixels[0], endpoint=False) + \
|
||||
0.5*(x_max - x_min)/pixels[0]
|
||||
y_coords = np.linspace(y_max, y_min, pixels[1], endpoint=False) - \
|
||||
0.5*(y_max - y_min)/pixels[1]
|
||||
|
||||
colors = {}
|
||||
img = np.zeros(pixels + (4,)) # Use RGBA form
|
||||
for i, x in enumerate(x_coords):
|
||||
for j, y in enumerate(y_coords):
|
||||
if basis == 'xy':
|
||||
path = self.find((x, y, center[2]))
|
||||
elif basis == 'yz':
|
||||
path = self.find((center[0], x, y))
|
||||
elif basis == 'xz':
|
||||
path = self.find((x, center[1], y))
|
||||
|
||||
if len(path) > 0:
|
||||
try:
|
||||
if color_by == 'cell':
|
||||
uid = path[-1].id
|
||||
elif color_by == 'material':
|
||||
if path[-1].fill_type == 'material':
|
||||
uid = path[-1].fill.id
|
||||
else:
|
||||
continue
|
||||
except AttributeError:
|
||||
continue
|
||||
if uid not in colors:
|
||||
colors[uid] = (random.random(), random.random(),
|
||||
random.random(), 1.0)
|
||||
img[j,i,:] = colors[uid]
|
||||
|
||||
plt.imshow(img, extent=(x_min, x_max, y_min, y_max))
|
||||
plt.show()
|
||||
|
||||
def add_cell(self, cell):
|
||||
"""Add a cell to the universe.
|
||||
|
||||
|
|
|
|||
|
|
@ -350,7 +350,6 @@ class InputSet(object):
|
|||
# Define fuel lattices.
|
||||
l100 = openmc.RectLattice(name='Fuel assembly (lower half)',
|
||||
lattice_id=100)
|
||||
l100.dimension = (17, 17)
|
||||
l100.lower_left = (-10.71, -10.71)
|
||||
l100.pitch = (1.26, 1.26)
|
||||
l100.universes = [
|
||||
|
|
@ -384,7 +383,6 @@ class InputSet(object):
|
|||
|
||||
l101 = openmc.RectLattice(name='Fuel assembly (upper half)',
|
||||
lattice_id=101)
|
||||
l101.dimension = (17, 17)
|
||||
l101.lower_left = (-10.71, -10.71)
|
||||
l101.pitch = (1.26, 1.26)
|
||||
l101.universes = [
|
||||
|
|
@ -444,7 +442,6 @@ class InputSet(object):
|
|||
# Define core lattices
|
||||
l200 = openmc.RectLattice(name='Core lattice (lower half)',
|
||||
lattice_id=200)
|
||||
l200.dimension = (21, 21)
|
||||
l200.lower_left = (-224.91, -224.91)
|
||||
l200.pitch = (21.42, 21.42)
|
||||
l200.universes = [
|
||||
|
|
@ -472,7 +469,6 @@ class InputSet(object):
|
|||
|
||||
l201 = openmc.RectLattice(name='Core lattice (lower half)',
|
||||
lattice_id=201)
|
||||
l201.dimension = (21, 21)
|
||||
l201.lower_left = (-224.91, -224.91)
|
||||
l201.pitch = (21.42, 21.42)
|
||||
l201.universes = [
|
||||
|
|
|
|||
|
|
@ -24,7 +24,6 @@ class AsymmetricLatticeTestHarness(PyAPITestHarness):
|
|||
|
||||
# Construct a 3x3 lattice of fuel assemblies
|
||||
core_lat = openmc.RectLattice(name='3x3 Core Lattice', lattice_id=202)
|
||||
core_lat.dimension = (3, 3)
|
||||
core_lat.lower_left = (-32.13, -32.13)
|
||||
core_lat.pitch = (21.42, 21.42)
|
||||
core_lat.universes = [[fuel, water, water],
|
||||
|
|
|
|||
|
|
@ -53,7 +53,6 @@ class DistribmatTestHarness(PyAPITestHarness):
|
|||
fuel_univ.add_cells((c11, c12))
|
||||
|
||||
lat = openmc.RectLattice(lattice_id=101)
|
||||
lat.dimension = [2, 2]
|
||||
lat.lower_left = [-2.0, -2.0]
|
||||
lat.pitch = [2.0, 2.0]
|
||||
lat.universes = [[fuel_univ]*2]*2
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue