Generalize RotationalPeriodicBC for X-, Y-, or Z-axis

include/openmc/boundary_condition.h

@@ -138,18 +138,26 @@ class TranslationalPeriodicBC : public PeriodicBC {
 //==============================================================================
 //! A BC that rotates particles about a global axis.
 //
-//! Currently only rotations about the z-axis are supported.
+//! Only rotations about the x, y, and z axes are supported.
 //==============================================================================
 
 class RotationalPeriodicBC : public PeriodicBC {
 public:
-  RotationalPeriodicBC(int i_surf, int j_surf);
-
+  enum PeriodicAxis { x, y, z };
+  RotationalPeriodicBC(int i_surf, int j_surf, PeriodicAxis axis = z);
+  float compute_periodic_rotation(
+    float rise_1, float run_1, float rise_2, float run_2) const;
   void handle_particle(Particle& p, const Surface& surf) const override;
 
 protected:
   //! Angle about the axis by which particle coordinates will be rotated
   double angle_;
+  //! 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
+  int zero_axis_idx;
+  int axis_1_idx_;
+  int axis_2_idx_;
 };
 
 } // namespace openmc

openmc/surface.py

@@ -123,9 +123,8 @@ class Surface(IDManagerMixin, ABC):
     boundary_type : {'transmission', 'vacuum', 'reflective', 'periodic', 'white'}, optional
         Boundary condition that defines the behavior for particles hitting the
         surface. Defaults to transmissive boundary condition where particles
-        freely pass through the surface. Note that periodic boundary conditions
-        can only be applied to x-, y-, and z-planes, and only axis-aligned
-        periodicity is supported.
+        freely pass through the surface. Note that only axis-aligned
+        periodicity is supported periodic around the x-, y-, and z-axes.
     albedo : float, optional
         Albedo of the surfaces as a ratio of particle weight after interaction
         with the surface to the initial weight. Values must be positive. Only

src/boundary_condition.cpp

@@ -158,63 +158,39 @@ void TranslationalPeriodicBC::handle_particle(
 // RotationalPeriodicBC implementation
 //==============================================================================
 
-RotationalPeriodicBC::RotationalPeriodicBC(int i_surf, int j_surf)
+RotationalPeriodicBC::RotationalPeriodicBC(
+  int i_surf, int j_surf, PeriodicAxis axis)
   : 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_));
-  }
-
-  // Check the type of the second surface
-  bool surf2_is_xyplane;
-  if (const auto* ptr = dynamic_cast<const SurfaceXPlane*>(&surf2)) {
-    surf2_is_xyplane = true;
-  } else if (const auto* ptr = dynamic_cast<const SurfaceYPlane*>(&surf2)) {
-    surf2_is_xyplane = true;
-  } else if (const auto* ptr = dynamic_cast<const SurfacePlane*>(&surf2)) {
-    surf2_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.",
-      surf2.id_));
+  // below convention for right handed coordinate system
+  switch (axis) {
+  case x:
+    zero_axis_idx = 0; // x component of plane must be zero
+    axis_1_idx_ = 1;   // y component independent
+    axis_2_idx_ = 2;   // z component dependent
+    break;
+  case y:
+    zero_axis_idx = 1; // y component of plane must be zero
+    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
+    axis_1_idx_ = 0;   // x component independent
+    axis_2_idx_ = 1;   // y component dependent
+    break;
+  default:
+    throw std::invalid_argument(
+      fmt::format("You've specified an axis that is not x, y, or z."));
   }
 
   // Compute the surface normal vectors and make sure they are perpendicular
-  // to the z-axis
+  // to the correct axis
   Direction norm1 = surf1.normal({0, 0, 0});
   Direction norm2 = surf2.normal({0, 0, 0});
