Fix a bug in rotational periodic boundary conditions

include/openmc/boundary_condition.h

@@ -152,6 +152,8 @@ class RotationalPeriodicBC : public PeriodicBC {
 protected:
   //! Angle about the axis by which particle coordinates will be rotated
   double angle_;
+  //! Do we need to flip surfaces senses when applying the transformation?
+  bool flip_sense_;
   //! Ensure that choice of axes is right handed. axis_1_idx_ corresponds to the
   //! independent axis and axis_2_idx_ corresponds to the dependent axis in the
   //! 2D plane  perpendicular to the planes' axis of rotation

src/boundary_condition.cpp

@@ -173,14 +173,9 @@ RotationalPeriodicBC::RotationalPeriodicBC(
     axis_2_idx_ = 2;    // z component dependent
     break;
   case y:
-    // for a right handed coordinate system, z should be the independent axis
-    // but this would cause the y-rotation case to be different than the other
-    // two. using a left handed coordinate system and a negative rotation the
-    // compute angle and rotation matrix behavior mimics that of the x and z
-    // cases
     zero_axis_idx_ = 1; // y component of plane must be zero
-    axis_1_idx_ = 0;    // x component independent
-    axis_2_idx_ = 2;    // z component dependent
+    axis_1_idx_ = 2;    // z component independent
+    axis_2_idx_ = 0;    // x component dependent
     break;
   case z:
     zero_axis_idx_ = 2; // z component of plane must be zero
@@ -192,6 +187,9 @@ RotationalPeriodicBC::RotationalPeriodicBC(
       fmt::format("You've specified an axis that is not x, y, or z."));
   }
 
+  Direction ax = {0.0, 0.0, 0.0};
+  ax[zero_axis_idx_] = 1.0;
+
   // Compute the surface normal vectors and make sure they are perpendicular
   // to the correct axis
   Direction norm1 = surf1.normal({0, 0, 0});
@@ -212,8 +210,20 @@ RotationalPeriodicBC::RotationalPeriodicBC(
       surf2.id_));
   }
 
-  angle_ = compute_periodic_rotation(norm1[axis_2_idx_], norm1[axis_1_idx_],
-    norm2[axis_2_idx_], norm2[axis_1_idx_]);
+  // Compute the signed rotation angle about the periodic axis. Note that
+  // (n1×n2)·a = |n1||n2|sin(θ) and n1·n2 = |n1||n2|cos(θ), where a is the axis
+  // of rotation.
+  auto c = norm1.cross(norm2);
+  angle_ = std::atan2(c.dot(ax), norm1.dot(norm2));
+
+  // If the normals point in the same general direction, the surface sense
+  // should change when crossing the boundary
+  flip_sense_ = (norm1.dot(norm2) > 0.0);
+  if (!flip_sense_)
+    angle_ += PI;
+
+  // Normalize range of angle to [-PI,PI].
+  angle_ = std::remainder(angle_, 2 * PI);
 
   // Warn the user if the angle does not evenly divide a circle
   double rem = std::abs(std::remainder((2 * PI / angle_), 1.0));
@@ -225,47 +235,18 @@ RotationalPeriodicBC::RotationalPeriodicBC(
   }
 }
 
