diff --git a/openmc/model/__init__.py b/openmc/model/__init__.py index 9fa999dd4e..8ec96845ec 100644 --- a/openmc/model/__init__.py +++ b/openmc/model/__init__.py @@ -1,3 +1,4 @@ from .triso import * from .model import * from .funcs import * +from .pin import * diff --git a/openmc/model/pin.py b/openmc/model/pin.py new file mode 100644 index 0000000000..0ec0c4dbdc --- /dev/null +++ b/openmc/model/pin.py @@ -0,0 +1,249 @@ +""" +Helper class for building a pin from concentric cylinders +""" + + +from math import sqrt +from numbers import Real +from operator import attrgetter + +import openmc +import openmc.checkvalue as cv +from . import subdivide + +__all__ = ["Pin"] + + +class Pin(openmc.Universe): + """Special universe used to model pins + + Parameters + ---------- + surfaces: iterable of :class:`openmc.Cylinder` + Cylinders used to define boundaries + between materials. All cylinders must be + concentric and of the same orientation, e.g. + all :class:`openmc.ZCylinders` + materials: iterable of :class:`openmc.Material` + Materials to go between ``surfaces``. There must be one + more material than surfaces, corresponding to the material + that spans all space outside the final ring. + universe_id: None or int + Unique identifier for this universe + name: str + Name for this universe + + See Also + -------- + + :meth:`openmc.Pin.from_radii` - Convinience function to build + from radii of cylinders + """ + + def __init__(self, surfaces, materials, universe_id=None, name=""): + cv.check_iterable_type("materials", materials, openmc.Material) + cv.check_length("surfaces", surfaces, len(materials) - 1) + + # Ensure that all surfaces are same type of cylinder + self._check_surfaces(surfaces) + regions = subdivide(surfaces) + + cells = [ + openmc.Cell(fill=m, region=r) for m, r in zip(materials, regions) + ] + + super().__init__(universe_id=universe_id, name=name, cells=cells) + self._list_cells = cells # need ordered by radial position + + def _check_surfaces(self, surfaces): + cv.check_type( + "surface 0", + surfaces[0], + (openmc.ZCylinder, openmc.YCylinder, openmc.XCylinder), + ) + if isinstance(surfaces[0], openmc.ZCylinder): + center_getter = attrgetter("x0", "y0") + elif isinstance(surfaces[0], openmc.YClylinder): + center_getter = attrgetter("x0", "y0") + elif isinstance(surfaces[0], openmc.XClylinder): + center_getter = attrgetter("z0", "y0") + else: + raise TypeError( + "Not configured to interpret {} surfaces".format( + surfaces[0].__class__.__name__ + ) + ) + cv.check_iterable_type("surfaces", surfaces[1:], type(surfaces[0])) + # Check for concentric-ness and increasing radii + centers = set() + radii = [] + rad = 0 + for ix, surf in enumerate(surfaces): + cur_rad = surf.r + if cur_rad <= rad: + raise ValueError( + "Surfaces do not appear to be increasing in radius. " + "Surface {} at index {} has radius {:7.3E} compared to " + "previous radius of {:7.5E}".format( + surf.id, ix, cur_rad, rad + ) + ) + rad = cur_rad + radii.append(cur_rad) + centers.add(center_getter(surf)) + + if len(centers) > 1: + raise ValueError( + "Surfaces do not appear to be concentric. The following " + "centers were found: {}".format(centers) + ) + self._radii = radii + self._surfaces = surfaces + self._surf_type = type(surfaces[0]) + + @classmethod + def from_radii( + cls, + radii, + materials, + universe_id=None, + name="", + orientation="z", + center=(0.0, 0.0), + ): + """Construct using radii of concentric cylinders and materials + + Parameters + ---------- + radii: iterable of float + Radii of the intended cylinders. Must be all positive + values and increasing + materials: iterable of :class:`openmc.Material` + Materials used to fill cylinders created by ``radii``, + starting from inside to the outside of the pin. + There must be one extra material corresponding to all + area outside the last ring + universe_id: None or int + Unique identifier to give this pin + name: str + Name of the created universe + orientation: {"x", "y", "z"} + Axis along which to orient the pin. Default is ``"z"`` + center: iterable of float + Center of the pin in the plane perpendicular + """ + cv.check_iterable_type("materials", materials, openmc.Material) + cv.check_length("radii", radii, len(materials) - 1) + for ix, rad in enumerate(radii): + if rad < 0: + raise ValueError( + "Radius {:7.3E} at index {} is non-positive".format( + rad, ix + ) + ) + if ix and rad <= radii[ix - 1]: + raise ValueError( + "Radii must be increasing values. Radius {:7.3E} at index " + "{} is not greater than previous value of {:7.3E}".format( + rad, ix, radii[ix - 1] + ) + ) + + if orientation == "z": + surfCls = openmc.ZCylinder + basis = ["x0", "y0"] + elif orientation == "y": + surfCls = openmc.YCylinder + basis = ["x0", "z0"] + elif orientation == "z": + surfCls = openmc.XCylinder + basis = ["y0", "z0"] + else: + raise ValueError( + "Orientation of {} not understood".format(orientation) + ) + + centerKwargs = dict(zip(basis, center)) + + surfaces = [surfCls(r=rad, **centerKwargs) for rad in radii] + + return cls(surfaces, materials, universe_id=universe_id, name=name) + + def subdivide_ring(self, ring_index, n_divs): + """Divide one ring of the pin into equal-area rings + + Each new ring will be added to the model, and filled with + a unique material copied from the original. + + Parameters + ring_index: int + Index of the ring to be divided where 0 is the innermost + ring. Will not divide the outermost region as there is no + upper bound + n_divs: int + Number of equal area divisions to make in this ring + """ + # Don't allow subdivision of outer, infinite region + cv.check_less_than("ring_index", ring_index, len(self._radii)) + cv.check_type("n_divs", n_divs, Real) + cv.check_greater_than("n_divs", n_divs, 1) + + if ring_index < 0: + ring_index = len(self._list_cells) + ring_index + + # Get all the information we need to replicate this + # region with unique insides + orig_cell = self._list_cells[ring_index] + + lower_rad = self._radii[ring_index - 1] if ring_index else 0.0 + area_term = (self._radii[ring_index] ** 2 - lower_rad ** 2) / n_divs + + new_radii = [] + new_surfaces = [] + new_cells = [] + + # Adding N - 1 new regions + # N - 2 surfaces are made + # Original cell is not removed, but not occupies last ring + + for i in range(n_divs - 1): + r = sqrt(area_term + lower_rad ** 2) + lower_rad = r + new_radii.append(r) + surf = self._surf_type(r=r) + new_surfaces.append(surf) + if i == 0: + if ring_index: + region = ( + -surf & +self._surfaces[ring_index - 1] + ) + else: + region = -surf + else: + region = -surf & +new_surfaces[-2] + new_cells.append( + openmc.Cell(region=region, fill=orig_cell.fill.clone()) + ) + + orig_cell.region = -self._surfaces[ring_index] & +surf + + self.add_cells(new_cells) + + self._list_cells = ( + self._list_cells[:ring_index] + + new_cells + + self._list_cells[ring_index:] + ) + self._radii = ( + self._radii[:ring_index] + new_radii + self._radii[ring_index:] + ) + self._surfaces = ( + self._surfaces[:ring_index] + + new_surfaces + + self._surfaces[ring_index:] + ) + + @property + def radii(self): + """Return a tuple of the radii in this :class:`Pin`""" + return tuple(self._radii) diff --git a/tests/unit_tests/test_pin.py b/tests/unit_tests/test_pin.py new file mode 100644 index 0000000000..132e5f026a --- /dev/null +++ b/tests/unit_tests/test_pin.py @@ -0,0 +1,86 @@ +""" +Tests for constructing Pin universes +""" + +import numpy +import pytest + +import openmc +from openmc.model import Pin + + +@pytest.fixture +def pin_mats(): + fuel = openmc.Material(name="UO2") + clad = openmc.Material(name="zirc") + water = openmc.Material(name="water") + return fuel, clad, water + + +@pytest.fixture +def good_radii(): + return (0.4, 0.42) + + +def test_failure(pin_mats, good_radii): + """Check for various failure modes""" + # Bad material type + with pytest.raises(TypeError): + Pin.from_radii(good_radii, [mat.name for mat in pin_mats]) + + # Incorrect lengths + with pytest.raises(ValueError) as exec_info: + Pin.from_radii(good_radii[: len(pin_mats) - 2], pin_mats) + assert "length" in str(exec_info) + + # Non-positive radii + rad = (-0.1,) + good_radii[1:] + with pytest.raises(ValueError) as exec_info: + Pin.from_radii(rad, pin_mats) + assert "index 0" in str(exec_info) + + # Non-increasing radii + rad = tuple(reversed(good_radii)) + with pytest.raises(ValueError) as exec_info: + Pin.from_radii(rad, pin_mats) + assert "index 1" in str(exec_info) + + # Bad orientation + with pytest.raises(ValueError) as exec_info: + Pin.from_radii(good_radii, pin_mats, orientation="fail") + assert "Orientation" in str(exec_info) + + +def test_from_radii(pin_mats, good_radii): + name = "test pin" + p = Pin.from_radii(good_radii, pin_mats, name=name) + assert len(p.cells) == len(pin_mats) + assert p.name == name + assert p.radii == good_radii + +def test_subdivide(pin_mats, good_radii): + surfs = [openmc.ZCylinder(r=r) for r in good_radii] + pin = Pin(surfs, pin_mats) + assert pin.radii == good_radii + assert len(pin.cells) == len(pin_mats) + + # subdivide inner region + N = 5 + pin.subdivide_ring(0, N) + assert len(pin.radii) == len(good_radii) + N - 1 + assert len(pin.cells) == len(pin_mats) + N - 1 + # check volumes of new rings + bounds = (0,) + pin.radii[:N] + sqrs = numpy.square(bounds) + assert sqrs[1:] - sqrs[:-1] == pytest.approx(good_radii[0] ** 2 / N) + + # subdivide non-inner most region + new_pin = Pin.from_radii(good_radii, pin_mats) + new_pin.subdivide_ring(1, N) + assert len(new_pin.radii) == len(good_radii) + N - 1 + assert len(new_pin.cells) == len(pin_mats) + N - 1 + # check volumes of new rings + bounds = new_pin.radii[:N + 1] + sqrs = numpy.square(bounds) + assert sqrs[1:] - sqrs[:-1] == pytest.approx((good_radii[1] ** 2 - good_radii[0] ** 2) / N) +