Changed state point routine names and fixed Python scripts.

src/criticality.F90

@@ -7,7 +7,7 @@ module criticality
   use output,      only: write_message, header, print_columns
   use physics,     only: transport
   use source,      only: get_source_particle, write_source_binary
-  use state_point, only: create_state_point, replay_batch_history
+  use state_point, only: write_state_point, replay_batch_history
   use string,      only: to_str
   use tally,       only: synchronize_tallies
   use timing,      only: timer_start, timer_stop
@@ -152,7 +152,7 @@ contains
     do i = 1, n_state_points
        if (current_batch == statepoint_batch(i)) then
           ! Create state point file
-          call create_state_point()
+          call write_state_point()
           exit
        end if
     end do

src/fixed_source.F90

@@ -7,7 +7,7 @@ module fixed_source
   use random_lcg,  only: set_particle_seed
   use source,      only: initialize_particle, sample_external_source, &
                          copy_source_attributes
-  use state_point, only: create_state_point
+  use state_point, only: write_state_point
   use string,      only: to_str
   use tally,       only: synchronize_tallies
   use timing,      only: timer_start, timer_stop
@@ -115,7 +115,7 @@ contains
     if (master) then
        do i = 1, n_state_points
           if (current_batch == statepoint_batch(i)) then
-             call create_state_point()
+             call write_state_point()
              exit
           end if
        end do

src/state_point.F90

@@ -16,13 +16,13 @@ module state_point
 contains
 
 !===============================================================================
-! CREATE_STATE_POINT creates a state point binary file that can be used for
+! WRITE_STATE_POINT creates a state point binary file that can be used for
 ! restarting a run or for getting intermediate tally results
 !===============================================================================
 
-  subroutine create_state_point()
+  subroutine write_state_point()
 
-    integer :: i ! loo pindex
+    integer :: i ! loop index
     integer :: n ! temporary array length
     type(TallyObject), pointer :: t => null()
 
@@ -41,7 +41,7 @@ contains
 
     ! Write message
     message = "Creating state point " // trim(path_state_point) // "..."
-    call write_message()
+    call write_message(1)
 
 #ifdef MPI
     ! ==========================================================================
@@ -55,7 +55,7 @@ contains
        ! =======================================================================
        ! RUN INFORMATION AND TALLY METADATA
 
-       call state_point_header(fh)
+       call write_state_point_header(fh)
 
        ! =======================================================================
        ! TALLY RESULTS
@@ -243,14 +243,14 @@ contains
     close(UNIT_STATE)
 #endif
 
-  end subroutine create_state_point
+  end subroutine write_state_point
 
 #ifdef MPI
 !===============================================================================
-! STATE_POINT_HEADER
+! WRITE_STATE_POINT_HEADER
 !===============================================================================
 
-  subroutine state_point_header(fh)
+  subroutine write_state_point_header(fh)
 
     integer, intent(inout) :: fh ! file handle
 
@@ -405,7 +405,7 @@ contains
             MPI_INTEGER, MPI_STATUS_IGNORE, mpi_err)
     end do TALLY_METADATA
 
-  end subroutine state_point_header
+  end subroutine write_state_point_header
 #endif
 
 !===============================================================================

src/utils/statepoint.py

@@ -181,14 +181,19 @@ class StatePoint(BinaryFile):
         if not self._metadata:
             self._read_metadata()
 
-        read_tallies = self._get_int()[0]
+        # Flag indicating if tallies are present
+        tallies_present = self._get_int()[0]
 
-        if read_tallies:
+        # Read tally results
+        if tallies_present:
             for t in self.tallies:
                 n = t.n_score_bins * t.n_filter_bins
                 t.values = np.array(self._get_double(2*n))
                 t.values.shape = (t.n_filter_bins, t.n_score_bins, 2)
 
+        # Indicate that tally values have been read
+        self._values = True
+
     def read_source(self):
         # Check whether tally values have been read
         if not self._values:

src/utils/statepoint_cmp.py

@@ -27,15 +27,18 @@ assert sp1.seed == sp2.seed
 assert sp1.run_mode == sp2.run_mode
 assert sp1.n_particles == sp2.n_particles
 assert sp1.n_batches == sp2.n_batches
-assert sp1.n_inactive == sp2.n_inactive
-assert sp1.gen_per_batch == sp2.gen_per_batch
-assert sp1.current_batch == sp2.current_batch
 
-# Compare keff results
-assert_allclose(sp1.k_batch, sp2.k_batch)
+if sp1.run_mode == 2:
+    # Compare criticality information
+    assert sp1.n_inactive == sp2.n_inactive
+    assert sp1.gen_per_batch == sp2.gen_per_batch
+    assert sp1.current_batch == sp2.current_batch
 
-# Compare entropy results
-assert_allclose(sp1.entropy, sp2.entropy)
+    # Compare keff results
+    assert_allclose(sp1.k_batch, sp2.k_batch)
+
+    # Compare entropy results
+    assert_allclose(sp1.entropy, sp2.entropy)
 
 # Compare global tallies
 assert_allclose(sp1.global_tallies, sp2.global_tallies)
@@ -69,3 +72,15 @@ for t1, t2 in zip(sp1.tallies, sp2.tallies):
 
     # Compare tally results
     assert_allclose(t1.values, t2.values)
+
+# If criticality, compare source sites
+if sp1.run_mode == 2:
+    sp1.read_source()
+    sp2.read_source()
+
+    assert len(sp1.source) == len(sp2.source)
+    for s1, s2 in zip(sp1.source, sp2.source):
+        assert s1.weight == s2.weight
+        assert s1.xyz == s2.xyz
+        assert s1.uvw == s2.uvw
+        assert s1.E == s2.E

src/utils/statepoint_histogram.py

@@ -1,7 +1,6 @@
 #!/usr/bin/env python
 
 from sys import argv
-from struct import unpack
 from math import sqrt
 
 import scipy.stats
@@ -19,6 +18,10 @@ score = int(argv[3]) if len(argv) > 3 else 1
 # Create StatePoint object
 sp = StatePoint(filename)
 
+# Check if tallies are present
+if not sp._get_int()[0]:
+    print("No tally data in state point!")
+
 # Calculate t-value for 95% two-sided CI
 n = sp.current_batch - sp.n_inactive
 t_value = scipy.stats.t.ppf(0.975, n - 1)
@@ -35,7 +38,7 @@ for i, t in enumerate(sp.tallies):
     # Loop over filter/score bins
     for j in range(n_bins):
         # Read sum and sum-squared
-        s, s2 = unpack('=2d', sp._f.read(16))
+        s, s2 = sp._get_double(2)
         s /= n
         if s != 0.0:
             relative_error = t_value*sqrt((s2/n - s*s)/(n-1))/s

src/utils/statepoint_meshplot.py

@@ -1,7 +1,6 @@
 #!/usr/bin/env python
 
 from sys import argv
-from struct import unpack
 from math import sqrt
 
 import numpy as np
@@ -20,6 +19,10 @@ score = int(argv[3]) if len(argv) > 3 else 1
 # Create StatePoint object
 sp = StatePoint(filename)
 
+# Check if tallies are present
+if not sp._get_int()[0]:
+    print("No tally data in state point!")
+
 # Calculate t-value for 95% two-sided CI
 n = sp.current_batch - sp.n_inactive
 t_value = scipy.stats.t.ppf(0.975, n - 1)
@@ -51,7 +54,7 @@ for t in sp.tallies:
             sp._f.seek(start + x*ny*nz*ns*16 + y*nz*ns*16)
 
             # Read sum and sum-squared
-            s, s2 = unpack('=2d', sp._f.read(16))
+            s, s2 = sp._get_double(2)
             s /= n
             mean[x,y] = s
             if s != 0.0: