Add periodic boundary conditions

[paulromano]

This pull request adds support for periodic boundary conditions on x-, y-, and z-planes. The user simply specifies boundary="periodic" on pairs of surfaces and OpenMC will automatically match surfaces together.

No special support is needed in the Python API.

[samuelshaner]

Do you think there should also be a check in this block to make sure the complementary boundary of a periodic boundary be periodic (i.e. if the x-min surface is set to periodic, the x-max surface must also be periodic)? We do this in OpenMOC since it doesn't make much sense for only one of the surfaces to be periodic.

[paulromano]

Good suggestion; I've just added such a check.

[smharper]

This sounds weird. Maybe "Cannot transfer particle N across periodic surface in a lower universe."

To be really clear, maybe we should also add "Boundary conditions must be applied to universe 0"

[paulromano]

Oh wow, don't know how this error message got messed up. That's definitely not what I intended to write. I'll fix that.

[smharper]

The search for periodic surfaces can be fooled by a geometry that uses translations. For example:

<?xml version='1.0' encoding='utf-8'?>
<geometry>
  <surface id="1" type="x-plane" coeffs="-3" />
  <cell id="1" material="1" region="-1" universe="1" />
  <cell id="2" material="2" region="1" universe="1" />

  <surface id="11" type="z-cylinder" coeffs="0.0 0.0 0.5" />
  <surface id="12" type="x-plane" coeffs="-2." boundary="periodic"  />
  <surface id="13" type="x-plane" coeffs="2." boundary="periodic"  />
  <cell id="11" fill="1" region="-11" translation="3 0 0"/>
  <cell id="12" material="2" region="11 12 -13" />
</geometry>

However, this is a small error that I think people are very unlikely to encounter, and I can't imagine any geometries that are impossible to build with this limitation so I think we can merge it anyway.

[smharper]

I think this could use a little more explanation in the docs. Something like "Periodic boundary conditions can only be applied to x-, y-, and z-planes. Only axis-aligned periodicity is supported, i.e. x-planes can only be paired with x-planes. Specify which planes are periodic and the code will automatically identify which planes are paired together."

[smharper]

I think we will eventually want to support more complex types of periodicity like hexagonal prisms or rotationally periodic octants. This system won't easily scale to those types of geometries. For that reason, I propose we use separate objects to specify BCs. Like this:

<boundary_conditions>
  <reflective surfaces="1 2" />
  <translation_periodic surfaces="3 4" />
  <rotation_periodic surfaces="5 6" />
  <vacuum surfaces="7" />
<boundary_conditions>

[paulromano]

@smharper Let me cogitate on your proposal. I personally like having boundary conditions specified as attributes of a surface, since that's how we ultimately store them. For rotational periodic conditions, I could imagine something like:

<surface id="1" type="x-plane" coeffs="0.0" boundary="periodic" opposite="2" />
<surface id="2" type="y-plane" coeffs="0.0" boundary="periodic" opposite="1" />

The Python side would have the user specify an opposite attribute that expects an instance of Surface. The name doesn't have to be opposite if you can think of something better ("opposite_boundary", "periodic_with", "periodic_surface", etc.).

[smharper]

Okay, I am fine with leaving the BCs as a surface attribute, but I think we should implement that opposite syntax now. That way, we know our syntax won't change when we generalize to non-axis aligned planes and rotational symmetry. And it will dodge that error I pointed out with translated surfaces. As for the name, maybe paired_bc is more specific?

Just a thought for the future, I think we may eventually want to make BCs a separate object, both in the input and the code. I think there's a lot of research that could be done with really complex BCs (e.g. energy-dependent albedos or maybe another coupled neutronics solver), and that would be easier if the BCs are distinct objects.

[wbinventor]

I don't have a strong opinion about this BCs design debate, but do know that basic cartesian axis-aligned periodic BCs as are implemented here will be useful to a few us of looking to graduate soon. Longer term I probably see @smharper's point that BCs might make sense as an abstraction separate from a surface property. In any case, looking to using simple periodic BCs soon!

[paulromano]

@smharper This PR has been updated per your suggestions. The Plane class and its subclasses now have a periodic_surface attribute, and on the XML side a <surface> has a new periodic_surface_id="..." attribute. I thought this naming scheme best represents what the user is actually supposed to supply (a Surface on the Python side and an id on the XML side). Extending what we have here to rotationally symmetric periodic BCs is now just an implementation exercise--the user input shouldn't change.

[paulromano]

I should also mention I updated the test to use the Python API so that it is tested as well.

[smharper]

an -> can

[smharper]

periodic_surface_id also needs to be added to the docs

[paulromano]

@smharper Up-to-date with your comments

[smharper]

Looks good. Merged

[paulromano]

Thanks for the reviews today @smharper!