crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 21:55:41 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Function based pin construction

Browse source 
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
This commit is contained in:
Andrew Johnson 2019-07-05 16:51:10 -05:00
parent 48b8b8b69a
commit 50ea1ed02a
No known key found for this signature in database
GPG key ID: 253418E91B7F6FEB
4 changed files with 167 additions and 288 deletions
Download patch file Download diff file Expand all files Collapse all files

1
openmc/model/__init__.py
View file

@ -1,4 +1,3 @@
from .triso import *
from .model import *
from .funcs import *
from .pin import *

130
openmc/model/funcs.py
View file

@ -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)

249
openmc/model/pin.py
View file

@ -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)

75
tests/unit_tests/test_pin.py
View file

@ -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)