-  if (std::abs(norm1.z) > FP_PRECISION) {
-    throw std::invalid_argument(fmt::format(
-      "Rotational periodic BCs are only "
-      "supported for rotations about the z-axis, but surface {} is not "
-      "perpendicular to the z-axis.",
-      surf1.id_));
-  }
-  if (std::abs(norm2.z) > FP_PRECISION) {
-    throw std::invalid_argument(fmt::format(
-      "Rotational periodic BCs are only "
-      "supported for rotations about the z-axis, but surface {} is not "
-      "perpendicular to the z-axis.",
-      surf2.id_));
-  }
-
   // Make sure both surfaces intersect the origin
   if (std::abs(surf1.evaluate({0, 0, 0})) > FP_COINCIDENT) {
     throw std::invalid_argument(fmt::format(
@@ -231,15 +207,8 @@ RotationalPeriodicBC::RotationalPeriodicBC(int i_surf, int j_surf)
       surf2.id_));
   }
 
-  // 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(norm1.y, norm1.x);
-  double theta2 = std::atan2(norm2.y, norm2.x) + PI;
-  angle_ = theta2 - theta1;
+  angle_ = compute_periodic_rotation(norm1[axis_2_idx_], norm1[axis_1_idx_],
+    norm2[axis_2_idx_], norm2[axis_1_idx_]);
 
   // Warn the user if the angle does not evenly divide a circle
   double rem = std::abs(std::remainder((2 * PI / angle_), 1.0));
@@ -251,6 +220,20 @@ RotationalPeriodicBC::RotationalPeriodicBC(int i_surf, int j_surf)
   }
 }
 
