Generalize RotationalPeriodicBC for X-, Y-, or Z-axis (#3591)

Co-authored-by: GuySten <62616591+GuySten@users.noreply.github.com>

Author: lewisgross1296

docs/source/io_formats/geometry.rst

@@ -38,11 +38,9 @@ Each ``<surface>`` element can have the following attributes or sub-elements:
 
   :boundary:
      The boundary condition for the surface. This can be "transmission",
-     "vacuum", "reflective", or "periodic". 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.
+     "vacuum", "reflective", or "periodic". Specify which planes are
+     periodic and the code will automatically identify which planes are
+     paired together.
 
     *Default*: "transmission"
 

docs/source/usersguide/geometry.rst

@@ -192,7 +192,7 @@ Otherwise it is necessary to specify pairs explicitly using the
 Both rotational and translational periodic boundary conditions are specified in
 the same fashion. If both planes have the same normal vector, a translational
 periodicity is assumed; rotational periodicity is assumed otherwise. Currently,
-only rotations about the :math:`z`-axis are supported.
+rotations must be about the :math:`x`-, :math:`y`-, or :math:`z`-axis.
 
 For a rotational periodic BC, the normal vectors of each surface must point
 inwards---towards the valid geometry. For example, a :class:`XPlane` and

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);
+  double compute_periodic_rotation(
+    double rise_1, double run_1, double rise_2, double 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 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
@@ -822,8 +821,7 @@ class XPlane(PlaneMixin, Surface):
     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. Only axis-aligned periodicity is
-        supported, i.e., x-planes can only be paired with x-planes.
+        freely pass through the surface.
     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
@@ -887,8 +885,7 @@ class YPlane(PlaneMixin, Surface):
     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. Only axis-aligned periodicity is
-        supported, i.e., y-planes can only be paired with y-planes.
+        freely pass through the surface.
     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
@@ -952,8 +949,7 @@ class ZPlane(PlaneMixin, Surface):
     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. Only axis-aligned periodicity is
-        supported, i.e., z-planes can only be paired with z-planes.
+        freely pass through the surface.
     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,44 @@ 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:
+    // 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
+    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 +212,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 +225,20 @@ RotationalPeriodicBC::RotationalPeriodicBC(int i_surf, int j_surf)
   }
 }
 
+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
 {
@@ -278,10 +266,16 @@ 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;
+  new_r[zero_axis_idx_] = r[zero_axis_idx_];
+  new_r[axis_1_idx_] = cos_theta * r[axis_1_idx_] - sin_theta * r[axis_2_idx_];
+  new_r[axis_2_idx_] = sin_theta * r[axis_1_idx_] + cos_theta * r[axis_2_idx_];
+
+  Direction new_u;
+  new_u[zero_axis_idx_] = u[zero_axis_idx_];
+  new_u[axis_1_idx_] = cos_theta * u[axis_1_idx_] - sin_theta * u[axis_2_idx_];
+  new_u[axis_2_idx_] = sin_theta * u[axis_1_idx_] + cos_theta * u[axis_2_idx_];
 
   // Handle the effects of the surface albedo on the particle's weight.
   BoundaryCondition::handle_albedo(p, surf);

src/surface.cpp

@@ -1334,8 +1334,44 @@ void read_surfaces(pugi::xml_node node)
       surf1.bc_ = make_unique<TranslationalPeriodicBC>(i_surf, j_surf);
       surf2.bc_ = make_unique<TranslationalPeriodicBC>(i_surf, j_surf);
     } else {
-      surf1.bc_ = make_unique<RotationalPeriodicBC>(i_surf, j_surf);
-      surf2.bc_ = make_unique<RotationalPeriodicBC>(i_surf, j_surf);
+      // check that both normals have at least one 0 component
+      if (std::abs(norm1.x) > FP_PRECISION &&
+          std::abs(norm1.y) > FP_PRECISION &&
+          std::abs(norm1.z) > FP_PRECISION) {
+        fatal_error(fmt::format(
+          "The normal ({}) of the periodic surface ({}) does not contain any "
+          "component with a zero value. A RotationalPeriodicBC requires one "
+          "component which is zero for both plane normals.",
+          norm1, i_surf));
+      }
+      if (std::abs(norm2.x) > FP_PRECISION &&
+          std::abs(norm2.y) > FP_PRECISION &&
+          std::abs(norm2.z) > FP_PRECISION) {
+        fatal_error(fmt::format(
+          "The normal ({}) of the periodic surface ({}) does not contain any "
+          "component with a zero value. A RotationalPeriodicBC requires one "
+          "component which is zero for both plane normals.",
+          norm2, j_surf));
+      }
+      // find common zero component, which indicates the periodic axis
+      RotationalPeriodicBC::PeriodicAxis axis;
+      if (std::abs(norm1.x) <= FP_PRECISION &&
+          std::abs(norm2.x) <= FP_PRECISION) {
+        axis = RotationalPeriodicBC::PeriodicAxis::x;
+      } else if (std::abs(norm1.y) <= FP_PRECISION &&
+                 std::abs(norm2.y) <= FP_PRECISION) {
+        axis = RotationalPeriodicBC::PeriodicAxis::y;
+      } else if (std::abs(norm1.z) <= FP_PRECISION &&
+                 std::abs(norm2.z) <= FP_PRECISION) {
+        axis = RotationalPeriodicBC::PeriodicAxis::z;
+      } else {
+        fatal_error(fmt::format(
+          "There is no component which is 0.0 in both normal vectors. This "
+          "indicates that the two planes are not periodic about the X, Y, or Z "
+          "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);
     }
 
     // If albedo data is present in albedo map, set the boundary albedo.

tests/regression_tests/periodic_cyls/__init__.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()