From 1475e25bf344ac397358c872c1d2e5d19d98dccc Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sun, 30 Jun 2019 13:26:43 -0500 Subject: [PATCH] Add Universe subclass openmc.model.Pin Designed to facilitate building universe that represent pins, or concentric cylinders of materials. Can be built by passing surfaces and materials, or radii and materials. Places last material spanning region out to infinity from the last surface. Supports subdividing rings into equal volume slices of unique [cloned] materials, useful for depletion. Unit tests added for failure modes as well as basic operation --- openmc/model/__init__.py | 1 + openmc/model/pin.py | 249 +++++++++++++++++++++++++++++++++++ tests/unit_tests/test_pin.py | 86 ++++++++++++ 3 files changed, 336 insertions(+) create mode 100644 openmc/model/pin.py create mode 100644 tests/unit_tests/test_pin.py 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) +