Feature request: Arbitrary rotation angle for periodic boundary conditions

[lewisgross1296]

Description

I've gone back in time a bit too see the development of periodic BCs in OpenMC (#459, #647, #1129, #1820 ). It seems like there is only support for increments of 90 degrees WRT rotational periodic boundary conditions. Unfortunately, the model I'm working on has 3 way symmetry. I've included an image

Alternatives

It seems like there are no alternatives, which can be especially hampering for hexagonal core shapes and their natural tendency to be symmetric in increments of 30 or 60 degrees

Compatibility

I think the API could stay the same, but that burdens the developers to extend periodic BC code to arbitrary angles. It might be faster to add another API call that accomplishes this and doesn't break the existing rotational symmetry, but it would likely be more confusing for the user. I'm assuming this is a medium skill/effort change, though I'm not exactly sure what needs to happen to the existing code and the reason for current limitations

[lewisgross1296]

A colleague brought up an alternative that doesn't use a periodic BC, but does accomplish a similar goal of lowering the model size. This alternative only works when you are using clones though. If someone was only interested in transport and not depletion, it wouldn't help. (i.e. if you don't have any universe clones in your model)

I'll share, in case others are curious. For my case, since I want to exploit 3-way symmetry, I divide the model in 3 and look at which compacts are exactly 120 degrees from each other. You can "mimic" rotational periodic BCs by defining a cloned universe for the first compact and re-using that same object in the other two positions. You do this for each ring. I wrote out how it should look for the image above.

For me, I reduce 54 unique materials per layer to 18, which definitely helps with my geometry.xml and materials.xml size.

[gridley]

Are you sure this can't be done with TranslationalPeriodicBC? You can see in the code there is support for a general Plane to be listed.

[lewisgross1296]

I did just glance at include/openmc/boundary_condition.h. I'm seeing both TranslationalPeriodicBC and RotationalPeriodicBC, but it was my understanding that you can only use the Translational BC to go from XPlanes to XPlanes, YPlanes to YPlanes, or ZPlanes to ZPlanes.

Since I needed 120 degree symmetry, I thought I was unable to leverage this. I'm not sure I totally understand the source code here and if it confirms this. I could make an MWE to test, but what I did above is working for now.

[gridley]

Closing because this is supported--the HexagonalPrism class shows how to set this up in its implementation. The periodic_surface attribute just has to be set. In fact, for the example here, you could have used a HexagonalPrism and just had set the periodic attribute to true.

[lewisgross1296]

Closing because this is supported--the HexagonalPrism class shows how to set this up in its implementation. The periodic_surface attribute just has to be set. In fact, for the example here, you could have used a HexagonalPrism and just had set the periodic attribute to true.

I'm not totally sure we're on the same page. Apologies for not articulating what I meant well enough. I've attached a rough sketch that hopefully helps. The periodic boundary works on the surface of the hexagonal prism if you have a full prism, but I was hoping to cutting the model above into a third of the size and apply a periodic boundary to the planes that slice the full hexagon.

In the attached, I want to simulate only the shaded region. Since it is a 120 degree sector, no matter how you slice it, there is a plane oriented with a cardinal direction and a plane that is not an XPlane or YPlane. AFAIK those two planes cannot be paired as a periodic boundary. Do you see what I mean?

[gridley]

Yeah Lewis! This actually should work if you just set the periodic surfaces to each other. The C++ side will automatically build a RotationalPeriodic BC in this case, which will rotate the particle about the Z axis. So in this case the intersection of the two inner planes just has to be on the z axis.

[lewisgross1296]

I was able to get around the need for this before, but I think this issue should be re-opened (or I can make a new issue that's more generally about periodic BCs).

I attempted to pair an openmc.XPlane with a generic openmc.Plane in order to create a 1/6th core model for a horizontal reactor (i.e. coolant flows in the y-direction). Here's what I tried

# hex_lattice definition
# ...
z_cyl = openmc.ZCylinder(r=6,boundary_type="vacuum")
z_min = openmc.ZPlane(z0=0, boundary_type="vacuum")
z_max = openmc.ZPlane(z0=10, boundary_type="vacuum")
z_core_region = -z_cyl & +z_min & -z_max
rotated_region = z_core_region.rotate((-90.0, 0.0, 0.0))
periodic_bounding_xplane = openmc.XPlane(
    x0=0, boundary_type="periodic"
)
periodic_bounding_plane = openmc.Plane(
    a=-np.sqrt(3) / 3,
    b=0.0,
    c=1,
    boundary_type="periodic"
)
periodic_bounding_plane.periodic_surface = periodic_bounding_xplane
periodic_bounding_xplane.periodic_surface = periodic_bounding_plane

cell = openmc.Cell(region=rotated_region &+periodic_bounding_xplane &+periodic_bounding_plane, fill=hex_lattice)
cell.rotation = (-90.0, 0.0, 0.0)

but when I try to load this into the plotter, it complains

Reading geometry XML file...
terminate called after throwing an instance of 'std::invalid_argument'
  what():  Rotational periodic BCs are only supported for rotations about the z-axis, but surface 955436 is not perpendicular to the z-axis.

The above suggestion would work if my reactor was z-aligned and I tried it sucessfully, but it needs to be in the y-direction for use in a multiphysics (Cardinal) simulation. This is due to these lines

openmc/src/boundary_condition.cpp

Lines 161 to 181 in d1df80a

	RotationalPeriodicBC::RotationalPeriodicBC(int i_surf, int j_surf)
	: PeriodicBC(i_surf, j_surf)
	{
	Surface& surf1 {*model::surfaces[i_surf_]};
	Surface& surf2 {*model::surfaces[j_surf_]};
	// Check the type of the first surface
	bool surf1_is_xyplane;
	if (const auto* ptr = dynamic_cast<const SurfaceXPlane*>(&surf1)) {
	surf1_is_xyplane = true;
	} else if (const auto* ptr = dynamic_cast<const SurfaceYPlane*>(&surf1)) {
	surf1_is_xyplane = true;
	} else if (const auto* ptr = dynamic_cast<const SurfacePlane*>(&surf1)) {
	surf1_is_xyplane = false;
	} else {
	throw std::invalid_argument(fmt::format(
	"Surface {} is an invalid type for "
	"rotational periodic BCs. Only x-planes, y-planes, or general planes "
	"(that are perpendicular to z) are supported for these BCs.",
	surf1.id_));
	}

which show that the above recommendation would only work for paring XPlanes, YPlanes, or Plane perpendicular to the z direction as periodic BCs.

Here's a visualization the above "reactor" (changed the BC to make the plot)

As far as I can tell, there's no way to do a periodic BC for a reactor like this (not z-axis aligned). IMO there's no physical reason why my two planes should not be allowed to be a periodic boundary. I think an extension to allow this case and similar ones is reasonable. Curious to hear other takes, but I'd argue that we should support periodicity in any orientation as long the axis of rotations for the periodicity is either x, y or ,z.

I can get started on a PR for this, but let me know if I should open a new issue.