Logical evaluation of surfaces for non-simple cells

[paulromano]

OpenMC currently only supports simple cells, i.e. cells defined as the intersection of multiple half-spaces. For ease of geometric modeling, we may want to add support for the union and difference operators.

[paulromano]

@wbinventor @smharper Reviving discussion on an ancient issue. Curious to know what both of you think about whether we ought to support non-simple cells, and if so how you would envision the implementation working. As a reference point, MCNP allows you to put a union (:) or complement (#) operator in the surface notation, so you could have something like 1 -2 : 3 #(4 -5) which would read as "the union of (1 intersected with -2) and (the intersection of 3 with the complement of (4 intersected with -5))". Note that the order of operations is complement, intersection, union. And, as in OpenMC, intersection is implicit via whitespace.

I can see two ways of going about this:

1. Full support in OpenMC with a notation akin to that used in MCNP. Pros: gives the greatest flexibility in modeling and defining complex regions over which tallies can be collected. Cons: Make distance-to-surface calculations more complex.
2. Support in OpenCG with conversion of non-simple cells to simple cells. I started thinking about the logic for this, and it becomes a bit complex. For example, if you had the union of five half-spaces there are up to 31 unique simple cells -- an easy way is to determine this is to look at a five-way Venn diagram:

[wbinventor]

@paulromano - union/complement operations would be very nice to have. From a reactor modeling standpoint, some of the objects that would be much easier to construct with these features include the following:

• SFR duct, interstitial coolant
• BWR duct, crucifix CRs
• PWR spacer grids, baffle

Some of these objects can be roughly modeled using the lattice outer universe option. However, I don't find this to be a very satisfying or flexible solution. This is particularly true for some scenarios where one may want to tally in those complex regions with distribcells (i.e., for MGXS or dose tallies you would want the spatially distributed tallies rather than spatially-averaged).

You may have noticed that the OpenCG compatibility module contains the logic to do #2 from your description above for certain types of cells/surfaces. In particular, we have a "parallelepiped" surface in OpenCG with versions aligned along the x, y, and z axes (akin to cylinders in OpenMC). If one designs a Cell that contains the "outer" halfspace of a parallelepiped, the OpenCG compatibility module will convert that single Cell into 8 different cells which together comprise the outside of the parallelepiped (together with any other bounding surfaces there may be). Although this works and is used for the grid spacers in our OpenCG BEAVRS benchmark, I certainly don't think it is a robust path forward due to the complexities that you alluded to in your Venn diagram.

Hence, I would vote in favor of #1 and incorporate cell unions/complements in OpenMC directly. One alternative (or perhaps complementary) approach would be to add more complex surface types to OpenMC, such as parallelepipeds and hexagons. Of course these would not be quadratic surfaces, but that shouldn't matter as long as one can determine the distance to each surface type given (x,y,z). This would resolve most if not all of the scenarios I envisioned above.

That being said, for more complex geometries, I think union/complements surfaces would be necessary. I have thought about how to do this for OpenMOC and OpenCG before and I think it would be fairly simple with the right data structures. Even so, it is worth considering that it may be more worthwhile to focus on geometric abstraction per #438 to enable the use of more sophisticated (and perhaps more optimized, though I would be skeptical) geometry modeling libraries to do this and more for us.

[paulromano]

I was leaning in that direction as well. If we're going to go through the hassle and complexity of supporting it one way or another, might as well have it directly in OpenMC. Now all we need to do is convince @smharper to implement it after he's don'e with #438 and #214! :wink:

[smharper]

I agree that we should probably have support for these complex cells implemented directly in OpenMC. I think if we did that, a natural evolution would be to offer structures like hex ducts which would use the same complex cell code on the backend.

That said, this is pretty low on my priority list. Maybe it would make a good UROP project?

[wbinventor]

You would be the right person to find one with an interest in learning Fortran :-) On Aug 17, 2015 5:30 PM, "Sterling Harper" notifications@github.com wrote:

I agree that we should probably have support for these complex cells implemented directly in OpenMC. I think if we did that, a natural evolution would be to offer structures like hex ducts which would use the same complex cell code on the backend.

That said, this is pretty low on my priority list. Maybe it would make a good UROP project?

— Reply to this email directly or view it on GitHub https://github.com/mit-crpg/openmc/issues/60#issuecomment-132004557.

[paulromano]

I took some time this week to implement this and I now have what I think to be a working prototype. Turned out to be easier than I thought it would be. Along the way, I also refactored how surfaces are handled (each type of surface now has its own derived type that inherits from an abstract Surface type) and it has considerably cleaned up the geometry. For user input, I've chosen ~ as the complement operator and ^ for the union operator, which I think are a bit more intuitive than the MCNP-equivalents # and :.

