Refactor and generalize the boundary condition implementation

[smharper]

This PR gives an overhaul of the core boundary condition software. Mostly, this is to clean-up and generalize the periodic BCs.

The old periodic BC implementation was kind of messy. The same PeriodicSurface::periodic_translate function was used for both rotational and translational BCs, and this function would redundantly check during particle transport which type of periodicity was relevant. The old implementation also made some unchecked assumptions about the types of surfaces that were paired via these BC's, and it relied on several dynamic_cast calls during particle transport.

This PR cleans up much of that by moving most of the BC logic into a new BoundaryCondition class. Surfaces are now assigned a pointer to a BoundaryCondition object, and the different BC types are handled by different virtual functions rather than an enum switch-case statement. The periodic BCs are split into the distinct TranslationalPeriodicBC and RotationalPeriodicBC classes.

These changes generalize the periodic BCs a bit. We can now handle rotations with the general SurfacePlane object instead of just SurfaceXPlane/SurfaceYPlane pairs. This means we can now model 3-fold, 6-fold, etc. rotational periodicity as well as 4-fold. (Although it is still limited to rotations about the z-axis at x=0, y=0). For periodic BCs, there is also now a lot more checking done during initialization to make sure the user hasn't specified some un-expected/non-implemented case.

[gridley]

Tiny note: you don't need to explicitly initialize any of the C++ smart pointers to nullptr, since that is automatic and already in their default constructors.

[smharper]

Tiny note: you don't need to explicitly initialize any of the C++ smart pointers to nullptr, since that is automatic and already in their default constructors.

Fair, but I was thinking it'd be good to be explicit about this and give a sort of hint to other other devs that nullptr is an expected value for that pointer and must be handled. I guess I'm drawing a bit on the Zen of Python here.

[gridley]

Fair enough. It's stylistic, end of the day. There is a fairly different zen written by a former lead developer of MCNP I have found myself subscribing to lately, although I can't remember where I once read it. The essence of it is that if you write your code concisely such that the whole of the algorithm you're writing fits on the screen all at once, it's more clear to see the overall flow of the code. But then again... is anything written in fortran worth looking up to?

[paulromano]

Very nice reorganization @smharper, and it's great that this also adds support for generic rotational symmetry. Here's my first cut at a review -- I'm going to do some sanity checks with simple models as well so I might have another round of feedback.

[paulromano]

Using a lambda might make this a bit cleaner

[paulromano]

Please add some comments explaining the logic here

[paulromano]

If we were to store the pairs of values such that the lesser (or greater) value comes first, you could use a set and thereby avoid the elimination of duplicates further down

[paulromano]

Probably can remove this?

[paulromano]

It looks like this doesn't work for all cases. I tried doing a 1/6th rotational symmetry problem with the following:

plane1 = openmc.YPlane(boundary_type='periodic')

angle = pi/3
p1 = (0., 0., 0.)
p2 = (cos(angle), sin(angle), 0.)
p3 = (0., 0., 1.)
plane2 = openmc.Plane.from_points(p1, p2, p3, boundary_type='periodic')
plane2.periodic_surface = plane1

Note that for plane 2, the normal is (p2 - p1) x (p3 - p1) which results in a direction that is down and to the right. Here is what you end up with (planes in green, normals in blue):

As you can see, we end up with an angle between the normals that is 2π/3 instead of π/3 as desired. The end result is that when you try to run a model like this, all particles that cross the periodic boundary get lost. If you swap p2 and p3 in the definition of the plane, the normal will go in the other direction, the calculated angle is π/3, and everything works as expected.

[paulromano]

Hmm, thinking about this a little more, if instead of switching the planes, I use the region below the y-plane (which gives me 1/3rd rotational symmetry), the simulation works fine because in that case the calculated angle matches the region. So this implies that how the angle is calculated depends on how the surfaces are actually used to create a cell, which seems awkward. Maybe we need some sort of convention for how a model must be constructed if rotational periodic BCs are used?

[smharper]

I haven't figured out a good solution to that problem yet. When I wrote this code, I was operating off the assumption that one normal vector would normally point inward, and the other outward. (This is likely what happens if you have an azimuthally segmented problem and you define all the planes from the general Plane class in a loop, as was done in the problem that inspired this PR.) But this assumption seems to be false often enough.

Perhaps we'll have to specify a convention, but I'll keep thinking to see if we can come up with something easy. Maybe we can compute the bounding box of the root universe and check it against the normals of the bounding planes.

[paulromano]

Sounds good -- let me know if you come up with something clever. The bounding box idea sounds a little complicated, so if it's too difficult we can always just adopt some convention that is expected to be obeyed. I think as long as it is well documented, then at least we have something we can point to if (and when) a user messes it up.

[smharper]

I've abandoned the attempt now to try and resolve where the geometry lies with respect to the periodic normal vectors. Here it'll always be assumed that the normals point inwards towards the geometry (like with an XPlane and YPlane in the first quadrant).

I've added some notes to the docs to explain this, and I've also added a note about normals to the lost particle message to hint users on that point.

[paulromano]

@smharper Sorry for the delay here -- I tried my 1/6th rotational periodic case again and it still doesn't work (even though I believe my normals are inward facing). It looks like the calculation of the angle is still the same, in which case it would result in an angle of 2π/3 instead of π/3 for this case.

[smharper]

Phew, this geometry problem is surprisingly tricky. I think those last commits should do the trick.

[paulromano]

Thanks for bearing with me. I tested once more on my own 1/6th periodic model and it worked great this time. Answers match an equivalent problem without periodicity.

[smharper]

No worries on the delay. I know how it is! Glad I finally implemented this though. I've wanted to do something like it from the very beginning of the C++ rewrite.

Now I think the only little quibble from the C++ rewrite that I have left is the surface indexing scheme. It's unsatisfying that we always add one to get a signed surface index. I wonder if we can make an efficient class that essentially does the same thing (uses the lower 31 bits for the actual index and then the 32nd bit indicates halfspace) but wraps it in nicer syntax.

[paulromano]

@smharper Yes, I think a simple wrapper class for the surface indexing like what you propose would be great! The off-by-one issue bothers me too.