Add periodic boundary conditions

docs/source/usersguide/input.rst

@@ -898,7 +898,7 @@ Each ``<surface>`` element can have the following attributes or sub-elements:
 
   :boundary:
     The boundary condition for the surface. This can be "transmission",
-    "vacuum", or "reflective".
+    "vacuum", "reflective", or "periodic".
 
     *Default*: "transmission"
 

src/geometry.F90

@@ -475,6 +475,68 @@ contains
              &// trim(to_str(surf%id)))
       end if
       return
+    elseif (surf % bc == BC_PERIODIC .and. run_mode /= MODE_PLOTTING) then
+      ! =======================================================================
+      ! PERIODIC BOUNDARY
+
+      ! Do not handle periodic boundary conditions on lower universes
+      if (p % n_coord /= 1) then
+        call handle_lost_particle(p, "Cannot period particle " &
+             // trim(to_str(p % id)) // " off surface in a lower universe.")
+        return
+      end if
+
+      ! Score surface currents since reflection causes the direction of the
+      ! particle to change -- artificially move the particle slightly back in
+      ! case the surface crossing in coincident with a mesh boundary
+
+      if (active_current_tallies % size() > 0) then
+        p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw
+        call score_surface_current(p)
+        p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw
+      end if
+
+      select type (surf)
+      type is (SurfaceXPlane)
+        select type (opposite => surfaces(surf % opposite) % obj)
+        type is (SurfaceXPlane)
+          p % coord(1) % xyz(1) = opposite % x0
+        end select
+
+      type is (SurfaceYPlane)
+        select type (opposite => surfaces(surf % opposite) % obj)
+        type is (SurfaceYPlane)
+          p % coord(1) % xyz(2) = opposite % y0
+        end select
+
+      type is (SurfaceZPlane)
+        select type (opposite => surfaces(surf % opposite) % obj)
+        type is (SurfaceZPlane)
+          p % coord(1) % xyz(3) = opposite % z0
+        end select
+      end select
+
+      ! Reassign particle's surface
+      p % surface = sign(surfaces(surf % opposite) % obj % id, p % surface)
+
+      ! Figure out what cell particle is in now
+      p % n_coord = 1
+      call find_cell(p, found)
+      if (.not. found) then
+        call handle_lost_particle(p, "Couldn't find particle after hitting &
+             &periodic boundary on surface " // trim(to_str(surf%id)) // ".")
+        return
+      end if
+
+      ! Set previous coordinate going slightly past surface crossing
+      p % last_xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw
+
+      ! Diagnostic message
+      if (verbosity >= 10 .or. trace) then
+        call write_message("    Hit periodic boundary on surface " &
+             // trim(to_str(surf%id)))
+      end if
+      return
     end if
 
     ! ==========================================================================

src/initialize.F90

@@ -580,6 +580,20 @@ contains
     class(Lattice),    pointer :: lat => null()
     type(TallyObject), pointer :: t => null()
 
+    ! Adjust opposite surfaces for periodic boundaries
+    do i = 1, size(surfaces)
+      associate (surf => surfaces(i) % obj)
+        if (surf % bc == BC_PERIODIC) then
+          if (surface_dict % has_key(surf % opposite)) then
+            surf % opposite = surface_dict % get_key(surf % opposite)
+          else
+            call fatal_error("Could not find opposite surface " // &
+                 trim(to_str(surf % opposite)) // ".")
+          end if
+        end if
+      end associate
+    end do
+
     do i = 1, n_cells
       ! =======================================================================
       ! ADJUST REGION SPECIFICATION FOR EACH CELL

src/input_xml.F90

@@ -1109,6 +1109,8 @@ contains
     integer :: universe_num
     integer :: n_cells_in_univ
     integer :: coeffs_reqd
+    integer :: i_xmin, i_xmax, i_ymin, i_ymax, i_zmin, i_zmax
+    real(8) :: xmin, xmax, ymin, ymax, zmin, zmax
     integer, allocatable :: temp_int_array(:)
     real(8) :: phi, theta, psi
     real(8), allocatable :: coeffs(:)
@@ -1387,6 +1389,13 @@ contains
       call fatal_error("No surfaces found in geometry.xml!")
     end if
 
+    xmin = INFINITY
+    xmax = -INFINITY
+    ymin = INFINITY
+    ymax = -INFINITY
+    zmin = INFINITY
+    zmax = -INFINITY
+
     ! Allocate cells array
     allocate(surfaces(n_surfaces))
 
@@ -1478,10 +1487,28 @@ contains
       select type(s)
       type is (SurfaceXPlane)
         s%x0 = coeffs(1)
+
+        ! Determine outer surfaces
+        xmin = min(xmin, s % x0)
+        xmax = max(xmax, s % x0)
+        if (xmin == s % x0) i_xmin = i
+        if (xmax == s % x0) i_xmax = i
       type is (SurfaceYPlane)
         s%y0 = coeffs(1)
+
+        ! Determine outer surfaces
+        ymin = min(ymin, s % y0)
+        ymax = max(ymax, s % y0)
+        if (ymin == s % y0) i_ymin = i
+        if (ymax == s % y0) i_ymax = i
       type is (SurfaceZPlane)
         s%z0 = coeffs(1)
+
+        ! Determine outer surfaces
+        zmin = min(zmin, s % z0)
+        zmax = max(zmax, s % z0)
+        if (zmin == s % z0) i_zmin = i
+        if (zmax == s % z0) i_zmax = i
       type is (SurfacePlane)
         s%A = coeffs(1)
         s%B = coeffs(2)
@@ -1548,11 +1575,13 @@ contains
       case ('reflective', 'reflect', 'reflecting')
         s%bc = BC_REFLECT
         boundary_exists = .true.
+      case ('periodic')
+        s%bc = BC_PERIODIC
+        boundary_exists = .true.
       case default
         call fatal_error("Unknown boundary condition '" // trim(word) // &
              &"' specified on surface " // trim(to_str(s%id)))
       end select
-
       ! Add surface to dictionary
       call surface_dict % add_key(s%id, i)
     end do
@@ -1563,6 +1592,47 @@ contains
       call fatal_error("No boundary conditions were applied to any surfaces!")
     end if
 
+    ! Determine opposite side for periodic boundaries
+    do i = 1, size(surfaces)
+      if (surfaces(i) % obj % bc == BC_PERIODIC) then
+        select type (surf => surfaces(i) % obj)
+        type is (SurfaceXPlane)
+          if (i == i_xmin) then
+            surf % opposite = i_xmax
+          elseif (i == i_xmax) then
+            surf % opposite = i_xmin
+          else
+            call fatal_error("Periodic boundary condition applied to &
+                 &interior surface.")
+          end if
+
+        type is (SurfaceYPlane)
+          if (i == i_ymin) then
+            surf % opposite = i_ymax
+          elseif (i == i_ymax) then
+            surf % opposite = i_ymin
+          else
+            call fatal_error("Periodic boundary condition applied to &
+                 &interior surface.")
+          end if
+
+        type is (SurfaceZPlane)
+          if (i == i_zmin) then
+            surf % opposite = i_zmax
+          elseif (i == i_zmax) then
+            surf % opposite = i_zmin
+          else
+            call fatal_error("Periodic boundary condition applied to &
+                 &interior surface.")
+          end if
+
+        class default
+          call fatal_error("Periodic boundary condition applied to &
+               &non-planar surface.")
+        end select
+      end if
+    end do
+
     ! ==========================================================================
     ! READ LATTICES FROM GEOMETRY.XML
 

src/relaxng/geometry.rnc

@@ -21,8 +21,9 @@ element geometry {
     (element type { xsd:string { maxLength = "15" } } |
      attribute type { xsd:string { maxLength = "15" } }) &
     (element coeffs { list { xsd:double+ } } | attribute coeffs { list { xsd:double+ } }) &
-    (element boundary { ( "transmit" | "reflective" | "vacuum" ) } |
-     attribute boundary { ( "transmit" | "reflective" | "vacuum" ) })?
+    (element boundary { ( "transmit" | "reflective" | "vacuum" | "periodic" ) } |
+     attribute boundary { ( "transmit" | "reflective" | "vacuum" | "periodic" ) })? &
+    (element opposite { xsd:int } | attribute opposite { xsd:int })?
   }*
 
   & element lattice {

src/relaxng/geometry.rng

@@ -173,6 +173,7 @@
                   <value>transmit</value>
                   <value>reflective</value>
                   <value>vacuum</value>
+                  <value>periodic</value>
                 </choice>
               </element>
               <attribute name="boundary">
@@ -180,10 +181,21 @@
                   <value>transmit</value>
                   <value>reflective</value>
                   <value>vacuum</value>
+                  <value>periodic</value>
                 </choice>
               </attribute>
             </choice>
           </optional>
+          <optional>
+            <choice>
+              <element name="opposite">
+                <data type="int"/>
+              </element>
+              <attribute name="opposite">
+                <data type="int"/>
+              </attribute>
+            </choice>
+          </optional>
         </interleave>
       </element>
     </zeroOrMore>

src/surface_header.F90

@@ -15,6 +15,7 @@ module surface_header
          neighbor_pos(:), &           ! List of cells on positive side
          neighbor_neg(:)              ! List of cells on negative side
     integer :: bc                     ! Boundary condition
+    integer :: opposite               ! Opposite surface for periodic boundary
     character(len=104) :: name = ""    ! User-defined name
   contains
     procedure :: sense

tests/test_periodic/geometry.xml

@@ -0,0 +1,12 @@
+<?xml version='1.0' encoding='utf-8'?>
+<geometry>
+  <surface id="1" type="x-plane" coeffs="-5." boundary="periodic"  />
+  <surface id="2" type="x-plane" coeffs="5." boundary="periodic"  />
+  <surface id="3" type="y-plane" coeffs="-5." boundary="periodic"  />
+  <surface id="4" type="y-plane" coeffs="5." boundary="periodic"  />
+  <surface id="5" type="z-plane" coeffs="-5." boundary="reflective"  />
+  <surface id="6" type="z-plane" coeffs="5." boundary="reflective"  />
+  <surface id="7" type="z-cylinder" coeffs="-2.5 2.5 2.0" />
+  <cell id="1" material="1" region="1 -2 3 -4 5 -6 7" />
+  <cell id="2" material="2" region="5 -6 -7" />
+</geometry>

tests/test_periodic/materials.xml

@@ -0,0 +1,13 @@
+<?xml version='1.0' encoding='utf-8'?>
+<materials>
+  <material id="1">
+    <density units="g/cc" value="1.0" />
+    <nuclide name="H-1" xs="71c" ao="2.0" />
+    <nuclide name="O-16" xs="71c" ao="2.0" />
+    <sab name="HH2O" xs="71t" />
+  </material>
+  <material id="2">
+    <density units="g/cc" value="4.5" />
+    <nuclide name="U-235" xs="71c" ao="1.0" />
+  </material>
+</materials>

tests/test_periodic/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+1.542742E+00 4.410461E-02

tests/test_periodic/settings.xml

@@ -0,0 +1,13 @@
+<?xml version='1.0' encoding='utf-8'?>
+<settings>
+  <eigenvalue>
+    <particles>1000</particles>
+    <batches>4</batches>
+    <inactive>0</inactive>
+  </eigenvalue>
+  <source strength="1.0">
+    <space type="box">
+      <parameters>-5. -5. -5. 5. 5. 5.</parameters>
+    </space>
+  </source>
+</settings>

tests/test_periodic/tallies.xml

@@ -0,0 +1,14 @@
+<?xml version="1.0"?>
+<tallies>
+  <mesh id="1">
+    <type>regular</type>
+    <lower_left>-200. -1e50</lower_left>
+    <upper_right>200. 1e50</upper_right>
+    <dimension>50 1</dimension>
+  </mesh>
+  <tally id="1">
+    <estimator>collision</estimator>
+    <filter type="mesh" bins="1" />
+    <scores>fission</scores>
+  </tally>
+</tallies>

tests/test_periodic/test_periodic.py

@@ -0,0 +1,11 @@
+#!/usr/bin/env python
+
+import os
+import sys
+sys.path.insert(0, os.pardir)
+from testing_harness import TestHarness
+
+
+if __name__ == '__main__':
+    harness = TestHarness('statepoint.4.h5')
+    harness.main()

tests/testing_harness.py

@@ -133,9 +133,10 @@ class TestHarness(object):
 
     def _cleanup(self):
         """Delete statepoints, tally, and test files."""
-        output = glob.glob(os.path.join(os.getcwd(), 'statepoint.*.*'))
+        output = glob.glob(os.path.join(os.getcwd(), 'statepoint.*.h5'))
         output.append(os.path.join(os.getcwd(), 'tallies.out'))
         output.append(os.path.join(os.getcwd(), 'results_test.dat'))
+        output.append(os.path.join(os.getcwd(), 'summary.h5'))
         for f in output:
             if os.path.exists(f):
                 os.remove(f)