Merge pull request #236 from mit-crpg/source_interval

More control for source bank writing

docs/source/usersguide/input.rst

@@ -348,8 +348,10 @@ attributes/sub-elements:
 
 The ``<state_point>`` element indicates at what batches a state point file
 should be written. A state point file can be used to restart a run or to get
-tally results at any batch. This element has the following
-attributes/sub-elements:
+tally results at any batch. The default behavior when using this tag is to 
+write out the source bank in the state_point file. This behavior can be 
+customized by using the ``<source_point>`` element. This element has the
+following attributes/sub-elements:
 
   :batches:
     A list of integers separated by spaces indicating at what batches a state
@@ -364,19 +366,52 @@ attributes/sub-elements:
 
     *Default*: None
 
+``<source_point>`` Element
+--------------------------
+
+The ``<source_point>`` element indicates at what batches the source bank 
+should be written. The source bank can be either written out within a state  
+point file or separately in a source point file. This element has the following
+attributes/sub-elements:
+
+  :batches:
+    A list of integers separated by spaces indicating at what batches a state
+    point file should be written. It should be noted that if source_separate
+    tag is not set to "true", this list must be a subset of state point batches.
+
+    *Default*: Last batch only
+
+  :interval:
+    A single integer :math:`n` indicating that a state point should be written
+    every :math:`n` batches. This option can be given in lieu of listing
+    batches explicitly. It should be noted that if source_separate tag is not
+    set to "true", this value should produce a list of batches that is a subset
+    of state point batches.
+
+    *Default*: None
+
   :source_separate:
-    If this element is set to "true", a separate binary source file will be
+    If this element is set to "true", a separate binary source point file will be
     written. Otherwise, the source sites will be written in the state point
     directly.
 
     *Default*: false
 
-  :source_write: If this element is set to "false", source sites are not written
-    to the state point file. This can substantially reduce the size of state
-    points if large numbers of particles per batch are used.
+  :source_write:
+    If this element is set to "false", source sites are not written
+    to the state point or source point file. This can substantially reduce the 
+    size of state points if large numbers of particles per batch are used.
 
     *Default*: true
 
+  :overwrite_latest:
+    If this element is set to "true", a source point file containing
+    the source bank will be written out to a separate file named 
+    ``source.binary`` or ``source.h5`` depending on if HDF5 is enabled. 
+    This file will be overwritten at every single batch so that the latest
+    source bank will be available. It should be noted that a user can set both 
+    this element to "true" and specify batches to write a permanent source bank.
+
 ``<survival_biasing>`` Element
 ------------------------------
 

src/constants.F90

@@ -11,13 +11,14 @@ module constants
   integer, parameter :: VERSION_RELEASE = 3
 
   ! Revision numbers for binary files
-  integer, parameter :: REVISION_STATEPOINT       = 10
+  integer, parameter :: REVISION_STATEPOINT       = 11
   integer, parameter :: REVISION_PARTICLE_RESTART = 1
 
   ! Binary file types
   integer, parameter :: &
        FILETYPE_STATEPOINT       = -1, &
-       FILETYPE_PARTICLE_RESTART = -2
+       FILETYPE_PARTICLE_RESTART = -2, &
+       FILETYPE_SOURCE           = -3
 
   ! ============================================================================
   ! ADJUSTABLE PARAMETERS

src/eigenvalue.F90

@@ -17,7 +17,7 @@ module eigenvalue
   use random_lcg,   only: prn, set_particle_seed, prn_skip
   use search,       only: binary_search
   use source,       only: get_source_particle
-  use state_point,  only: write_state_point
+  use state_point,  only: write_state_point, write_source_point
   use string,       only: to_str
   use tally,        only: synchronize_tallies, setup_active_usertallies, &
                           reset_result
@@ -222,6 +222,12 @@ contains
       call write_state_point()
     end if
 
+    ! Write out source point if it's been specified for this batch
+    if ((sourcepoint_batch % contains(current_batch) .or. source_latest) .and. &
+        source_write) then
+      call write_source_point()
+    end if
+
     if (master .and. current_batch == n_batches) then
       ! Make sure combined estimate of k-effective is calculated at the last
       ! batch in case no state point is written

src/global.F90

@@ -189,6 +189,7 @@ module global
   ! Write source at end of simulation
   logical :: source_separate = .false.
   logical :: source_write = .true.
+  logical :: source_latest = .false.
 
   ! ============================================================================
   ! PARALLEL PROCESSING VARIABLES
@@ -260,6 +261,7 @@ module global
   character(MAX_FILE_LEN) :: path_cross_sections   ! Path to cross_sections.xml
   character(MAX_FILE_LEN) :: path_source = ''      ! Path to binary source
   character(MAX_FILE_LEN) :: path_state_point      ! Path to binary state point
+  character(MAX_FILE_LEN) :: path_source_point     ! Path to binary source point
   character(MAX_FILE_LEN) :: path_particle_restart ! Path to particle restart
   character(MAX_FILE_LEN) :: path_output = ''      ! Path to output directory
 
@@ -363,6 +365,10 @@ module global
   integer :: n_state_points = 0
   type(SetInt) :: statepoint_batch
 
+  ! Information about source points to be written
+  integer :: n_source_points = 0
+  type(SetInt) :: sourcepoint_batch
+
   ! Various output options
   logical :: output_summary = .false.
   logical :: output_xs      = .false.

src/initialize.F90

@@ -363,8 +363,50 @@ contains
           case (FILETYPE_PARTICLE_RESTART)
             path_particle_restart = argv(i)
             particle_restart_run = .true.
+          case default
+            message = "Unrecognized file after restart flag."
+            call fatal_error()
           end select
 
+          ! If its a restart run check for additional source file
+          if (restart_run .and. i + 1 <= argc) then
+
+            ! Increment arg
+            i = i + 1
+
+            ! Check if it has extension we can read
+            if ((ends_with(argv(i), '.binary') .or. &
+                 ends_with(argv(i), '.h5'))) then
+
+              ! Check file type is a source file
+              call sp % file_open(argv(i), 'r', serial = .false.)
+              call sp % read_data(filetype, 'filetype')
+              call sp % file_close()
+              if (filetype /= FILETYPE_SOURCE) then
+                message = "Second file after restart flag must be a source file"
+                call fatal_error()
+              end if
+
+              ! It is a source file
+              path_source_point = argv(i)
+
+            else ! Different option is specified not a source file
+
+              ! Source is in statepoint file
+              path_source_point = path_state_point
+
+              ! Set argument back
+              i = i - 1
+
+            end if
+
+          else ! No command line arg after statepoint
+
+            ! Source is assumed to be in statepoint file
+            path_source_point = path_state_point
+
+          end if
+
         case ('-g', '-geometry-debug', '--geometry-debug')
           check_overlaps = .true.
 

src/input_xml.F90

@@ -625,20 +625,6 @@ contains
         n_state_points = 1
         call statepoint_batch % add(n_batches)
       end if
-
-      ! Check if the user has specified to write binary source file
-      if (check_for_node(node_sp, "source_separate")) then
-        call get_node_value(node_sp, "source_separate", temp_str)
-        call lower_case(temp_str)
-        if (trim(temp_str) == 'true' .or. &
-             trim(temp_str) == '1') source_separate = .true.
-      end if
-      if (check_for_node(node_sp, "source_write")) then
-        call get_node_value(node_sp, "source_write", temp_str)
-        call lower_case(temp_str)
-        if (trim(temp_str) == 'false' .or. &
-             trim(temp_str) == '0') source_write = .false.
-      end if
     else
       ! If no <state_point> tag was present, by default write state point at
       ! last batch only
@@ -646,6 +632,85 @@ contains
       call statepoint_batch % add(n_batches)
     end if
 
+    ! Check if the user has specified to write source points
+    if (check_for_node(doc, "source_point")) then
+
+      ! Get pointer to source_point node
+      call get_node_ptr(doc, "source_point", node_sp)
+
+      ! Determine number of batches at which to store source points
+      if (check_for_node(node_sp, "batches")) then
+        n_source_points = get_arraysize_integer(node_sp, "batches") 
+      else
+        n_source_points = 0
+      end if
+
+      if (n_source_points > 0) then
+        ! User gave specific batches to write source points
+        allocate(temp_int_array(n_source_points))
+        call get_node_array(node_sp, "batches", temp_int_array)
+        do i = 1, n_source_points
+          call sourcepoint_batch % add(temp_int_array(i))
+        end do
+        deallocate(temp_int_array)
+      elseif (check_for_node(node_sp, "interval")) then
+        ! User gave an interval for writing source points
+        call get_node_value(node_sp, "interval", temp_int)
+        n_source_points = n_batches / temp_int
+        do i = 1, n_source_points
+          call sourcepoint_batch % add(temp_int * i)
+        end do
+      else
+        ! If neither were specified, write source points with state points
+        n_source_points = n_state_points
+        do i = 1, n_state_points
+          call sourcepoint_batch % add(statepoint_batch % get_item(i))
+        end do
+      end if
+
+      ! Check if the user has specified to write binary source file
+      if (check_for_node(node_sp, "separate")) then
+        call get_node_value(node_sp, "separate", temp_str)
+        call lower_case(temp_str)
+        if (trim(temp_str) == 'true' .or. &
+             trim(temp_str) == '1') source_separate = .true.
+      end if
+      if (check_for_node(node_sp, "write")) then
+        call get_node_value(node_sp, "write", temp_str)
+        call lower_case(temp_str)
+        if (trim(temp_str) == 'false' .or. &
+             trim(temp_str) == '0') source_write = .false.
+      end if
+      if (check_for_node(node_sp, "overwrite_latest")) then
+        call get_node_value(node_sp, "overwrite_latest", temp_str)
+        call lower_case(temp_str)
+        if (trim(temp_str) == 'true' .or. &
+             trim(temp_str) == '1') source_latest = .true.
+      end if
+    else
+      ! If no <source_point> tag was present, by default we keep source bank in
+      ! statepoint file and write it out at statepoints intervals 
+      source_separate = .false.
+      n_source_points = n_state_points
+      do i = 1, n_state_points
+        call sourcepoint_batch % add(statepoint_batch % get_item(i))
+      end do
+    end if
+
+    ! If source is not seperate and is to be written out in the statepoint file,
+    ! make sure that the sourcepoint batch numbers are contained in the
+    ! statepoint list
+    if (.not. source_separate) then
+      do i = 1, n_source_points
+        if (.not. statepoint_batch % contains(sourcepoint_batch % &
+            get_item(i))) then
+          message = 'Sourcepoint batches are not a subset&
+                    & of statepoint batches.'
+          call fatal_error()
+        end if
+      end do
+    end if
+
     ! Check if the user has specified to not reduce tallies at the end of every
     ! batch
     if (check_for_node(doc, "no_reduce")) then