+float RotationalPeriodicBC::compute_periodic_rotation(
+  float rise_1, float run_1, float rise_2, float 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
 {
@@ -278,10 +261,10 @@ void RotationalPeriodicBC::handle_particle(
   Direction u = p.u();
   double cos_theta = std::cos(theta);
   double sin_theta = std::sin(theta);
-  Position new_r = {
-    cos_theta * r.x - sin_theta * r.y, sin_theta * r.x + cos_theta * r.y, r.z};
-  Direction new_u = {
-    cos_theta * u.x - sin_theta * u.y, sin_theta * u.x + cos_theta * u.y, u.z};
+  Position new_r = {cos_theta * r[axis_1_idx_] - sin_theta * r[axis_2_idx_],
+    sin_theta * r[axis_1_idx_] + cos_theta * r[axis_2_idx_], r[zero_axis_idx]};
+  Direction new_u = {cos_theta * u[axis_1_idx_] - sin_theta * u[axis_2_idx_],
+    sin_theta * u[axis_1_idx_] + cos_theta * u[axis_2_idx_], u[zero_axis_idx]};
 
   // Handle the effects of the surface albedo on the particle's weight.
   BoundaryCondition::handle_albedo(p, surf);

tests/regression_tests/periodic_cyls/test.py

@@ -0,0 +1,91 @@
+import openmc
+import numpy as np
+import pytest
+from openmc.utility_funcs import change_directory
+from tests.testing_harness import PyAPITestHarness
+
+
+@pytest.fixture
+def xcyl_model():
+    model = openmc.Model()
+    # Define materials
+    fuel = openmc.Material()
+    fuel.add_nuclide('U235', 0.2)
+    fuel.add_nuclide('U238', 0.8)
+    fuel.set_density('g/cc', 19.1)
+    model.materials = openmc.Materials([fuel])
+
+    # Define geometry
+    # finite cylinder
+    x_min = openmc.XPlane(x0=0.0, boundary_type='reflective')
+    x_max = openmc.XPlane(x0=20.0, boundary_type='reflective')
+    x_cyl = openmc.XCylinder(r=20.0,boundary_type='vacuum')
+    # slice cylinder for periodic BC
+    periodic_bounding_yplane = openmc.YPlane(y0=0, boundary_type='periodic')
+    periodic_bounding_plane = openmc.Plane(
+        a=0.0, b=-np.sqrt(3) / 3, c=1, boundary_type='periodic',
+    )
+    sixth_cyl_cell = openmc.Cell(1, fill=fuel, region = 
+        +x_min &- x_max & -x_cyl & +periodic_bounding_yplane & +periodic_bounding_plane)
+    periodic_bounding_yplane.periodic_surface = periodic_bounding_plane
+    periodic_bounding_plane.periodic_surface = periodic_bounding_yplane
+
+    model.geometry = openmc.Geometry([sixth_cyl_cell])
+
+
+    # Define settings
+    model.settings.particles = 1000
+    model.settings.batches = 4
+    model.settings.inactive = 0
+    model.settings.source = openmc.IndependentSource(space=openmc.stats.Box(
+        (0, 0, 0), (20, 20, 20))
+    )
+    return model
+
+@pytest.fixture
+def ycyl_model():
+    model = openmc.Model()
+    # Define materials
+    fuel = openmc.Material()
+    fuel.add_nuclide('U235', 0.2)
+    fuel.add_nuclide('U238', 0.8)
+    fuel.set_density('g/cc', 19.1)
+    model.materials = openmc.Materials([fuel])
+
+    # Define geometry
+    # finite cylinder
+    y_min = openmc.YPlane(y0=0.0, boundary_type='reflective')
+    y_max = openmc.YPlane(y0=20.0, boundary_type='reflective')
+    y_cyl = openmc.YCylinder(r=20.0,boundary_type='vacuum')
+    # slice cylinder for periodic BC
+    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',
+    )
+    sixth_cyl_cell = openmc.Cell(1, fill=fuel, region = 
+        +y_min &- y_max & -y_cyl & +periodic_bounding_xplane & +periodic_bounding_plane)
+    periodic_bounding_xplane.periodic_surface = periodic_bounding_plane
+    periodic_bounding_plane.periodic_surface = periodic_bounding_xplane
+    model.geometry = openmc.Geometry([sixth_cyl_cell])
+
+
+    # Define settings
+    model.settings.particles = 1000
+    model.settings.batches = 4
+    model.settings.inactive = 0
+    model.settings.source = openmc.IndependentSource(space=openmc.stats.Box(
+        (0, 0, 0), (20, 20, 20))
+    )
+    return model
+
+def test_xcyl(xcyl_model):
+    with change_directory("xcyl_model"):
+        openmc.reset_auto_ids()
+        harness = PyAPITestHarness('statepoint.4.h5', xcyl_model)
+        harness.main()
+
+def test_ycyl(ycyl_model):
+    with change_directory("ycyl_model"):
+        openmc.reset_auto_ids()
+        harness = PyAPITestHarness('statepoint.4.h5', ycyl_model)
+        harness.main()

tests/regression_tests/periodic_cyls/xcyl_model/inputs_true.dat

@@ -0,0 +1,29 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1" depletable="true">
+      <density value="19.1" units="g/cc"/>
+      <nuclide name="U235" ao="0.2"/>
+      <nuclide name="U238" ao="0.8"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="1 -2 -3 4 5" universe="1"/>
+    <surface id="1" type="x-plane" boundary="reflective" coeffs="0.0"/>
+    <surface id="2" type="x-plane" boundary="reflective" coeffs="20.0"/>
+    <surface id="3" type="x-cylinder" boundary="vacuum" coeffs="0.0 0.0 20.0"/>
+    <surface id="4" type="y-plane" boundary="periodic" coeffs="0" periodic_surface_id="5"/>
+    <surface id="5" type="plane" boundary="periodic" coeffs="0.0 -0.5773502691896257 1 0.0" periodic_surface_id="4"/>
+  </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 20 20 20</parameters>
+      </space>
+    </source>
+  </settings>
+  </model>

tests/regression_tests/periodic_cyls/ycyl_model/inputs_true.dat

@@ -0,0 +1,29 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1" depletable="true">
+      <density value="19.1" units="g/cc"/>
+      <nuclide name="U235" ao="0.2"/>
+      <nuclide name="U238" ao="0.8"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="1 -2 -3 4 5" universe="1"/>
+    <surface id="1" type="y-plane" boundary="reflective" coeffs="0.0"/>
+    <surface id="2" type="y-plane" boundary="reflective" coeffs="20.0"/>
+    <surface id="3" type="y-cylinder" boundary="vacuum" coeffs="0.0 0.0 20.0"/>
+    <surface id="4" type="x-plane" boundary="periodic" coeffs="0" periodic_surface_id="5"/>
+    <surface id="5" type="plane" boundary="periodic" coeffs="-0.5773502691896257 0.0 1 0.0" periodic_surface_id="4"/>
+  </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 20 20 20</parameters>
+      </space>
+    </source>
+  </settings>
+  </model>