Add support for periodic boundary conditions on general planes

openmc/surface.py

@@ -2110,6 +2110,11 @@ def get_hexagonal_prism(edge_length=1., orientation='y',
         prism = -right & +left & -upper_right & -upper_left & \
                 +lower_right & +lower_left
 
+        if boundary_type == 'periodic':
+            right.periodic_surface = left
+            upper_right.periodic_surface = lower_left
+            lower_right.periodic_surface = upper_left
+
     elif orientation == 'x':
         top = YPlane(y0=sqrt(3.)/2*l, boundary_type=boundary_type)
         bottom = YPlane(y0=-sqrt(3.)/2*l, boundary_type=boundary_type)
@@ -2129,6 +2134,11 @@ def get_hexagonal_prism(edge_length=1., orientation='y',
         prism = -top & +bottom & -upper_right & +lower_right & \
                             +lower_left & -upper_left
 
+        if boundary_type == 'periodic':
+            top.periodic_surface = bottom
+            upper_right.periodic_surface = lower_left
+            lower_right.periodic_surface = upper_left
+
     # Handle rounded corners if given
     if corner_radius > 0.:
         c = sqrt(3.)/2

src/geometry.F90

@@ -392,6 +392,7 @@ contains
     real(8) :: v         ! y-component of direction
     real(8) :: w         ! z-component of direction
     real(8) :: norm      ! "norm" of surface normal
+    real(8) :: d         ! distance between point and plane
     real(8) :: xyz(3)    ! Saved global coordinate
     integer :: i_surface ! index in surfaces
     logical :: found     ! particle found in universe?
@@ -532,6 +533,20 @@ contains
         type is (SurfaceZPlane)
           p % coord(1) % xyz(3) = opposite % z0
         end select
+
+      type is (SurfacePlane)
+        select type (opposite => surfaces(surf % i_periodic) % obj)
+        type is (SurfacePlane)
+          ! Get surface normal for opposite plane
+          xyz(:) = opposite % normal(p % coord(1) % xyz)
+
+          ! Determine distance to plane
+          norm = xyz(1)*xyz(1) + xyz(2)*xyz(2) + xyz(3)*xyz(3)
+          d = opposite % evaluate(p % coord(1) % xyz) / norm
+
+          ! Move particle along normal vector based on distance
+          p % coord(1) % xyz(:) = p % coord(1) % xyz(:) - d*xyz
+        end select
       end select
 
       ! Reassign particle's surface

src/input_xml.F90

@@ -1658,6 +1658,15 @@ contains
             surf % i_periodic = surface_dict % get_key(surf % i_periodic)
           end if
 
+        type is (SurfacePlane)
+          if (surf % i_periodic == NONE) then
+            call fatal_error("No matching periodic surface specified for &
+                 &periodic boundary condition on surface " // &
+                 trim(to_str(surf % id)) // ".")
+          else
+            surf % i_periodic = surface_dict % get_key(surf % i_periodic)
+          end if
+
         class default
           call fatal_error("Periodic boundary condition applied to &
                &non-planar surface.")