Added no_reduce option and introduced n_realizations and total_weight global variables. No_reduce does not yet work properly.

src/global.F90

@@ -117,6 +117,10 @@ module global
   integer :: n_tracklength_tallies = 0 ! # of track-length tallies
   integer :: n_current_tallies     = 0 ! # of surface current tallies
 
+  ! Normalization for statistics
+  integer :: n_realizations ! # of independent realizations
+  real(8) :: total_weight   ! total starting particle weight in realization
+
   ! Flag for turning tallies on
   logical :: tallies_on
 
@@ -174,6 +178,9 @@ module global
   integer :: mpi_err     ! MPI error code
   integer :: MPI_BANK    ! MPI datatype for fission bank
 
+  ! No reduction at end of batch
+  logical :: no_reduce = .false.
+
   ! ============================================================================
   ! TIMING VARIABLES
 

src/input_xml.F90

@@ -68,6 +68,7 @@ contains
     energy_grid_ = "union"
     seed_ = 0_8
     write_source_ = ""
+    no_reduce_ = ""
     source_ % type = ""
 
     ! Parse settings.xml file
@@ -293,6 +294,18 @@ contains
     ! Check if the user has specified to write binary source file
     if (trim(write_source_) == 'on') write_source = .true.
 
+    ! Check if the user has specified to not reduce tallies at the end of every
+    ! batch
+    if (trim(no_reduce_) == 'on') no_reduce = .true.
+
+    ! Determine number of realizations
+    if (no_reduce) then
+       n_realizations = n_active * n_procs
+    else
+       n_realizations = n_active
+    end if
+       
+
   end subroutine read_settings_xml
 
 !===============================================================================

src/main.F90

@@ -72,6 +72,9 @@ contains
        message = "Simulating batch " // trim(to_str(current_batch)) // "..."
        call write_message(8)
 
+       ! Reset total starting weight
+       total_weight = ZERO
+
        ! =======================================================================
        ! LOOP OVER GENERATIONS
        GENERATION_LOOP: do current_gen = 1, gen_per_batch

src/source.F90

@@ -172,6 +172,9 @@ contains
     if (current_batch == trace_batch .and. current_gen == trace_gen .and. &
          p % id == trace_particle) trace = .true.
 
+    ! Add paricle's starting weight to count for normalizing tallies later
+    total_weight = total_weight + src % wgt
+
   end subroutine get_source_particle
 
 !===============================================================================

src/tally.F90

@@ -1310,8 +1310,12 @@ contains
     type(TallyObject), pointer :: t => null()
 
 #ifdef MPI
-    call reduce_tallies()
-    if (.not. master) return
+    if (no_reduce) then
+       if (current_batch == n_batches) call reduce_tallies()
+    else
+       call reduce_tallies()
+       if (.not. master) return
+    end if
 #endif
 
     do i = 1, n_tallies
@@ -1369,6 +1373,19 @@ contains
 
     end do
 
+    ! We also need to determine the total starting weight of particles from the
+    ! last realization
+    if (.not. no_reduce) then
+       if (master) then
+          call MPI_REDUCE(MPI_IN_PLACE, total_weight, 1, MPI_REAL8, MPI_SUM, &
+               0, MPI_COMM_WORLD, mpi_err)
+       else
+          call MPI_REDUCE(total_weight, total_weight, 1, MPI_REAL8, MPI_SUM, &
+               0, MPI_COMM_WORLD, mpi_err)
+       end if
+    end if
+
+
   end subroutine reduce_tallies
 #endif
 
@@ -1808,7 +1825,7 @@ contains
     ! within a cycle to the variables sum and sum_sq. This will later allow us
     ! to calculate a variance on the tallies
 
-    val = score % value/(n_particles*gen_per_batch)
+    val = score % value/total_weight
     score % sum    = score % sum    + val
     score % sum_sq = score % sum_sq + val*val
 
@@ -1830,9 +1847,9 @@ contains
     ! have used Bessel's correction so that the estimator of the variance of the
     ! sample mean is unbiased.
 
-    score % sum    = score % sum/n_active
-    score % sum_sq = sqrt((score % sum_sq/n_active - score % sum**2) / &
-         (n_active - 1))
+    score % sum    = score % sum/n_realizations
+    score % sum_sq = sqrt((score % sum_sq/n_realizations - score % sum * &
+         score % sum) / (n_realizations - 1))
 
   end subroutine statistics_score
 

src/xml-fortran/templates/settings_t.xml

@@ -31,6 +31,7 @@
   <variable name="cutoff_" tag="cutoff" type="cutoff_xml" dimension="1" />
   <variable name="energy_grid_" tag="energy_grid" type="word" length="7" />
   <variable name="entropy_" tag="entropy" type="mesh_xml" dimension="1" />
+  <variable name="no_reduce_" tag="no_reduce" type="word" length="3" />
   <variable name="ptables_" tag="ptables" type="word" length="3" />
   <variable name="seed_" tag="seed" type="integer" />
   <variable name="source_" tag="source" type="source_xml" />