mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-25 12:35:29 -04:00
560 lines
17 KiB
Python
560 lines
17 KiB
Python
from abc import ABCMeta
|
|
from collections.abc import Iterable
|
|
from numbers import Real, Integral
|
|
from xml.etree import ElementTree as ET
|
|
import sys
|
|
import warnings
|
|
|
|
import numpy as np
|
|
|
|
import openmc.checkvalue as cv
|
|
import openmc
|
|
from openmc._xml import get_text
|
|
from openmc.mixin import EqualityMixin, IDManagerMixin
|
|
|
|
|
|
class MeshBase(IDManagerMixin, metaclass=ABCMeta):
|
|
"""A mesh that partitions geometry for tallying purposes.
|
|
|
|
Parameters
|
|
----------
|
|
mesh_id : int
|
|
Unique identifier for the mesh
|
|
name : str
|
|
Name of the mesh
|
|
|
|
Attributes
|
|
----------
|
|
id : int
|
|
Unique identifier for the mesh
|
|
name : str
|
|
Name of the mesh
|
|
|
|
"""
|
|
|
|
next_id = 1
|
|
used_ids = set()
|
|
|
|
def __init__(self, mesh_id=None, name=''):
|
|
# Initialize Mesh class attributes
|
|
self.id = mesh_id
|
|
self.name = name
|
|
|
|
@property
|
|
def name(self):
|
|
return self._name
|
|
|
|
@name.setter
|
|
def name(self, name):
|
|
if name is not None:
|
|
cv.check_type('name for mesh ID="{0}"'.format(self._id),
|
|
name, str)
|
|
self._name = name
|
|
else:
|
|
self._name = ''
|
|
|
|
@classmethod
|
|
def from_hdf5(cls, group):
|
|
"""Create mesh from HDF5 group
|
|
|
|
Parameters
|
|
----------
|
|
group : h5py.Group
|
|
Group in HDF5 file
|
|
|
|
Returns
|
|
-------
|
|
openmc.MeshBase
|
|
Instance of a MeshBase subclass
|
|
|
|
"""
|
|
|
|
mesh_type = group['type'][()].decode()
|
|
if mesh_type == 'regular':
|
|
return RegularMesh.from_hdf5(group)
|
|
elif mesh_type == 'rectilinear':
|
|
return RectilinearMesh.from_hdf5(group)
|
|
else:
|
|
raise ValueError('Unrecognized mesh type: "' + mesh_type + '"')
|
|
|
|
|
|
class RegularMesh(MeshBase):
|
|
"""A regular Cartesian mesh in one, two, or three dimensions
|
|
|
|
Parameters
|
|
----------
|
|
mesh_id : int
|
|
Unique identifier for the mesh
|
|
name : str
|
|
Name of the mesh
|
|
|
|
Attributes
|
|
----------
|
|
id : int
|
|
Unique identifier for the mesh
|
|
name : str
|
|
Name of the mesh
|
|
dimension : Iterable of int
|
|
The number of mesh cells in each direction.
|
|
n_dimension : int
|
|
Number of mesh dimensions.
|
|
lower_left : Iterable of float
|
|
The lower-left corner of the structured mesh. If only two coordinate are
|
|
given, it is assumed that the mesh is an x-y mesh.
|
|
upper_right : Iterable of float
|
|
The upper-right corner of the structrued mesh. If only two coordinate
|
|
are given, it is assumed that the mesh is an x-y mesh.
|
|
width : Iterable of float
|
|
The width of mesh cells in each direction.
|
|
indices : Iterable of tuple
|
|
An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1),
|
|
(2, 1, 1), ...]
|
|
|
|
"""
|
|
|
|
def __init__(self, mesh_id=None, name=''):
|
|
super().__init__(mesh_id, name)
|
|
|
|
self._dimension = None
|
|
self._lower_left = None
|
|
self._upper_right = None
|
|
self._width = None
|
|
|
|
@property
|
|
def dimension(self):
|
|
return self._dimension
|
|
|
|
@property
|
|
def n_dimension(self):
|
|
return len(self._dimension)
|
|
|
|
@property
|
|
def lower_left(self):
|
|
return self._lower_left
|
|
|
|
@property
|
|
def upper_right(self):
|
|
return self._upper_right
|
|
|
|
@property
|
|
def width(self):
|
|
return self._width
|
|
|
|
@property
|
|
def num_mesh_cells(self):
|
|
return np.prod(self._dimension)
|
|
|
|
@property
|
|
def indices(self):
|
|
ndim = len(self._dimension)
|
|
if ndim == 3:
|
|
nx, ny, nz = self.dimension
|
|
return ((x, y, z)
|
|
for z in range(1, nz + 1)
|
|
for y in range(1, ny + 1)
|
|
for x in range(1, nx + 1))
|
|
elif ndim == 2:
|
|
nx, ny = self.dimension
|
|
return ((x, y)
|
|
for y in range(1, ny + 1)
|
|
for x in range(1, nx + 1))
|
|
else:
|
|
nx, = self.dimension
|
|
return ((x,) for x in range(1, nx + 1))
|
|
|
|
@dimension.setter
|
|
def dimension(self, dimension):
|
|
cv.check_type('mesh dimension', dimension, Iterable, Integral)
|
|
cv.check_length('mesh dimension', dimension, 1, 3)
|
|
self._dimension = dimension
|
|
|
|
@lower_left.setter
|
|
def lower_left(self, lower_left):
|
|
cv.check_type('mesh lower_left', lower_left, Iterable, Real)
|
|
cv.check_length('mesh lower_left', lower_left, 1, 3)
|
|
self._lower_left = lower_left
|
|
|
|
@upper_right.setter
|
|
def upper_right(self, upper_right):
|
|
cv.check_type('mesh upper_right', upper_right, Iterable, Real)
|
|
cv.check_length('mesh upper_right', upper_right, 1, 3)
|
|
self._upper_right = upper_right
|
|
|
|
@width.setter
|
|
def width(self, width):
|
|
cv.check_type('mesh width', width, Iterable, Real)
|
|
cv.check_length('mesh width', width, 1, 3)
|
|
self._width = width
|
|
|
|
def __repr__(self):
|
|
string = 'RegularMesh\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('\tType', '=\t', self._type)
|
|
string += '{0: <16}{1}{2}\n'.format('\tBasis', '=\t', self._dimension)
|
|
string += '{0: <16}{1}{2}\n'.format('\tWidth', '=\t', self._lower_left)
|
|
string += '{0: <16}{1}{2}\n'.format('\tOrigin', '=\t', self._upper_right)
|
|
string += '{0: <16}{1}{2}\n'.format('\tPixels', '=\t', self._width)
|
|
return string
|
|
|
|
@classmethod
|
|
def from_hdf5(cls, group):
|
|
mesh_id = int(group.name.split('/')[-1].lstrip('mesh '))
|
|
|
|
# Read and assign mesh properties
|
|
mesh = cls(mesh_id)
|
|
mesh.dimension = group['dimension'][()]
|
|
mesh.lower_left = group['lower_left'][()]
|
|
mesh.upper_right = group['upper_right'][()]
|
|
mesh.width = group['width'][()]
|
|
|
|
return mesh
|
|
|
|
@classmethod
|
|
def from_rect_lattice(cls, lattice, division=1, mesh_id=None, name=''):
|
|
"""Create mesh from an existing rectangular lattice
|
|
|
|
Parameters
|
|
----------
|
|
lattice : openmc.RectLattice
|
|
Rectangular lattice used as a template for this mesh
|
|
division : int
|
|
Number of mesh cells per lattice cell.
|
|
If not specified, there will be 1 mesh cell per lattice cell.
|
|
mesh_id : int
|
|
Unique identifier for the mesh
|
|
name : str
|
|
Name of the mesh
|
|
|
|
Returns
|
|
-------
|
|
openmc.RegularMesh
|
|
RegularMesh instance
|
|
|
|
"""
|
|
cv.check_type('rectangular lattice', lattice, openmc.RectLattice)
|
|
|
|
shape = np.array(lattice.shape)
|
|
width = lattice.pitch*shape
|
|
|
|
mesh = cls(mesh_id, name)
|
|
mesh.lower_left = lattice.lower_left
|
|
mesh.upper_right = lattice.lower_left + width
|
|
mesh.dimension = shape*division
|
|
|
|
return mesh
|
|
|
|
def to_xml_element(self):
|
|
"""Return XML representation of the mesh
|
|
|
|
Returns
|
|
-------
|
|
element : xml.etree.ElementTree.Element
|
|
XML element containing mesh data
|
|
|
|
"""
|
|
|
|
element = ET.Element("mesh")
|
|
element.set("id", str(self._id))
|
|
|
|
if self._dimension is not None:
|
|
subelement = ET.SubElement(element, "dimension")
|
|
subelement.text = ' '.join(map(str, self._dimension))
|
|
|
|
subelement = ET.SubElement(element, "lower_left")
|
|
subelement.text = ' '.join(map(str, self._lower_left))
|
|
|
|
if self._upper_right is not None:
|
|
subelement = ET.SubElement(element, "upper_right")
|
|
subelement.text = ' '.join(map(str, self._upper_right))
|
|
|
|
if self._width is not None:
|
|
subelement = ET.SubElement(element, "width")
|
|
subelement.text = ' '.join(map(str, self._width))
|
|
|
|
return element
|
|
|
|
@classmethod
|
|
def from_xml_element(cls, elem):
|
|
"""Generate mesh from an XML element
|
|
|
|
Parameters
|
|
----------
|
|
elem : xml.etree.ElementTree.Element
|
|
XML element
|
|
|
|
Returns
|
|
-------
|
|
openmc.Mesh
|
|
Mesh generated from XML element
|
|
|
|
"""
|
|
mesh_id = int(get_text(elem, 'id'))
|
|
mesh = cls(mesh_id)
|
|
|
|
mesh_type = get_text(elem, 'type')
|
|
if mesh_type is not None:
|
|
mesh.type = mesh_type
|
|
|
|
dimension = get_text(elem, 'dimension')
|
|
if dimension is not None:
|
|
mesh.dimension = [int(x) for x in dimension.split()]
|
|
|
|
lower_left = get_text(elem, 'lower_left')
|
|
if lower_left is not None:
|
|
mesh.lower_left = [float(x) for x in lower_left.split()]
|
|
|
|
upper_right = get_text(elem, 'upper_right')
|
|
if upper_right is not None:
|
|
mesh.upper_right = [float(x) for x in upper_right.split()]
|
|
|
|
width = get_text(elem, 'width')
|
|
if width is not None:
|
|
mesh.width = [float(x) for x in width.split()]
|
|
|
|
return mesh
|
|
|
|
def build_cells(self, bc=['reflective'] * 6):
|
|
"""Generates a lattice of universes with the same dimensionality
|
|
as the mesh object. The individual cells/universes produced
|
|
will not have material definitions applied and so downstream code
|
|
will have to apply that information.
|
|
|
|
Parameters
|
|
----------
|
|
bc : iterable of {'reflective', 'periodic', 'transmission', or 'vacuum'}
|
|
Boundary conditions for each of the four faces of a rectangle
|
|
(if aplying to a 2D mesh) or six faces of a parallelepiped
|
|
(if applying to a 3D mesh) provided in the following order:
|
|
[x min, x max, y min, y max, z min, z max]. 2-D cells do not
|
|
contain the z min and z max entries.
|
|
|
|
Returns
|
|
-------
|
|
root_cell : openmc.Cell
|
|
The cell containing the lattice representing the mesh geometry;
|
|
this cell is a single parallelepiped with boundaries matching
|
|
the outermost mesh boundary with the boundary conditions from bc
|
|
applied.
|
|
cells : iterable of openmc.Cell
|
|
The list of cells within each lattice position mimicking the mesh
|
|
geometry.
|
|
|
|
"""
|
|
|
|
cv.check_length('bc', bc, length_min=4, length_max=6)
|
|
for entry in bc:
|
|
cv.check_value('bc', entry, ['transmission', 'vacuum',
|
|
'reflective', 'periodic'])
|
|
|
|
n_dim = len(self.dimension)
|
|
|
|
# Build the cell which will contain the lattice
|
|
xplanes = [openmc.XPlane(self.lower_left[0], bc[0]),
|
|
openmc.XPlane(self.upper_right[0], bc[1])]
|
|
if n_dim == 1:
|
|
yplanes = [openmc.YPlane(-1e10, 'reflective'),
|
|
openmc.YPlane(1e10, 'reflective')]
|
|
else:
|
|
yplanes = [openmc.YPlane(self.lower_left[1], bc[2]),
|
|
openmc.YPlane(self.upper_right[1], bc[3])]
|
|
|
|
if n_dim <= 2:
|
|
# Would prefer to have the z ranges be the max supported float, but
|
|
# these values are apparently different between python and Fortran.
|
|
# Choosing a safe and sane default.
|
|
# Values of +/-1e10 are used here as there seems to be an
|
|
# inconsistency between what numpy uses as the max float and what
|
|
# Fortran expects for a real(8), so this avoids code complication
|
|
# and achieves the same goal.
|
|
zplanes = [openmc.ZPlane(-1e10, 'reflective'),
|
|
openmc.ZPlane(1e10, 'reflective')]
|
|
else:
|
|
zplanes = [openmc.ZPlane(self.lower_left[2], bc[4]),
|
|
openmc.ZPlane(self.upper_right[2], bc[5])]
|
|
root_cell = openmc.Cell()
|
|
root_cell.region = ((+xplanes[0] & -xplanes[1]) &
|
|
(+yplanes[0] & -yplanes[1]) &
|
|
(+zplanes[0] & -zplanes[1]))
|
|
|
|
# Build the universes which will be used for each of the (i,j,k)
|
|
# locations within the mesh.
|
|
# We will concurrently build cells to assign to these universes
|
|
cells = []
|
|
universes = []
|
|
for index in self.indices:
|
|
cells.append(openmc.Cell())
|
|
universes.append(openmc.Universe())
|
|
universes[-1].add_cell(cells[-1])
|
|
|
|
lattice = openmc.RectLattice()
|
|
lattice.lower_left = self.lower_left
|
|
|
|
# Assign the universe and rotate to match the indexing expected for
|
|
# the lattice
|
|
if n_dim == 1:
|
|
universe_array = np.array([universes])
|
|
elif n_dim == 2:
|
|
universe_array = np.empty(self.dimension[::-1],
|
|
dtype=openmc.Universe)
|
|
i = 0
|
|
for y in range(self.dimension[1] - 1, -1, -1):
|
|
for x in range(self.dimension[0]):
|
|
universe_array[y][x] = universes[i]
|
|
i += 1
|
|
else:
|
|
universe_array = np.empty(self.dimension[::-1],
|
|
dtype=openmc.Universe)
|
|
i = 0
|
|
for z in range(self.dimension[2]):
|
|
for y in range(self.dimension[1] - 1, -1, -1):
|
|
for x in range(self.dimension[0]):
|
|
universe_array[z][y][x] = universes[i]
|
|
i += 1
|
|
lattice.universes = universe_array
|
|
|
|
if self.width is not None:
|
|
lattice.pitch = self.width
|
|
else:
|
|
dx = ((self.upper_right[0] - self.lower_left[0]) /
|
|
self.dimension[0])
|
|
|
|
if n_dim == 1:
|
|
lattice.pitch = [dx]
|
|
elif n_dim == 2:
|
|
dy = ((self.upper_right[1] - self.lower_left[1]) /
|
|
self.dimension[1])
|
|
lattice.pitch = [dx, dy]
|
|
else:
|
|
dy = ((self.upper_right[1] - self.lower_left[1]) /
|
|
self.dimension[1])
|
|
dz = ((self.upper_right[2] - self.lower_left[2]) /
|
|
self.dimension[2])
|
|
lattice.pitch = [dx, dy, dz]
|
|
|
|
# Fill Cell with the Lattice
|
|
root_cell.fill = lattice
|
|
|
|
return root_cell, cells
|
|
|
|
|
|
def Mesh(*args, **kwargs):
|
|
warnings.warn("Mesh has been renamed RegularMesh. Future versions of "
|
|
"OpenMC will not accept the name Mesh.")
|
|
return RegularMesh(*args, **kwargs)
|
|
|
|
|
|
class RectilinearMesh(MeshBase):
|
|
"""A 3D rectilinear Cartesian mesh
|
|
|
|
Parameters
|
|
----------
|
|
mesh_id : int
|
|
Unique identifier for the mesh
|
|
name : str
|
|
Name of the mesh
|
|
|
|
Attributes
|
|
----------
|
|
id : int
|
|
Unique identifier for the mesh
|
|
name : str
|
|
Name of the mesh
|
|
n_dimension : int
|
|
Number of mesh dimensions (always 3 for a RectilinearMesh).
|
|
x_grid : Iterable of float
|
|
Mesh boundary points along the x-axis.
|
|
y_grid : Iterable of float
|
|
Mesh boundary points along the y-axis.
|
|
z_grid : Iterable of float
|
|
Mesh boundary points along the z-axis.
|
|
indices : Iterable of tuple
|
|
An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1),
|
|
(2, 1, 1), ...]
|
|
|
|
"""
|
|
|
|
def __init__(self, mesh_id=None, name=''):
|
|
super().__init__(mesh_id, name)
|
|
|
|
self._x_grid = None
|
|
self._y_grid = None
|
|
self._z_grid = None
|
|
|
|
@property
|
|
def n_dimension(self):
|
|
return 3
|
|
|
|
@property
|
|
def x_grid(self):
|
|
return self._x_grid
|
|
|
|
@property
|
|
def y_grid(self):
|
|
return self._y_grid
|
|
|
|
@property
|
|
def z_grid(self):
|
|
return self._z_grid
|
|
|
|
@property
|
|
def indices(self):
|
|
nx = len(self.x_grid) - 1
|
|
ny = len(self.y_grid) - 1
|
|
nz = len(self.z_grid) - 1
|
|
return ((x, y, z)
|
|
for z in range(1, nz + 1)
|
|
for y in range(1, ny + 1)
|
|
for x in range(1, nx + 1))
|
|
|
|
@x_grid.setter
|
|
def x_grid(self, grid):
|
|
cv.check_type('mesh x_grid', grid, Iterable, Real)
|
|
self._x_grid = grid
|
|
|
|
@y_grid.setter
|
|
def y_grid(self, grid):
|
|
cv.check_type('mesh y_grid', grid, Iterable, Real)
|
|
self._y_grid = grid
|
|
|
|
@z_grid.setter
|
|
def z_grid(self, grid):
|
|
cv.check_type('mesh z_grid', grid, Iterable, Real)
|
|
self._z_grid = grid
|
|
|
|
@classmethod
|
|
def from_hdf5(cls, group):
|
|
mesh_id = int(group.name.split('/')[-1].lstrip('mesh '))
|
|
|
|
# Read and assign mesh properties
|
|
mesh = cls(mesh_id)
|
|
mesh.x_grid = group['x_grid'][()]
|
|
mesh.y_grid = group['y_grid'][()]
|
|
mesh.z_grid = group['z_grid'][()]
|
|
|
|
return mesh
|
|
|
|
def to_xml_element(self):
|
|
"""Return XML representation of the mesh
|
|
|
|
Returns
|
|
-------
|
|
element : xml.etree.ElementTree.Element
|
|
XML element containing mesh data
|
|
|
|
"""
|
|
|
|
element = ET.Element("mesh")
|
|
element.set("id", str(self._id))
|
|
element.set("type", "rectilinear")
|
|
|
|
subelement = ET.SubElement(element, "x_grid")
|
|
subelement.text = ' '.join(map(str, self.x_grid))
|
|
|
|
subelement = ET.SubElement(element, "y_grid")
|
|
subelement.text = ' '.join(map(str, self.y_grid))
|
|
|
|
subelement = ET.SubElement(element, "z_grid")
|
|
subelement.text = ' '.join(map(str, self.z_grid))
|
|
|
|
return element
|