@wbinventor @smharper I have two questions I'd like input on:

1. Should we support the MCNP shorthand where the complement directly followed by a number x means "the complement of cell x's surface specification". Normally, the complement operator is followed by an expression in parentheses.
2. How should the Python API handle this newfound flexibility? Right now, we only have Cell.add_surface(surface, halfspace); nowhere is there any concept of operators.

[wbinventor]

@paulromano it is cool to hear that you've had some time to work on this, this functionality will make certain types of geometries much easier to build. I've never used complements and unions in MCNP so I don't have some historical preference to follow their operator notation - I agree that it is best to follow the standard mathematical notation of ^ and ~.

Before replying to your questions, I have a question of my own - it seems that the union operator could not only be applied between two or more surface halfspaces, but also between cells themselves. Thoughts?

1. I like the idea of putting the expression in parentheses. The shorthand notation doesn't seem relevant when one has a Python API to generate the input files, but that's just my opinion.
2. What about adding some new surface types to the API (perhaps SurfaceUnion and SurfaceComplement) akin to the CrossFilter, CrossScore and CrossNuclide used for tally arithmetic? For example, a user might define the outside of a square (i.e., for a spacer grid) using a SurfaceUnion as follows:

# Instantiate planes to bound a pin cell akin to a spacer grid
left = openmc.XPlane(x0=-0.62992, name='left')
right = openmc.XPlane(x0=0.62992, name='right')
bottom = openmc.YPlane( y0=-0.62992, name='bottom')
top = openmc.YPlane(y0=0.62992, name='top')

# Instantiate a SurfaceUnion and add surface halfspaces 
# to represent the outside of the box
box = openmc.SurfaceUnion(name='spacer grid box')
box.add_surface(surface=left, halfspace=-1)
box.add_surface(surface=right, halfspace=+1)
box.add_surface(surface=top, halfspace=+1)
box.add_surface(surface=bottom, halfspace=-1)

# Instantiate cells for fuel, clad, water, what have you
...

# Create cell for spacer grid box
grid = openmc.Cell(cell_id=1, name='grid')
grid.fill = zircaloy

# Add the surface union representing the outside of the box
grid.add_surface(surface=box, halfspace=+1)

# Do more stuff
...

