From 02890579e6760860c73d965cb0d8386f2c186b33 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 23 Aug 2012 14:50:34 -0400 Subject: [PATCH] Changed state point routine names and fixed Python scripts. --- src/criticality.F90 | 4 ++-- src/fixed_source.F90 | 4 ++-- src/state_point.F90 | 18 +++++++++--------- src/utils/statepoint.py | 9 +++++++-- src/utils/statepoint_cmp.py | 29 ++++++++++++++++++++++------- src/utils/statepoint_histogram.py | 7 +++++-- src/utils/statepoint_meshplot.py | 7 +++++-- 7 files changed, 52 insertions(+), 26 deletions(-) diff --git a/src/criticality.F90 b/src/criticality.F90 index 5c833b72f4..229ba48659 100644 --- a/src/criticality.F90 +++ b/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 diff --git a/src/fixed_source.F90 b/src/fixed_source.F90 index ead8676e64..dd2e60601c 100644 --- a/src/fixed_source.F90 +++ b/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 diff --git a/src/state_point.F90 b/src/state_point.F90 index 9c6350feec..2b29247f82 100644 --- a/src/state_point.F90 +++ b/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 !=============================================================================== diff --git a/src/utils/statepoint.py b/src/utils/statepoint.py index 45db158b8a..bd549771a1 100644 --- a/src/utils/statepoint.py +++ b/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: diff --git a/src/utils/statepoint_cmp.py b/src/utils/statepoint_cmp.py index cd55cbfd22..42ac1e7194 100755 --- a/src/utils/statepoint_cmp.py +++ b/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 diff --git a/src/utils/statepoint_histogram.py b/src/utils/statepoint_histogram.py index c1144d7251..dae1e81dee 100755 --- a/src/utils/statepoint_histogram.py +++ b/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 diff --git a/src/utils/statepoint_meshplot.py b/src/utils/statepoint_meshplot.py index 0b2278c96c..d68b3b6cfa 100755 --- a/src/utils/statepoint_meshplot.py +++ b/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: