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

include/openmc/surface.h

@@ -2,6 +2,7 @@
 #define OPENMC_SURFACE_H
 
 #include <limits> // For numeric_limits
+#include <set>
 #include <string>
 #include <unordered_map>
 
@@ -378,7 +379,39 @@ class SurfaceZTorus : public Surface {
 // Non-member functions
 //==============================================================================
 
-void read_surfaces(pugi::xml_node node);
+//! Read surface definitions from XML and populate the global surfaces vector.
+//!
+//! This function parses surface elements from the XML input, creates the
+//! appropriate surface objects, and identifies periodic surfaces along with
+//! their albedo values and sense information.
+//!
+//! \param node XML node containing surface definitions
+//! \param[out] periodic_pairs Set of surface ID pairs representing periodic
+//!   boundary conditions
+//! \param[out] albedo_map Map of surface IDs to albedo values for periodic
+//!   surfaces
+//! \param[out] periodic_sense_map Map of surface IDs to their sense values
+//!   (used to determine orientation for periodic BCs)
+void read_surfaces(pugi::xml_node node,
+  std::set<std::pair<int, int>>& periodic_pairs,
+  std::unordered_map<int, double>& albedo_map,
+  std::unordered_map<int, int>& periodic_sense_map);
+
+//! Resolve periodic surface pairs and assign boundary conditions.
+//!
+//! This function completes the setup of periodic boundary conditions by
+//! resolving unpaired periodic surfaces, determining whether each pair
+//! represents translational or rotational periodicity based on surface
+//! normals, and assigning the appropriate boundary condition objects.
+//!
+//! \param[inout] periodic_pairs Set of surface ID pairs representing periodic
+//!   boundary conditions; unpaired entries are resolved
+//! \param albedo_map Map of surface IDs to albedo values for periodic surfaces
+//! \param periodic_sense_map Map of surface IDs to their sense values (used to
+//!   determine orientation for periodic BCs)
+void prepare_boundary_conditions(std::set<std::pair<int, int>>& periodic_pairs,
+  std::unordered_map<int, double>& albedo_map,
+  std::unordered_map<int, int>& periodic_sense_map);
 
 void free_memory_surfaces();
 

src/boundary_condition.cpp

@@ -160,7 +160,7 @@ void TranslationalPeriodicBC::handle_particle(
 
 RotationalPeriodicBC::RotationalPeriodicBC(
   int i_surf, int j_surf, PeriodicAxis axis)
-  : PeriodicBC(i_surf, j_surf)
+  : PeriodicBC(std::abs(i_surf) - 1, std::abs(j_surf) - 1)
 {
   Surface& surf1 {*model::surfaces[i_surf_]};
   Surface& surf2 {*model::surfaces[j_surf_]};
@@ -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,10 +187,16 @@ 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;
+
+  auto i_sign = std::copysign(1, i_surf);
+  auto j_sign = -std::copysign(1, j_surf);
+
   // Compute the surface normal vectors and make sure they are perpendicular
   // to the correct axis
-  Direction norm1 = surf1.normal({0, 0, 0});
-  Direction norm2 = surf2.normal({0, 0, 0});
+  Direction norm1 = i_sign * surf1.normal({0, 0, 0});
+  Direction norm2 = j_sign * surf2.normal({0, 0, 0});
   // Make sure both surfaces intersect the origin
   if (std::abs(surf1.evaluate({0, 0, 0})) > FP_COINCIDENT) {
     throw std::invalid_argument(fmt::format(
@@ -212,8 +213,15 @@ 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_ = (i_sign * j_sign > 0.0);
 
   // 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 +233,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/geometry_aux.cpp

@@ -55,8 +55,13 @@ void read_geometry_xml()
 void read_geometry_xml(pugi::xml_node root)
 {
   // Read surfaces, cells, lattice
-  read_surfaces(root);
+  std::set<std::pair<int, int>> periodic_pairs;
+  std::unordered_map<int, double> albedo_map;
+  std::unordered_map<int, int> periodic_sense_map;
+
+  read_surfaces(root, periodic_pairs, albedo_map, periodic_sense_map);
   read_cells(root);
+  prepare_boundary_conditions(periodic_pairs, albedo_map, periodic_sense_map);
   read_lattices(root);
 
   // Check to make sure a boundary condition was applied to at least one

src/particle.cpp

@@ -744,9 +744,7 @@ void Particle::cross_periodic_bc(
   if (!neighbor_list_find_cell(*this)) {
     mark_as_lost("Couldn't find particle after hitting periodic "
                  "boundary on surface " +
-                 std::to_string(surf.id_) +
-                 ". The normal vector "
-                 "of one periodic surface may need to be reversed.");
+                 std::to_string(surf.id_) + ".");
     return;
   }
 

src/surface.cpp

@@ -9,6 +9,7 @@
 #include <fmt/core.h>
 
 #include "openmc/array.h"
+#include "openmc/cell.h"
 #include "openmc/container_util.h"
 #include "openmc/error.h"
 #include "openmc/external/quartic_solver.h"
@@ -1169,7 +1170,10 @@ Direction SurfaceZTorus::normal(Position r) const
 
 //==============================================================================
 
-void read_surfaces(pugi::xml_node node)
+void read_surfaces(pugi::xml_node node,
+  std::set<std::pair<int, int>>& periodic_pairs,
+  std::unordered_map<int, double>& albedo_map,
+  std::unordered_map<int, int>& periodic_sense_map)
 {
   // Count the number of surfaces
   int n_surfaces = 0;
@@ -1180,8 +1184,6 @@ void read_surfaces(pugi::xml_node node)
   // Loop over XML surface elements and populate the array.  Keep track of
   // periodic surfaces and their albedos.
   model::surfaces.reserve(n_surfaces);
-  std::set<std::pair<int, int>> periodic_pairs;
-  std::unordered_map<int, double> albedo_map;
   {
     pugi::xml_node surf_node;
     int i_surf;
@@ -1244,6 +1246,7 @@ void read_surfaces(pugi::xml_node node)
       if (check_for_node(surf_node, "boundary")) {
         std::string surf_bc = get_node_value(surf_node, "boundary", true, true);
         if (surf_bc == "periodic") {
+          periodic_sense_map[model::surfaces.back()->id_] = 0;
           // Check for surface albedo. Skip sanity check as it is already done
           // in the Surface class's constructor.
           if (check_for_node(surf_node, "albedo")) {
@@ -1275,6 +1278,28 @@ void read_surfaces(pugi::xml_node node)
         fmt::format("Two or more surfaces use the same unique ID: {}", id));
     }
   }
+}
+
+void prepare_boundary_conditions(std::set<std::pair<int, int>>& periodic_pairs,
+  std::unordered_map<int, double>& albedo_map,
+  std::unordered_map<int, int>& periodic_sense_map)
+{
+  // Fill the senses map for periodic surfaces
+  auto n_periodic = periodic_sense_map.size();
+  for (const auto& cell : model::cells) {
+    if (n_periodic == 0)
+      break; // Early exit once all periodic surfaces found
+
+    for (auto s : cell->surfaces()) {
+      auto surf_idx = std::abs(s) - 1;
+      auto id = model::surfaces[surf_idx]->id_;
+
+      if (periodic_sense_map.count(id)) {
+        periodic_sense_map[id] = std::copysign(1, s);
+        --n_periodic;
+      }
+    }
+  }
 
   // Resolve unpaired periodic surfaces.  A lambda function is used with
   // std::find_if to identify the unpaired surfaces.
@@ -1370,8 +1395,12 @@ void read_surfaces(pugi::xml_node node)
           "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);
+      auto i_sign = periodic_sense_map[periodic_pair.first];
+      auto j_sign = periodic_sense_map[periodic_pair.second];
+      surf1.bc_ = make_unique<RotationalPeriodicBC>(
+        i_sign * (i_surf + 1), j_sign * (j_surf + 1), axis);
+      surf2.bc_ = make_unique<RotationalPeriodicBC>(
+        j_sign * (j_surf + 1), i_sign * (i_surf + 1), 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