The above example sketches out how one might go about creating the outside of a box to represent a spacer grid (of course one would also need to represent the inside of a larger box using standard surface halfspace intersections or a lattice). The SurfaceUnion could be added to Cells just as a standard Surface might be. I'm not sure what the halfspace of a SurfaceUnion might represent in this scenario, but I'll give it some more thought (perhaps just use it to flip the signs for each of the surface halfspaces in the SurfaceUnion, but I'm not sure what good this would do).

Finally, I think a SurfaceComplement might follow the same sort of scheme. What do you think?

[paulromano]

I tend to agree that given a Python API, the complement shorthand doesn't add much. And internally, OpenMC would just expand the shorthand anyway so it doesn't give a computational advantage.

I think your proposal is a step in the right direction. To me, it seems like we're still missing an abstraction for "a region of space" which is not exactly the same as a cell, rather one level below it. It might be clearer if we have a Volume object that is created from surface half-spaces. We'd need a way to go from a Surface to a Volume, and a way to join multiple volumes together. For Surface -> Volume, we could:

• Have a factory @staticmethod on the Volume class: negative_halfspace = Volume.halfspace(surface1, side='-')
• Have a @property of Surface which generates halfspaces negative_halfspace = surface1.negative

For joining multiple volumes, we could:

• Have a method for each operator, e.g. v3 = v1.intersect(v2)
• Have a factory @staticmethod on the Volume class: v3 = Volume.intersection(v1, v2)

There are probably more ways you could go about it. Just throwing a few off the top of my head. We'd also need a new method/property on the Cell object to associate a volume with it.

So, to repeat your example but with something like what I'm proposing

# Instantiate planes to bound a pin cell akin to a spacer grid
left = openmc.XPlane(x0=-0.62992, name='left')
right = openmc.XPlane(x0=0.62992, name='right')
bottom = openmc.YPlane( y0=-0.62992, name='bottom')
top = openmc.YPlane(y0=0.62992, name='top')

# Create a volume from four half-spaces
box = openmc.Volume.union(left.negative, right.positive, top.negative, bottom.positive)

# Instantiate cells for fuel, clad, water, what have you
...

# Create cell for spacer grid box
grid = openmc.Cell(cell_id=1, name='grid')
grid.fill = zircaloy

# Add the surface union representing the outside of the box
grid.assign_region(box)

# Do more stuff
...

One potential downside is that the add_surface method that exists today won't really make sense. We could keep it around but make it deprecated, i.e., if a user tries to call add_surface, the default behavior will be to take the intersection of the currently assigned volume and whatever was passed.

One more thought -- maybe "Region" is less ambiguous "Volume". One might like to have a method Region.get_volume() that returns the volume in cm^3.

[wbinventor]

I need to give this some more thought, but I like the direction you're proposing here. I do think that Region would be more appropriate, but am still a bit unclear on how exactly this would be different than a Cell (and think many users might find this distinction confusing without very good documentation).

[paulromano]

The distinction is that a region specifically only refers to a region of space. A cell then associates with that region a material, fill universe, or lattice. Other physical properties may be associated with the cell too (e.g., temperature) in the future. I personally find this abstraction more compelling then the "add_surface" method where it is implicit that intersection between half-spaces occurs. The notation above forces the user to say that they want the intersection between surface half-spaces.

The other thing I like about it is that the defined region is stored as a CSG tree with Volume objects as its nodes. This makes traversing it very simple if we ever want to add capability in the Python API to do something with it. The implementation would look something like

class Region(object):
    def __init__(self, operator, nodes):
        self.operator = operator
        self.nodes = nodes

    @staticmethod
    def intersection(cls, *nodes):
        return cls('intersection', nodes)
    ...

[wbinventor]

Ok, good explanation, I agree that this distinction is a good one and would simplify the inputs for more complex geometries. Are you thinking that the Region would simply be a construct provided by the Python API, or also in the back-end with its own XML input spec? On Sep 25, 2015 4:21 AM, "Paul Romano" notifications@github.com wrote:

The distinction is that a region specifically only refers to a region of space. A cell then associates with that region a material, fill universe, or lattice. Other physical properties may be associated with the cell too (e.g., temperature) in the future. I personally find this abstraction more compelling then the "add_surface" method where it is implicit that intersection between half-spaces occurs. The notation above forces the user to say that they want the intersection between surface half-spaces.

The other thing I like about it is that the defined region is stored as a CSG tree with Volume objects as its nodes. This makes traversing it very simple if we ever want to add capability in the Python API to do something with it. The implementation would look something like

class Region(object): def init(self, operator, nodes): self.operator = operator self.nodes = nodes

@staticmethod
def intersection(cls, *nodes):
    return cls('intersection', nodes)
...

— Reply to this email directly or view it on GitHub https://github.com/mit-crpg/openmc/issues/60#issuecomment-143157368.

[paulromano]

I coded up an implementation quickly and now suggest the following layout:

class Region(object):
    __metaclass__ = ABCMeta

    @abstractmethod
    def __str__(self):
        return ''

class Intersection(Region):
    def __init__(self, *nodes):
        self.nodes = nodes

    def __str__(self):
        return '(' + ' '.join(map(str, self.nodes)) + ')'

class Union(Region):
    def __init__(self, *nodes):
        self.nodes = nodes

    def __str__(self):
        return '(' + ' ^ '.join(map(str, self.nodes)) + ')'

class Complement(Region):
    def __init__(self, node):
        self.node = node

    def __str__(self):
        return '~' + str(self.node)

class Halfspace(Region):
    def __init__(self, surface, side):
        self.surface = surface
        self.side = side

    def __str__(self):
        return '-' + str(self.surface.id) if self.side == '-' \
            else str(self.surface.id)

def Surface(object):
    ...
    @property
    def negative(self):
        return Halfspace(self, '-')
    ...

The beauty of the abstraction is that the __str__() methods can recursively go through the CSG tree to generate the proper infix notation to pass in <cell surfaces="..." />. Or, to give an example:

In [1]: import openmc

In [2]: from openmc.region import Intersection, Union, Complement

In [3]: x1 = openmc.XPlane(1, x0=-5)

In [4]: x2 = openmc.XPlane(2, x0=5)

In [5]: x3 = openmc.YPlane(3, y0=-5)

In [6]: x4 = openmc.YPlane(4, y0=5)

In [7]: x5 = openmc.ZPlane(5, z0=-5)

In [8]: vol = Union(Intersection(x1.positive, x2.negative),
   ...: Intersection(x3.positive, Complement(Intersection(x4.positive, x5.negative))))

In [9]: str(vol)
Out[9]: '((1 -2) ^ (3 ~(4 -5)))'

I'm thinking the Region construct would only be in the Python API, with that last string being given to OpenMC. It would probably make sense to rename <cell surfaces="..." /> to <cell region="..." />. We could still accept surfaces but just make it deprecated and eventually remove it at some point in the future.

[wbinventor]

I like your implementation sketch a lot @paulromano! And I agree that (at least for now) the XML notation should be kept as is and the region/volume delineation can simply be a Python construct.

[smharper]

Closed with #463