src/relaxng/settings.rnc

@@ -92,11 +92,22 @@ element settings {
         attribute batches { list { xsd:positiveInteger+ } }) |
       (element interval { xsd:positiveInteger } | 
         attribute interval { xsd:positiveInteger })
-    ) &
-    (element source_separate { xsd:boolean } | 
-      attribute source_separate { xsd:boolean })? &
-    (element source_write { xsd:boolean } | 
-      attribute source_write { xsd:boolean })?
+    )
+  }? &
+
+  element source_point {
+    (
+      (element batches { list { xsd:positiveInteger+ } } |
+        attribute batches { list { xsd:positiveInteger+ } }) |
+      (element interval { xsd:positiveInteger } | 
+        attribute interval { xsd:positiveInteger })
+    )? &
+    (element separate { xsd:boolean } | 
+      attribute separate { xsd:boolean })? &
+    (element write { xsd:boolean } | 
+      attribute write { xsd:boolean })? &
+    (element overwrite_latest { xsd:boolean} |
+      attribute overwrite_latest {xsd:boolean})?
   }? &
 
   element survival_biasing { xsd:boolean }? &

src/state_point.F90

@@ -232,6 +232,13 @@ contains
 
     end if
 
+    ! Indicate where source bank is stored in statepoint
+    if (source_separate) then
+      call sp % write_data(0, "source_present")
+    else
+      call sp % write_data(1, "source_present")
+    end if
+
     ! Check for the no-tally-reduction method
     if (.not. reduce_tallies) then
       ! If using the no-tally-reduction method, we need to collect tally
@@ -280,14 +287,45 @@ contains
 
     end if
 
-    ! Check for eigenvalue calculation
-    if (run_mode == MODE_EIGENVALUE .and. source_write) then
+  end subroutine write_state_point
 
-      ! Check for writing source out separately
+!===============================================================================
+! WRITE_SOURCE_POINT
+!===============================================================================
+
+  subroutine write_source_point()
+
+    type(BinaryOutput) :: sp
+    character(MAX_FILE_LEN) :: filename
+
+    ! Check to write out source for a specified batch
+    if (sourcepoint_batch % contains(current_batch)) then
+
+      ! Create or open up file
       if (source_separate) then
 
-        ! Set filename for source
-        filename = trim(path_output) // 'source.' // &
+        ! Set filename
+        filename = trim(path_output) // 'source.' // trim(to_str(current_batch))
+#ifdef HDF5
+        filename = trim(filename) // '.h5'
+#else
+        filename = trim(filename) // '.binary'
+#endif
+
+        ! Write message for new file creation
+        message = "Creating source file " // trim(filename) // "..."
+        call write_message(1)
+
+        ! Create separate source file
+        call sp % file_create(filename, serial = .false.)
+
+        ! Write file type
+        call sp % write_data(FILETYPE_SOURCE, "filetype")
+
+      else
+
+        ! Set filename for state point
+        filename = trim(path_output) // 'statepoint.' // &
                    trim(to_str(current_batch))
 #ifdef HDF5
         filename = trim(filename) // '.h5'
@@ -295,16 +333,7 @@ contains
         filename = trim(filename) // '.binary'
 #endif
 
-        ! Write message
-        message = "Creating source file " // trim(filename) // "..."
-        call write_message(1)
-
-        ! Create source file 
-        call sp % file_create(filename, serial = .false.)
-
-      else
-
-        ! Reopen state point file in parallel
+        ! Reopen statepoint file in parallel
         call sp % file_open(filename, 'w', serial = .false.)
 
       end if
@@ -317,7 +346,36 @@ contains
 
     end if
 