-double RotationalPeriodicBC::compute_periodic_rotation(
-  double rise_1, double run_1, double rise_2, double run_2) const
-{
-  // Compute the BC rotation angle.  Here it is assumed that both surface
-  // normal vectors point inwards---towards the valid geometry region.
-  // Consequently, the rotation angle is not the difference between the two
-  // normals, but is instead the difference between one normal and one
-  // anti-normal.  (An incident ray on one surface must be an outgoing ray on
-  // the other surface after rotation hence the anti-normal.)
-  double theta1 = std::atan2(rise_1, run_1);
-  double theta2 = std::atan2(rise_2, run_2) + PI;
-  return theta2 - theta1;
-}
-
 void RotationalPeriodicBC::handle_particle(
   Particle& p, const Surface& surf) const
 {
-  int i_particle_surf = p.surface_index();
-
-  // Figure out which of the two BC surfaces were struck to figure out if a
-  // forward or backward rotation is required.  Specify the other surface as
-  // the particle's new surface.
-  double theta;
-  int new_surface;
-  if (i_particle_surf == i_surf_) {
-    theta = angle_;
-    new_surface = p.surface() > 0 ? -(j_surf_ + 1) : j_surf_ + 1;
-  } else if (i_particle_surf == j_surf_) {
-    theta = -angle_;
-    new_surface = p.surface() > 0 ? -(i_surf_ + 1) : i_surf_ + 1;
-  } else {
-    throw std::runtime_error(
-      "Called BoundaryCondition::handle_particle after "
-      "hitting a surface, but that surface is not recognized by the BC.");
-  }
+  int new_surface = p.surface() > 0 ? -(j_surf_ + 1) : j_surf_ + 1;
+  if (flip_sense_)
+    new_surface = -new_surface;
 
-  // Rotate the particle's position and direction about the z-axis.
+  // Rotate the particle's position and direction.
   Position r = p.r();
   Direction u = p.u();
-  double cos_theta = std::cos(theta);
-  double sin_theta = std::sin(theta);
+  double cos_theta = std::cos(angle_);
+  double sin_theta = std::sin(angle_);
 
   Position new_r;
   new_r[zero_axis_idx_] = r[zero_axis_idx_];

src/surface.cpp

@@ -1371,7 +1371,7 @@ void read_surfaces(pugi::xml_node node)
           "axis, which is not supported."));
       }
       surf1.bc_ = make_unique<RotationalPeriodicBC>(i_surf, j_surf, axis);
-      surf2.bc_ = make_unique<RotationalPeriodicBC>(i_surf, j_surf, axis);
+      surf2.bc_ = make_unique<RotationalPeriodicBC>(j_surf, i_surf, axis);
     }
 
     // If albedo data is present in albedo map, set the boundary albedo.

tests/regression_tests/periodic_6fold/False-True/inputs_true.dat

@@ -0,0 +1,34 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="5">
+      <density value="1.0" units="g/cc"/>
+      <nuclide name="H1" ao="2.0"/>
+      <nuclide name="O16" ao="1.0"/>
+      <sab name="c_H_in_H2O"/>
+    </material>
+    <material id="6" depletable="true">
+      <density value="4.5" units="g/cc"/>
+      <nuclide name="U235" ao="1.0"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="5" region="9 -10 -11 12" universe="0"/>
+    <cell id="2" material="6" region="9 -10 -12" universe="0"/>
+    <surface id="9" type="plane" boundary="periodic" coeffs="0.4999999999999999 0.8660254037844387 0.0 0.0"/>
+    <surface id="10" type="plane" boundary="periodic" coeffs="-0.4999999999999999 0.8660254037844387 0.0 0.0"/>
+    <surface id="11" type="x-plane" boundary="reflective" coeffs="5.0"/>
+    <surface id="12" type="z-cylinder" coeffs="2.598076211353316 1.4999999999999998 2.0"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>1000</particles>
+    <batches>4</batches>
+    <inactive>0</inactive>
+    <source type="independent" strength="1.0" particle="neutron">
+      <space type="box">
+        <parameters>0 0 0 5 5 0</parameters>
+      </space>
+    </source>
+  </settings>
+  </model>

tests/regression_tests/periodic_6fold/False-True/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+1.848492E+00 2.933785E-03

tests/regression_tests/periodic_6fold/True-False/inputs_true.dat

