From 50ea1ed02af4a89b4a4b4ff3c4c5ccc922a16f36 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 5 Jul 2019 16:51:10 -0500 Subject: [PATCH] Function based pin construction Function openmc.model.funcs.pin has replaced the class-based openmc.model.pin.Pin Subdivision functionality is maintained by passing a dictionary of integer ring indexes -> number of divisions to the function Tests have been updated accordingly and are parametrized against surface type: X, Y, and Z cylinders --- openmc/model/__init__.py | 1 - openmc/model/funcs.py | 130 +++++++++++++++++- openmc/model/pin.py | 249 ----------------------------------- tests/unit_tests/test_pin.py | 75 ++++++----- 4 files changed, 167 insertions(+), 288 deletions(-) delete mode 100644 openmc/model/pin.py diff --git a/openmc/model/__init__.py b/openmc/model/__init__.py index 8ec96845ec..9fa999dd4e 100644 --- a/openmc/model/__init__.py +++ b/openmc/model/__init__.py @@ -1,4 +1,3 @@ from .triso import * from .model import * from .funcs import * -from .pin import * diff --git a/openmc/model/funcs.py b/openmc/model/funcs.py index 143033eca8..4b44b450ff 100644 --- a/openmc/model/funcs.py +++ b/openmc/model/funcs.py @@ -1,12 +1,16 @@ -from collections import OrderedDict from collections.abc import Iterable from math import sqrt from numbers import Real from functools import partial from warnings import warn +from operator import attrgetter -from openmc import XPlane, YPlane, Plane, ZCylinder, Quadric -from openmc.checkvalue import check_type, check_value +from openmc import ( + XPlane, YPlane, Plane, ZCylinder, Quadric, Cylinder, XCylinder, + YCylinder, Material, Universe, Cell) +from openmc.checkvalue import ( + check_type, check_value, check_length, check_less_than, + check_iterable_type) import openmc.data @@ -414,7 +418,7 @@ def cylinder_from_points(p1, p2, r, **kwargs): kwargs['j'] = cx*dy - cy*dx kwargs['k'] = -(dx*dx + dy*dy + dz*dz)*r*r - return openmc.Quadric(**kwargs) + return Quadric(**kwargs) def subdivide(surfaces): @@ -442,3 +446,121 @@ def subdivide(surfaces): regions.append(+s0 & -s1) regions.append(+surfaces[-1]) return regions + + +def pin(surfaces, materials, subdivisions=None, universe_id=None, name=""): + """Convenience function for building a fuel pin + + 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.ZCylinder` + 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. + subdivisions : None or dict of int to int + Dictionary describing which rings to subdivide and how + many times. Keys are indexes of the annular rings + to be divided. Will construct equal area rings + universe_id : None or int + Identifier for this universe + name : str + Name for this universe + + Returns + ------- + :class:`openmc.Universe` + Universe of concentric cylinders filled with the desired + materials + """ + check_type("materials", materials, Iterable, Material) + check_length("surfaces", surfaces, len(materials) - 1, len(materials) - 1) + # Check that all surfaces are of similar orientation + check_type("surface", surfaces[0], Cylinder) + surf_type = type(surfaces[0]) + check_iterable_type("surfaces", surfaces[1:], surf_type) + + # Check for increasing radii and equal centers + if surf_type is ZCylinder: + center_getter = attrgetter("x0", "y0") + elif surf_type is YCylinder: + center_getter = attrgetter("x0", "z0") + elif surf_type is XCylinder: + center_getter = attrgetter("z0", "y0") + else: + raise TypeError( + "Not configured to interpret {} surfaces".format( + surf_type.__name__)) + + centers = set() + prev_rad = 0 + for ix, surf in enumerate(surfaces): + cur_rad = surf.r + if cur_rad <= prev_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, prev_rad)) + prev_rad = 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)) + + if subdivisions is not None: + check_length("subdivisions", subdivisions, 1, len(surfaces)) + orig_indexes = list(subdivisions.keys()) + check_iterable_type("ring indexes", orig_indexes, int) + check_iterable_type( + "number of divisions", list(subdivisions.values()), int) + for ix in orig_indexes: + if ix < 0: + subdivisions[len(surfaces) + ix] = subdivisions.pop(ix) + # Dissallow subdivision on outer most, infinite region + check_less_than( + "outer ring", max(subdivisions), len(surfaces), equality=True) + + # ensure ability to concatenate + if not isinstance(materials, list): + materials = list(materials) + if not isinstance(surfaces, list): + surfaces = list(surfaces) + + # generate equal area divisions + # Adding N - 1 new regions + # N - 2 surfaces are made + # Original cell is not removed, but not occupies last ring + for ring_index in reversed(sorted(subdivisions.keys())): + nr = subdivisions[ring_index] + new_surfs = [] + + if ring_index == 0: + lower_rad = 0.0 + else: + lower_rad = surfaces[ring_index - 1].r + upper_rad = surfaces[ring_index].r + + area_term = (upper_rad ** 2 - lower_rad ** 2) / nr + + for new_index in range(nr - 1): + lower_rad = sqrt(area_term + lower_rad ** 2) + new_surfs.append(surf_type(r=lower_rad)) + + surfaces = ( + surfaces[:ring_index] + new_surfs + surfaces[ring_index:]) + materials = ( + materials[:ring_index] + + [materials[ring_index].clone() for _i in range(nr - 1)] + + materials[ring_index:]) + + # Build the universe + regions = subdivide(surfaces) + cells = [Cell(fill=f, region=r) for r, f in zip(regions, materials)] + return Universe(universe_id=universe_id, name=name, cells=cells) diff --git a/openmc/model/pin.py b/openmc/model/pin.py deleted file mode 100644 index 0ec0c4dbdc..0000000000 --- a/openmc/model/pin.py +++ /dev/null @@ -1,249 +0,0 @@ -""" -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 index 0940c5e30b..d2f7fa8016 100644 --- a/tests/unit_tests/test_pin.py +++ b/tests/unit_tests/test_pin.py @@ -6,7 +6,18 @@ import numpy import pytest import openmc -from openmc.model import Pin +from openmc.model import pin + + +def get_pin_radii(pin_univ): + """Return a sorted list of all radii from pin""" + rads = set() + + for cell in pin_univ.get_all_cells().values(): + surfs = cell.region.get_surfaces().values() + rads.update(set(s.r for s in surfs)) + + return list(sorted(rads)) @pytest.fixture @@ -24,60 +35,56 @@ def good_radii(): def test_failure(pin_mats, good_radii): """Check for various failure modes""" + good_surfaces = [openmc.ZCylinder(r=r) for r in good_radii] # Bad material type with pytest.raises(TypeError): - Pin.from_radii(good_radii, [mat.name for mat in pin_mats]) + pin(good_surfaces, [mat.name for mat in pin_mats]) # Incorrect lengths - with pytest.raises(ValueError, match="length") as exec_info: - Pin.from_radii(good_radii[: len(pin_mats) - 2], pin_mats) + with pytest.raises(ValueError, match="length"): + pin(good_surfaces[:len(pin_mats) - 2], pin_mats) # Non-positive radii - rad = (-0.1,) + good_radii[1:] - with pytest.raises(ValueError, match="index 0") as exec_info: - Pin.from_radii(rad, pin_mats) + rad = [openmc.ZCylinder(r=-0.1)] + good_surfaces[1:] + with pytest.raises(ValueError, match="index 0"): + pin(rad, pin_mats) # Non-increasing radii - rad = tuple(reversed(good_radii)) - with pytest.raises(ValueError, match="index 1") as exec_info: - Pin.from_radii(rad, pin_mats) + surfs = tuple(reversed(good_surfaces)) + with pytest.raises(ValueError, match="index 1"): + pin(surfs, pin_mats) # Bad orientation - with pytest.raises(ValueError, match="Orientation") as exec_info: - Pin.from_radii(good_radii, pin_mats, orientation="fail") + surfs = [openmc.XCylinder(r=good_surfaces[0].r)] + good_surfaces[1:] + with pytest.raises(TypeError, match="surfaces"): + pin(surfs, pin_mats) -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) +@pytest.mark.parametrize( + "surf_type", [openmc.ZCylinder, openmc.XCylinder, openmc.YCylinder]) +def test_subdivide(pin_mats, good_radii, surf_type): + """Test the subdivision with various orientations""" + surfs = [surf_type(r=r) for r in good_radii] + fresh = pin(surfs, pin_mats) + assert len(fresh.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 + div0 = pin(surfs, pin_mats, {0: N}) + assert len(div0.cells) == len(pin_mats) + N - 1 + # check volumes of new rings - bounds = (0,) + pin.radii[:N] + radii = get_pin_radii(div0) + bounds = [0] + 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 + new_pin = pin(surfs, pin_mats, {1: N}) 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) + radii = get_pin_radii(new_pin) + sqrs = numpy.square(radii[:N + 1]) assert sqrs[1:] - sqrs[:-1] == pytest.approx( (good_radii[1] ** 2 - good_radii[0] ** 2) / N)