-  end subroutine write_state_point
+    ! Also check to write source separately in overwritten file
+    if (source_latest) then
+
+      ! Set filename
+      filename = trim(path_output) // 'source'
+#ifdef HDF5
+      filename = trim(filename) // '.h5'
+#else
+      filename = trim(filename) // '.binary'
+#endif
+
+      ! Write message for new file creation
+      message = "Creating source file " // trim(filename) // "..."
+      call write_message(1)
+
+      ! Always create this file because it will be overwritten
+      call sp % file_create(filename, serial = .false.)
+
+      ! Write file type
+      call sp % write_data(FILETYPE_SOURCE, "filetype")
+
+      ! Write out source
+      call sp % write_source_bank()
+
+      ! Close file
+      call sp % file_close()
+
+    end if
+
+  end subroutine write_source_point
 
 !===============================================================================
 ! WRITE_TALLY_RESULTS_NR
@@ -467,6 +525,7 @@ contains
     integer                 :: length(4)
     integer                 :: int_array(3)
     integer, allocatable    :: temp_array(:)
+    logical                 :: source_present
     real(8)                 :: real_array(3) 
     type(TallyObject), pointer :: t => null()
 
@@ -671,6 +730,20 @@ contains
 
     end do TALLY_METADATA
 
+    ! Check for source in statepoint if needed
+    call sp % read_data(int_array(1), "source_present")
+    if (int_array(1) == 1) then
+      source_present = .true.
+    else
+      source_present = .false.
+    end if
+
+    ! Check to make sure source bank is present
+    if (path_source_point == path_state_point .and. .not. source_present) then
+      message = "Source bank must be contained in statepoint restart file"
+      call fatal_error()
+    end if 
+
     ! Read tallies to master
     if (master) then
 
@@ -711,26 +784,20 @@ contains
     if (run_mode == MODE_EIGENVALUE) then
 
       ! Check if source was written out separately
-      if (source_separate) then
+      if (.not. source_present) then
 
         ! Close statepoint file 
         call sp % file_close()
 
-        ! Set filename for source
-        filename = trim(path_output) // 'source.' // &
-                   trim(to_str(restart_batch))
-#ifdef HDF5
-        filename = trim(filename) // '.h5'
-#else
-        filename = trim(filename) // '.binary'
-#endif
-
         ! Write message
         message = "Loading source file " // trim(filename) // "..."
         call write_message(1)
 
         ! Open source file 
-        call sp % file_open(filename, 'r', serial = .false.)
+        call sp % file_open(path_source_point, 'r', serial = .false.)
+
+        ! Read file type
+        call sp % read_data(int_array(1), "filetype")
 
       end if
 

src/utils/statepoint.py

@@ -151,7 +151,7 @@ class StatePoint(object):
 
         # Read statepoint revision
         self.revision = self._get_int(path='revision')[0]
-        if self.revision != 10:
+        if self.revision != 11:
           raise Exception('Statepoint Revision is not consistent.')
 
         # Read OpenMC version
@@ -298,6 +298,13 @@ class StatePoint(object):
                 f.stride = stride
                 stride *= f.length
 
+        # Source bank present
+        source_present = self._get_int(path='source_present')[0]
+        if source_present == 1:
+            self.source_present = True       
+        else:
+            self.source_present = False
+
         # Set flag indicating metadata has already been read
         self._metadata = True
 
@@ -342,6 +349,11 @@ class StatePoint(object):
         if not self._results:
             self.read_results()
 
+        # Check if source bank is in statepoint
+        if not self.source_present:
+            print('Source not in statepoint file.')
+            return
+
         # For HDF5 state points, copy entire bank
         if self._hdf5:
             source_sites = self._f['source_bank'].value

tests/cleanup

@@ -0,0 +1,10 @@
+#!/bin/bash
+
+# This simple script ensures that all binary
+# output files have been deleted in all the
+# folders. This can occur if a previous error
+# occurred and the test suite was rerun without
+# deleting left over binary files. This will
+# cause an assertion error in some of the 
+# tests.
+find . \( -name "*.binary" -o -name "*.h5" \) -exec rm -f {} \;

tests/run_tests.py

@@ -117,6 +117,8 @@ if len(sys.argv) > 1:
                     if j.rfind(suffix) != -1:
                         if suffix == 'omp' and j == 'compile':
                             continue
+                        if j == 'compile':
+                            continue 
                         tests__.append(j)
         else:
             tests__.append(i)  # append specific test (e.g., mpi-debug)

tests/test_sourcepoint_batch/geometry.xml

@@ -0,0 +1,8 @@
+<?xml version="1.0"?>
+<geometry>
+
+  <!-- Sphere with radius 10 -->
+  <surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
+  <cell id="1" material="1" surfaces="-1" />
+    
+</geometry>

tests/test_sourcepoint_batch/materials.xml

@@ -0,0 +1,9 @@
+<?xml version="1.0"?>
+<materials>
+
+  <material id="1">
+    <density value="4.5" units="g/cc" />
+    <nuclide name="U-235" xs="70c" ao="1.0" />
+  </material>
+
+</materials>