@@ -0,0 +1,34 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="3">
+      <density value="1.0" units="g/cc"/>
+      <nuclide name="H1" ao="2.0"/>
+      <nuclide name="O16" ao="1.0"/>
+      <sab name="c_H_in_H2O"/>
+    </material>
+    <material id="4" depletable="true">
+      <density value="4.5" units="g/cc"/>
+      <nuclide name="U235" ao="1.0"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="3" region="-5 6 -7 8" universe="0"/>
+    <cell id="2" material="4" region="-5 6 -8" universe="0"/>
+    <surface id="5" type="plane" boundary="periodic" coeffs="-0.4999999999999999 -0.8660254037844387 0.0 0.0"/>
+    <surface id="6" type="plane" boundary="periodic" coeffs="0.4999999999999999 -0.8660254037844387 0.0 0.0"/>
+    <surface id="7" type="x-plane" boundary="reflective" coeffs="5.0"/>
+    <surface id="8" type="z-cylinder" coeffs="2.598076211353316 1.4999999999999998 2.0"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>1000</particles>
+    <batches>4</batches>
+    <inactive>0</inactive>
+    <source type="independent" strength="1.0" particle="neutron">
+      <space type="box">
+        <parameters>0 0 0 5 5 0</parameters>
+      </space>
+    </source>
+  </settings>
+  </model>

tests/regression_tests/periodic_6fold/True-False/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+1.848492E+00 2.933785E-03

tests/regression_tests/periodic_6fold/True-True/inputs_true.dat

@@ -0,0 +1,34 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="7">
+      <density value="1.0" units="g/cc"/>
+      <nuclide name="H1" ao="2.0"/>
+      <nuclide name="O16" ao="1.0"/>
+      <sab name="c_H_in_H2O"/>
+    </material>
+    <material id="8" depletable="true">
+      <density value="4.5" units="g/cc"/>
+      <nuclide name="U235" ao="1.0"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="7" region="-13 -14 -15 16" universe="0"/>
+    <cell id="2" material="8" region="-13 -14 -16" universe="0"/>
+    <surface id="13" type="plane" boundary="periodic" coeffs="-0.4999999999999999 -0.8660254037844387 0.0 0.0"/>
+    <surface id="14" type="plane" boundary="periodic" coeffs="-0.4999999999999999 0.8660254037844387 0.0 0.0"/>
+    <surface id="15" type="x-plane" boundary="reflective" coeffs="5.0"/>
+    <surface id="16" type="z-cylinder" coeffs="2.598076211353316 1.4999999999999998 2.0"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>1000</particles>
+    <batches>4</batches>
+    <inactive>0</inactive>
+    <source type="independent" strength="1.0" particle="neutron">
+      <space type="box">
+        <parameters>0 0 0 5 5 0</parameters>
+      </space>
+    </source>
+  </settings>
+  </model>

tests/regression_tests/periodic_6fold/True-True/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+1.848492E+00 2.933785E-03

tests/regression_tests/periodic_6fold/test.py

@@ -1,13 +1,15 @@
 from math import sin, cos, pi
 
 import openmc
+from openmc.utility_funcs import change_directory
 import pytest
 
 from tests.testing_harness import PyAPITestHarness
 
 
-@pytest.fixture
-def model():
+@pytest.mark.parametrize("flip1", [False, True])
+@pytest.mark.parametrize("flip2", [False, True])
+def test_periodic(flip1, flip2):
     model = openmc.model.Model()
 
     # Define materials
@@ -27,17 +29,40 @@ def model():
     # answers.
     theta1 = (-1/6 + 1/2) * pi
     theta2 = (1/6 - 1/2) * pi
-    plane1 = openmc.Plane(a=cos(theta1), b=sin(theta1), boundary_type='periodic')
-    plane2 = openmc.Plane(a=cos(theta2), b=sin(theta2), boundary_type='periodic')
+    if flip1:
+        plane1 = openmc.Plane(a=-cos(theta1), b=-sin(theta1), boundary_type='periodic')
+    else:
+        plane1 = openmc.Plane(a=cos(theta1), b=sin(theta1), boundary_type='periodic')
+    if flip2:
+        plane2 = openmc.Plane(a=-cos(theta2), b=-sin(theta2), boundary_type='periodic')
+    else:
+        plane2 = openmc.Plane(a=cos(theta2), b=sin(theta2), boundary_type='periodic')
 
     x_max = openmc.XPlane(5., boundary_type='reflective')
 
     z_cyl = openmc.ZCylinder(x0=3*cos(pi/6), y0=3*sin(pi/6), r=2.0)
-
-    outside_cyl = openmc.Cell(1, fill=water, region=(
-        +plane1 & +plane2 & -x_max & +z_cyl))
-    inside_cyl = openmc.Cell(2, fill=fuel, region=(
-        +plane1 & +plane2 & -z_cyl))
+
+    match (flip1,flip2):
+        case (False,False):
+            outside_cyl = openmc.Cell(1, fill=water, region=(
+                          +plane1 & +plane2 & -x_max & +z_cyl))
+            inside_cyl = openmc.Cell(2, fill=fuel, region=(
+                          +plane1 & +plane2 & -z_cyl))
+        case (False,True):
+            outside_cyl = openmc.Cell(1, fill=water, region=(
+                          +plane1 & -plane2 & -x_max & +z_cyl))
+            inside_cyl = openmc.Cell(2, fill=fuel, region=(
+                          +plane1 & -plane2 & -z_cyl))        
+        case (True,False):
+            outside_cyl = openmc.Cell(1, fill=water, region=(
+                          -plane1 & +plane2 & -x_max & +z_cyl))
+            inside_cyl = openmc.Cell(2, fill=fuel, region=(
+                          -plane1 & +plane2 & -z_cyl))        
+        case (True,True):
+            outside_cyl = openmc.Cell(1, fill=water, region=(
+                          -plane1 & -plane2 & -x_max & +z_cyl))
+            inside_cyl = openmc.Cell(2, fill=fuel, region=(
+                          -plane1 & -plane2 & -z_cyl))        
     root_universe = openmc.Universe(0, cells=(outside_cyl, inside_cyl))
     model.geometry = openmc.Geometry(root_universe)
 
@@ -49,9 +74,7 @@ def model():
     model.settings.source = openmc.IndependentSource(space=openmc.stats.Box(
         (0, 0, 0), (5, 5, 0))
     )
-    return model
-
+    with change_directory(f'{flip1}-{flip2}'):
+        harness = PyAPITestHarness('statepoint.4.h5', model)
+        harness.main()
 
-def test_periodic(model):
-    harness = PyAPITestHarness('statepoint.4.h5', model)
-    harness.main()

tests/unit_tests/test_periodic_bc.py

@@ -0,0 +1,47 @@
+from math import cos, sin, radians
+import random
+
+import openmc
+import pytest
+
+
+@pytest.mark.parametrize("angle", [30., 45., 60., 90., 120.])
+def test_rotational_periodic_bc(angle):
+    # Pick random starting angle
+    start = random.uniform(0., 360.)
+    degrees = angle
+    ang1 = radians(start)
+    ang2 = radians(start + degrees)
+
+    # Define three points on each plane and then randomly shuffle them
+    p1_points = [(0., 0., 0.), (cos(ang1), sin(ang1), 0.), (0., 0., 1.)]
+    p2_points = [(0., 0., 0.), (cos(ang2), sin(ang2), 0.), (0., 0., 1.)]
+    random.shuffle(p1_points)
+    random.shuffle(p2_points)
+
+    # Create periodic planes and a cylinder
+    p1 = openmc.Plane.from_points(*p1_points, boundary_type='periodic')
+    p2 = openmc.Plane.from_points(*p2_points, boundary_type='periodic')
+    p1.periodic_surface = p2
+    zcyl = openmc.ZCylinder(r=5., boundary_type='vacuum')
+
+    # Figure out which side of planes to use based on a point in the middle
+    ang_mid = radians(start + degrees/2.)
+    mid_point = (cos(ang_mid), sin(ang_mid), 0.)
+    r1 = -p1 if mid_point in -p1 else +p1
+    r2 = -p2 if mid_point in -p2 else +p2
+
+    # Create one cell bounded by the two planes and the cylinder
+    mat = openmc.Material(density=1.0, density_units='g/cm3', components={'U235': 1.0})
+    cell = openmc.Cell(fill=mat, region=r1 & r2 & -zcyl)
+
+    # Make the model complete
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([cell])
+    model.settings.source = openmc.IndependentSource(space=openmc.stats.Point(mid_point))
+    model.settings.particles = 1000
+    model.settings.batches = 10
+    model.settings.inactive = 5
+
+    # Run the model
+    model.run()