add error checking to determine which periodic axis is specified by the user

Author: lewisgross1296

include/openmc/boundary_condition.h

@@ -144,7 +144,7 @@ class TranslationalPeriodicBC : public PeriodicBC {
 class RotationalPeriodicBC : public PeriodicBC {
 public:
   enum PeriodicAxis { x, y, z };
-  RotationalPeriodicBC(int i_surf, int j_surf, PeriodicAxis axis = 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;

src/surface.cpp

@@ -1334,8 +1334,37 @@ 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 (norm1.x != 0.0 && norm1.y != 0.0 && norm1.z != 0.0) {
+        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 (norm2.x != 0.0 && norm2.y != 0.0 && norm2.z != 0.0) {
+        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 (norm1.x == 0.0 && norm2.x == 0.0) {
+        axis = RotationalPeriodicBC::PeriodicAxis::x;
+      } else if (norm1.y == 0.0 && norm2.y == 0.0) {
+        axis = RotationalPeriodicBC::PeriodicAxis::y;
+      } else if (norm1.z == 0.0 && norm2.z == 0.0) {
+        axis = RotationalPeriodicBC::PeriodicAxis::z;
+      } else {
+        fatal_error(
+          fmt::format("There is no component which is 0.0 in both normals. 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/xcyl_model/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+1.082283E+00 6.676373E-02

tests/regression_tests/periodic_cyls/ycyl_model/inputs_true.dat

@@ -1,19 +1,19 @@
 <?xml version='1.0' encoding='utf-8'?>
 <model>
   <materials>
-    <material id="1" depletable="true">
+    <material id="2" 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"/>
+    <cell id="1" material="2" region="6 -7 -8 9 10" universe="2"/>
+    <surface id="6" type="y-plane" boundary="reflective" coeffs="0.0"/>
+    <surface id="7" type="y-plane" boundary="reflective" coeffs="20.0"/>
+    <surface id="8" type="y-cylinder" boundary="vacuum" coeffs="0.0 0.0 20.0"/>
+    <surface id="9" type="x-plane" boundary="periodic" coeffs="0" periodic_surface_id="10"/>
+    <surface id="10" type="plane" boundary="periodic" coeffs="-0.5773502691896257 0.0 1 0.0" periodic_surface_id="9"/>
   </geometry>
   <settings>
     <run_mode>eigenvalue</run_mode>