tests/test_sourcepoint_batch/results.py

@@ -0,0 +1,32 @@
+#!/usr/bin/env python
+
+import sys
+
+# import statepoint
+sys.path.append('../../src/utils')
+import statepoint
+
+# read in statepoint file
+if len(sys.argv) > 1:
+    sp = statepoint.StatePoint(sys.argv[1])
+else:
+    sp = statepoint.StatePoint('statepoint.8.binary')
+sp.read_results()
+sp.read_source()
+
+# set up output string
+outstr = ''
+ 
+# write out k-combined
+outstr += 'k-combined:\n'
+outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
+
+# write out xyz
+xyz = sp.source[0].xyz
+for i in xyz:
+    outstr += "{0:12.6E} ".format(i)
+outstr += "\n"
+
+# write results to file
+with open('results_test.dat','w') as fh:
+    fh.write(outstr)

tests/test_sourcepoint_batch/results_true.dat

@@ -0,0 +1,3 @@
+k-combined:
+0.000000E+00 0.000000E+00
+-9.438655E-02 -4.436810E+00 -2.416825E+00 

tests/test_sourcepoint_batch/settings.xml

@@ -0,0 +1,19 @@
+<?xml version="1.0"?>
+<settings>
+
+  <state_point batches="2 3 4 5 8"/>
+  <source_point batches="2 5 8"/>
+
+  <eigenvalue>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <particles>1000</particles>
+  </eigenvalue>
+
+  <source>
+    <space type="box">
+      <parameters>-4 -4 -4  4  4  4</parameters>
+    </space>
+  </source>
+
+</settings>

tests/test_sourcepoint_batch/test_sourcepoint_batch.py

@@ -0,0 +1,44 @@
+#!/usr/bin/env python
+
+import os
+from subprocess import Popen, STDOUT, PIPE, call
+import filecmp
+from nose_mpi import NoseMPI
+import glob
+
+pwd = os.path.dirname(__file__)
+
+def setup(): 
+    os.putenv('PWD', pwd)
+    os.chdir(pwd)
+
+def test_run():
+    openmc_path = pwd + '/../../src/openmc'
+    if int(NoseMPI.mpi_np) > 0:
+        proc = Popen([NoseMPI.mpi_exec, '-np', NoseMPI.mpi_np, openmc_path],
+               stderr=STDOUT, stdout=PIPE)
+    else:
+        proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
+    print(proc.communicate()[0])
+    returncode = proc.returncode
+    assert returncode == 0
+
+def test_statepoint_exists():
+    statepoint = glob.glob(pwd + '/statepoint.*')
+    assert len(statepoint) == 5
+    assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5')
+
+def test_results():
+    statepoint = glob.glob(pwd + '/statepoint.8.*')
+    call(['python', 'results.py', statepoint[0]])
+    compare = filecmp.cmp('results_test.dat', 'results_true.dat')
+    if not compare:
+      os.rename('results_test.dat', 'results_error.dat')
+    assert compare
+
+def teardown():
+    output = glob.glob(pwd + '/statepoint.*')
+    output.append(pwd + '/results_test.dat')
+    for f in output:
+        if os.path.exists(f):
+            os.remove(f)

tests/test_sourcepoint_interval/geometry.xml

@@ -0,0 +1,8 @@
+<?xml version="1.0"?>
+<geometry>
+
+  <!-- Sphere with radius 10 -->
+  <surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
+  <cell id="1" material="1" surfaces="-1" />
+    
+</geometry>

tests/test_sourcepoint_interval/materials.xml

@@ -0,0 +1,9 @@
+<?xml version="1.0"?>
+<materials>
+
+  <material id="1">
+    <density value="4.5" units="g/cc" />
+    <nuclide name="U-235" xs="70c" ao="1.0" />
+  </material>
+
+</materials>

tests/test_sourcepoint_interval/results.py

@@ -0,0 +1,32 @@
+#!/usr/bin/env python
+
+import sys
+
+# import statepoint
+sys.path.append('../../src/utils')
+import statepoint
+
+# read in statepoint file
+if len(sys.argv) > 1:
+    sp = statepoint.StatePoint(sys.argv[1])
+else:
+    sp = statepoint.StatePoint('statepoint.8.binary')
+sp.read_results()
+sp.read_source()
+
+# set up output string
+outstr = ''
+ 
+# write out k-combined
+outstr += 'k-combined:\n'
+outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
+
+# write out xyz
+xyz = sp.source[0].xyz
+for i in xyz:
+    outstr += "{0:12.6E} ".format(i)
+outstr += "\n"
+
+# write results to file
+with open('results_test.dat','w') as fh:
+    fh.write(outstr)

tests/test_sourcepoint_interval/results_true.dat

@@ -0,0 +1,3 @@
+k-combined:
+0.000000E+00 0.000000E+00
+-9.438655E-02 -4.436810E+00 -2.416825E+00 

tests/test_sourcepoint_interval/settings.xml

@@ -0,0 +1,19 @@
+<?xml version="1.0"?>
+<settings>
+
+  <state_point interval="2"/>
+  <source_point interval="4"/>
+
+  <eigenvalue>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <particles>1000</particles>
+  </eigenvalue>
+
+  <source>
+    <space type="box">
+      <parameters>-4 -4 -4  4  4  4</parameters>
+    </space>
+  </source>
+
+</settings>

tests/test_sourcepoint_interval/test_sourcepoint_interval.py

@@ -0,0 +1,44 @@
+#!/usr/bin/env python
+
+import os
+from subprocess import Popen, STDOUT, PIPE, call
+import filecmp
+from nose_mpi import NoseMPI
+import glob
+
+pwd = os.path.dirname(__file__)
+
+def setup(): 
+    os.putenv('PWD', pwd)
+    os.chdir(pwd)
+
+def test_run():
+    openmc_path = pwd + '/../../src/openmc'
+    if int(NoseMPI.mpi_np) > 0:
+        proc = Popen([NoseMPI.mpi_exec, '-np', NoseMPI.mpi_np, openmc_path],
+               stderr=STDOUT, stdout=PIPE)
+    else:
+        proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
+    print(proc.communicate()[0])
+    returncode = proc.returncode
+    assert returncode == 0
+
+def test_statepoint_exists():
+    statepoint = glob.glob(pwd + '/statepoint.*')
+    assert len(statepoint) == 5
+    assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5')
+
+def test_results():
+    statepoint = glob.glob(pwd + '/statepoint.8.*')
+    call(['python', 'results.py', statepoint[0]])
+    compare = filecmp.cmp('results_test.dat', 'results_true.dat')
+    if not compare:
+      os.rename('results_test.dat', 'results_error.dat')
+    assert compare
+
+def teardown():
+    output = glob.glob(pwd + '/statepoint.*')
+    output.append(pwd + '/results_test.dat')
+    for f in output:
+        if os.path.exists(f):
+            os.remove(f)

tests/test_sourcepoint_latest/geometry.xml

@@ -0,0 +1,8 @@
+<?xml version="1.0"?>
+<geometry>
+
+  <!-- Sphere with radius 10 -->
+  <surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
+  <cell id="1" material="1" surfaces="-1" />
+    
+</geometry>

tests/test_sourcepoint_latest/materials.xml

@@ -0,0 +1,9 @@
+<?xml version="1.0"?>
+<materials>
+
+  <material id="1">
+    <density value="4.5" units="g/cc" />
+    <nuclide name="U-235" xs="70c" ao="1.0" />
+  </material>
+
+</materials>

tests/test_sourcepoint_latest/results.py

@@ -0,0 +1,25 @@
+#!/usr/bin/env python
+
+import sys
+
+# import statepoint
+sys.path.append('../../src/utils')
+import statepoint
+
+# read in statepoint file
+if len(sys.argv) > 1:
+    sp = statepoint.StatePoint(sys.argv[1])
+else:
+    sp = statepoint.StatePoint('statepoint.10.binary')
+sp.read_results()
+
+# set up output string
+outstr = ''
+ 
+# write out k-combined
+outstr += 'k-combined:\n'
+outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
+
+# write results to file
+with open('results_test.dat','w') as fh:
+    fh.write(outstr)

tests/test_sourcepoint_latest/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+3.011353E-01 2.854556E-03

tests/test_sourcepoint_latest/settings.xml

@@ -0,0 +1,18 @@
+<?xml version="1.0"?>
+<settings>
+
+  <source_point batches="0" separate="true" overwrite_latest="true"/>
+
+  <eigenvalue>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <particles>1000</particles>
+  </eigenvalue>
+
+  <source>
+    <space type="box">
+      <parameters>-4 -4 -4  4  4  4</parameters>
+    </space>
+  </source>
+
+</settings>

tests/test_sourcepoint_latest/test_sourcepoint_latest.py

@@ -0,0 +1,48 @@
+#!/usr/bin/env python
+
+import os
+from subprocess import Popen, STDOUT, PIPE, call
+import filecmp
+from nose_mpi import NoseMPI
+import glob
+
+pwd = os.path.dirname(__file__)
+
+def setup(): 
+    os.putenv('PWD', pwd)
+    os.chdir(pwd)
+
+def test_run():
+    openmc_path = pwd + '/../../src/openmc'
+    if int(NoseMPI.mpi_np) > 0:
+        proc = Popen([NoseMPI.mpi_exec, '-np', NoseMPI.mpi_np, openmc_path],
+               stderr=STDOUT, stdout=PIPE)
+    else:
+        proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
+    print(proc.communicate()[0])
+    returncode = proc.returncode
+    assert returncode == 0
+
+def test_statepoint_exists():
+    statepoint = glob.glob(pwd + '/statepoint.10.*')
+    assert len(statepoint) == 1
+    assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5')
+    source = glob.glob(pwd + '/source.*')
+    assert len(source) == 1
+    assert source[0].endswith('binary') or source[0].endswith('h5')
+
+def test_results():
+    statepoint = glob.glob(pwd + '/statepoint.10.*')
+    call(['python', 'results.py', statepoint[0]])
+    compare = filecmp.cmp('results_test.dat', 'results_true.dat')
+    if not compare:
+      os.rename('results_test.dat', 'results_error.dat')
+    assert compare
+
+def teardown():
+    output = glob.glob(pwd + '/statepoint.10.*')
+    output += glob.glob(pwd + '/source.*')
+    output.append(pwd + '/results_test.dat')
+    for f in output:
+        if os.path.exists(f):
+            os.remove(f)

tests/test_sourcepoint_restart/geometry.xml

@@ -0,0 +1,8 @@
+<?xml version="1.0"?>
+<geometry>
+
+  <!-- Sphere with radius 10 -->
+  <surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
+  <cell id="1" material="1" surfaces="-1" />
+    
+</geometry>

tests/test_sourcepoint_restart/materials.xml

@@ -0,0 +1,9 @@
+<?xml version="1.0"?>
+<materials>
+
+  <material id="1">
+    <density value="4.5" units="g/cc" />
+    <nuclide name="U-235" xs="70c" ao="1.0" />
+  </material>
+
+</materials>

tests/test_sourcepoint_restart/results.py

@@ -0,0 +1,44 @@
+#!/usr/bin/env python
+
+import sys
+import numpy as np
+
+# import statepoint
+sys.path.append('../../src/utils')
+import statepoint
+
+# read in statepoint file
+if len(sys.argv) > 1:
+    sp = statepoint.StatePoint(sys.argv[1])
+else:
+    sp = statepoint.StatePoint('statepoint.10.binary')
+sp.read_results()
+
+# extract tally results and convert to vector
+results1 = sp.tallies[0].results
+shape1 = results1.shape
+size1 = (np.product(shape1))
+results1 = np.reshape(results1, size1)
+results2 = sp.tallies[1].results
+shape2 = results2.shape
+size2 = (np.product(shape2))
+results2 = np.reshape(results2, size2)
+
+# set up output string
+outstr = ''
+ 
+# write out k-combined
+outstr += 'k-combined:\n'
+outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
+
+# write out tally results
+outstr += 'tally 1:\n'
+for item in results1:
+  outstr += "{0:12.6E}\n".format(item)
+outstr += 'tally 2:\n'
+for item in results2:
+  outstr += "{0:12.6E}\n".format(item)
+
+# write results to file
+with open('results_test.dat','w') as fh:
+    fh.write(outstr)

tests/test_sourcepoint_restart/settings.xml

@@ -0,0 +1,19 @@
+<?xml version="1.0"?>
+<settings>
+
+  <state_point batches="7 10" />
+  <source_point batches="7" separate="true" />
+
+  <eigenvalue>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <particles>1000</particles>
+  </eigenvalue>
+
+  <source>
+    <space type="box">
+      <parameters>-4 -4 -4  4  4  4</parameters>
+    </space>
+  </source>
+
+</settings>

tests/test_sourcepoint_restart/tallies.xml

@@ -0,0 +1,23 @@
+<?xml version="1.0"?>
+<tallies>
+
+  <mesh id="1">
+    <type>rectangular</type>
+    <dimension>5 3 4</dimension>
+    <lower_left>-10. -5. 0.</lower_left>
+    <upper_right>10. 4. 9.</upper_right>
+  </mesh>
+
+  <tally id="10">
+    <filter type="mesh" bins="1" />
+    <filter type="energy" bins="0.0 5.0 10.0" />
+    <filter type="energyout" bins="0.0 5.0 10.0" />
+    <scores>scatter-P3 nu-fission</scores>
+  </tally>
+
+  <tally id="5">
+    <filter type="cell" bins="1" />
+    <scores>fission absorption total flux</scores>
+  </tally>
+
+</tallies>

tests/test_sourcepoint_restart/test_sourcepoint_restart.py

@@ -0,0 +1,129 @@
+#!/usr/bin/env python
+
+import os
+from subprocess import Popen, STDOUT, PIPE, call
+import filecmp
+from nose_mpi import NoseMPI
+import glob
+
+pwd = os.path.dirname(__file__)
+
+def setup(): 
+    os.putenv('PWD', pwd)
+    os.chdir(pwd)
+
+def test_run():
+    openmc_path = pwd + '/../../src/openmc'
+    if int(NoseMPI.mpi_np) > 0:
+        proc = Popen([NoseMPI.mpi_exec, '-np', NoseMPI.mpi_np, openmc_path],
+               stderr=STDOUT, stdout=PIPE)
+    else:
+        proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
+    print(proc.communicate()[0])
+    returncode = proc.returncode
+    assert returncode == 0
+
+def test_created_statepoint():
+    statepoint = glob.glob(pwd + '/statepoint.*')
+    assert len(statepoint) == 2
+    assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5')
+    sourcepoint = glob.glob(pwd + '/source.7.*')
+    assert len(sourcepoint) == 1
+    assert sourcepoint[0].endswith('binary') or sourcepoint[0].endswith('h5')
+
+def test_results():
+    statepoint = glob.glob(pwd + '/statepoint.10.*')
+    call(['python', 'results.py', statepoint[0]])
+    compare = filecmp.cmp('results_test.dat', 'results_true.dat')
+    if not compare:
+      os.rename('results_test.dat', 'results_error.dat')
+    assert compare
+    os.remove(statepoint[0])
+
+def test_restart_form1():
+    statepoint = glob.glob(pwd + '/statepoint.7.*')
+    sourcepoint = glob.glob(pwd + '/source.7.*')
+    openmc_path = pwd + '/../../src/openmc'
+    if int(NoseMPI.mpi_np) > 0:
+        proc = Popen([NoseMPI.mpi_exec, '-np', NoseMPI.mpi_np, openmc_path, 
+                      '-r', statepoint[0], sourcepoint[0]],
+               stderr=STDOUT, stdout=PIPE)
+    else:
+        proc = Popen([openmc_path, '-r', statepoint[0], sourcepoint[0]], stderr=STDOUT, stdout=PIPE)
+    print(proc.communicate()[0])
+    returncode = proc.returncode
+    assert returncode == 0
+
+def test_created_statepoint_form1():
+    statepoint = glob.glob(pwd + '/statepoint.10.*')
+    assert len(statepoint) == 1
+    assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5')
+
+def test_results_form1():
+    statepoint = glob.glob(pwd + '/statepoint.10.*')
+    call(['python', 'results.py', statepoint[0]])
+    compare = filecmp.cmp('results_test.dat', 'results_true.dat')
+    if not compare:
+      os.rename('results_test.dat', 'results_error.dat')
+    assert compare
+    os.remove(statepoint[0])
+
+def test_restart_form2():
+    statepoint = glob.glob(pwd + '/statepoint.7.*')
+    sourcepoint = glob.glob(pwd + '/source.7.*')
+    openmc_path = pwd + '/../../src/openmc'
+    if int(NoseMPI.mpi_np) > 0:
+        proc = Popen([NoseMPI.mpi_exec, '-np', NoseMPI.mpi_np, openmc_path, 
+                      '--restart', statepoint[0], sourcepoint[0]],
+               stderr=STDOUT, stdout=PIPE)
+    else:
+        proc = Popen([openmc_path, '--restart', statepoint[0], sourcepoint[0]],
+                     stderr=STDOUT, stdout=PIPE)
+    print(proc.communicate()[0])
+    returncode = proc.returncode
+    assert returncode == 0
+
+def test_created_statepoint_form2():
+    statepoint = glob.glob(pwd + '/statepoint.10.*')
+    assert len(statepoint) == 1
+    assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5')
+
+def test_results_form2():
+    statepoint = glob.glob(pwd + '/statepoint.10.*')
+    call(['python', 'results.py', statepoint[0]])
+    compare = filecmp.cmp('results_test.dat', 'results_true.dat')
+    if not compare:
+      os.rename('results_test.dat', 'results_error.dat')
+    assert compare
+    os.remove(statepoint[0])
+
+def test_restart_serial():
+    statepoint = glob.glob(pwd + '/statepoint.7.*')
+    sourcepoint = glob.glob(pwd + '/source.7.*')
+    openmc_path = pwd + '/../../src/openmc'
+    proc = Popen([openmc_path, '--restart', statepoint[0], sourcepoint[0]], 
+                 stderr=STDOUT, stdout=PIPE)
+    print(proc.communicate()[0])
+    returncode = proc.returncode
+    assert returncode == 0
+
+def test_created_statepoint_serial():
+    statepoint = glob.glob(pwd + '/statepoint.10.*')
+    assert len(statepoint) == 1
+    assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5')
+
+def test_results_serial():
+    statepoint = glob.glob(pwd + '/statepoint.10.*')
+    call(['python', 'results.py', statepoint[0]])
+    compare = filecmp.cmp('results_test.dat', 'results_true.dat')
+    if not compare:
+      os.rename('results_test.dat', 'results_error.dat')
+    assert compare
+
+def teardown():
+    output = glob.glob(pwd + '/statepoint.*')
+    output += glob.glob(pwd + '/source.*')
+    output.append(pwd + '/results_test.dat')
+    for f in output:
+        if os.path.exists(f):
+            os.remove(f)

tests/test_statepoint_sourcesep/settings.xml

@@ -1,7 +1,8 @@
 <?xml version="1.0"?>
 <settings>
 
-  <state_point batches="10" source_separate="true" />
+  <state_point batches="10" />
+  <source_point separate="true" />
 
   <eigenvalue>
     <batches>10</batches>