Convert calculate_average_keff to C++

include/openmc/eigenvalue.h

@@ -1,6 +1,7 @@
 #ifndef OPENMC_EIGENVALUE_H
 #define OPENMC_EIGENVALUE_H
 
+#include <array>
 #include <cstdint> // for int64_t
 #include <vector>
 
@@ -15,6 +16,7 @@ namespace openmc {
 //==============================================================================
 
 extern double keff_generation; //!<  Single-generation k on each processor
+extern std::array<double, 2> k_sum; //!< Used to reduce sum and sum_sq
 extern std::vector<double> entropy; //!< Shannon entropy at each generation
 extern xt::xtensor<double, 1> source_frac; //!< Source fraction for UFS
 
@@ -28,6 +30,13 @@ extern "C" int64_t n_bank;
 //! Collect/normalize the tracklength keff from each process
 extern "C" void calculate_generation_keff();
 
+//! Calcaulte mean/standard deviation of keff during active generations
+//!
+//! This function sets the global variables keff and keff_std which represent
+//! the mean and standard deviation of the mean of k-effective over active
+//! generations. It also broadcasts the value from the master process.
+extern "C" void calculate_average_keff();
+
 //! Sample/redistribute source sites from accumulated fission sites
 extern "C" void synchronize_bank();
 

include/openmc/state_point.h

@@ -11,8 +11,8 @@ namespace openmc {
 
 extern "C" void write_source_bank(hid_t group_id, Bank* source_bank);
 extern "C" void read_source_bank(hid_t group_id, Bank* source_bank);
-
 extern "C" void write_tally_results_nr(hid_t file_id);
+extern "C" void restart_set_keff();
 
 } // namespace openmc
 #endif // OPENMC_STATE_POINT_H

src/api.F90

@@ -4,7 +4,7 @@ module openmc_api
 
   use bank_header,     only: openmc_source_bank
   use constants,       only: K_BOLTZMANN
-  use eigenvalue,      only: k_sum, openmc_get_keff
+  use eigenvalue,      only: openmc_get_keff
   use error
   use geometry,        only: find_cell
   use geometry_header
@@ -248,6 +248,10 @@ contains
     integer(C_INT) :: err
 
     integer :: i
+    interface
+      subroutine k_sum_reset() bind(C)
+      end subroutine
+    end interface
 
     if (allocated(tallies)) then
       do i = 1, size(tallies)
@@ -268,7 +272,7 @@ contains
     k_col_abs = ZERO
     k_col_tra = ZERO
     k_abs_tra = ZERO
-    k_sum(:) = ZERO
+    call k_sum_reset()
 
     ! Reset timers
     call time_total % reset()

src/eigenvalue.F90

@@ -2,79 +2,24 @@ module eigenvalue
 
   use, intrinsic :: ISO_C_BINDING
 
-  use algorithm,   only: binary_search
   use constants,   only: ZERO
-  use error,       only: fatal_error, warning
-  use math,        only: t_percentile
   use message_passing
-  use random_lcg,  only: prn, set_particle_seed, advance_prn_seed
   use settings
   use simulation_header
-  use string,      only: to_str
   use tally_header
-  use timer_header
 
   implicit none
 
-  real(8) :: k_sum(2) = ZERO ! Used to reduce sum and sum_sq
-
   interface
     subroutine calculate_generation_keff() bind(C)
     end subroutine
+
+    subroutine calculate_average_keff() bind(C)
+    end subroutine
   end interface
 
 contains
 
-!===============================================================================
-! CALCULATE_AVERAGE_KEFF calculates the mean and standard deviation of the mean
-! of k-effective during active generations and broadcasts the mean to all
-! processors
-!===============================================================================
-
-  subroutine calculate_average_keff()
-
-    integer :: i        ! overall generation within simulation
-    integer :: n        ! number of active generations
-    real(8) :: alpha    ! significance level for CI
-    real(8) :: t_value  ! t-value for confidence intervals
-
-    ! Determine overall generation and number of active generations
-    i = overall_generation()
-    if (current_batch > n_inactive) then
-      n = gen_per_batch*n_realizations + current_gen
-    else
-      n = 0
-    end if
-
-    if (n <= 0) then
-      ! For inactive generations, use current generation k as estimate for next
-      ! generation
-      keff = k_generation(i)
-
-    else
-      ! Sample mean of keff
-      k_sum(1) = k_sum(1) + k_generation(i)
-      k_sum(2) = k_sum(2) + k_generation(i)**2
-
-      ! Determine mean
-      keff = k_sum(1) / n
-
-      if (n > 1) then
-        if (confidence_intervals) then
-          ! Calculate t-value for confidence intervals
-          alpha = ONE - CONFIDENCE_LEVEL
-          t_value = t_percentile(ONE - alpha/TWO, n - 1)
-        else
-          t_value = ONE
-        end if
-
-        ! Standard deviation of the sample mean of k
-        keff_std = t_value * sqrt((k_sum(2)/n - keff**2) / (n - 1))
-      end if
-    end if
-
-  end subroutine calculate_average_keff
-
 !===============================================================================
 ! OPENMC_GET_KEFF calculates a minimum variance estimate of k-effective based on
 ! a linear combination of the collision, absorption, and tracklength

src/eigenvalue.cpp

@@ -8,6 +8,7 @@
 #include "openmc/constants.h"
 #include "openmc/error.h"
 #include "openmc/hdf5_interface.h"
+#include "openmc/math_functions.h"
 #include "openmc/mesh.h"
 #include "openmc/message_passing.h"
 #include "openmc/random_lcg.h"
@@ -18,6 +19,7 @@
 #include "openmc/tallies/tally.h"
 
 #include <algorithm> // for min
+#include <cmath> // for sqrt
 #include <string>
 
 namespace openmc {
@@ -27,6 +29,7 @@ namespace openmc {
 //==============================================================================
 
 double keff_generation;
+std::array<double, 2> k_sum;
 std::vector<double> entropy;
 xt::xtensor<double, 1> source_frac;
 
@@ -298,6 +301,46 @@ void synchronize_bank()
   time_bank.stop();
 }
 
+void calculate_average_keff()
+{
+  // Determine overall generation and number of active generations
+  int i = overall_generation() - 1;
+  int n;
+  if (simulation::current_batch > settings::n_inactive) {
+    n = settings::gen_per_batch*n_realizations + simulation::current_gen;
+  } else {
+    n = 0;
+  }
+
+  if (n <= 0) {
+    // For inactive generations, use current generation k as estimate for next
+    // generation
+    simulation::keff = simulation::k_generation[i];
+  } else {
+    // Sample mean of keff
+    k_sum[0] += simulation::k_generation[i];
+    k_sum[1] += std::pow(simulation::k_generation[i], 2);
+
+    // Determine mean
+    simulation::keff = k_sum[0] / n;
+
+    if (n > 1) {
+      double t_value;
+      if (settings::confidence_intervals) {
+        // Calculate t-value for confidence intervals
+        double alpha = 1.0 - CONFIDENCE_LEVEL;
+        t_value = t_percentile_c(1.0 - alpha/2.0, n - 1);
+      } else {
+        t_value = 1.0;
+      }
+
+      // Standard deviation of the sample mean of k
+      simulation::keff_std = t_value * std::sqrt((k_sum[1]/n -
+        std::pow(simulation::keff, 2)) / (n - 1));
+    }
+  }
+}
+
 void shannon_entropy()
 {
   // Get pointer to entropy mesh
@@ -430,6 +473,10 @@ extern "C" void read_eigenvalue_hdf5(hid_t group)
   read_dataset(group, "k_abs_tra", simulation::k_abs_tra);
 }
 
+//==============================================================================
+// Fortran compatibility
+//==============================================================================
+
 extern "C" double entropy_c(int i)
 {
   return entropy.at(i - 1);
@@ -440,5 +487,6 @@ extern "C" void entropy_clear()
   entropy.clear();
 }
 
+extern "C" void k_sum_reset() { k_sum.fill(0.0); }
 
 } // namespace openmc

src/simulation.F90

@@ -10,8 +10,7 @@ module simulation
   use cmfd_execute,    only: cmfd_init_batch, cmfd_tally_init, execute_cmfd
   use cmfd_header,     only: cmfd_on
   use constants,       only: ZERO
-  use eigenvalue,      only: calculate_average_keff, calculate_generation_keff, &
-                             k_sum
+  use eigenvalue,      only: calculate_average_keff, calculate_generation_keff
 #ifdef _OPENMP
   use eigenvalue,      only: join_bank_from_threads
 #endif

src/state_point.F90

@@ -16,7 +16,7 @@ module state_point
   use bank_header,        only: Bank
   use cmfd_header
   use constants
-  use eigenvalue,         only: openmc_get_keff, k_sum
+  use eigenvalue,         only: openmc_get_keff
   use endf,               only: reaction_name
   use error,              only: fatal_error, warning, write_message
   use hdf5_interface
@@ -488,7 +488,6 @@ contains
   subroutine load_state_point()
 
     integer :: i
-    integer :: n
     integer :: int_array(3)
     integer, allocatable :: array(:)
     integer(C_INT64_T) :: seed
@@ -504,6 +503,9 @@ contains
         import HID_T
         integer(HID_T), value :: group
       end subroutine
+
+      subroutine restart_set_keff() bind(C)
+      end subroutine
     end interface
 
     ! Write message
@@ -612,17 +614,7 @@ contains
 
     ! Set k_sum, keff, and current_batch based on whether restart file is part
     ! of active cycle or inactive cycle
-    if (restart_batch > n_inactive) then
-      do i = n_inactive + 1, restart_batch
-        k_sum(1) = k_sum(1) + k_generation(i)
-        k_sum(2) = k_sum(2) + k_generation(i)**2
-      end do
-      n = gen_per_batch*n_realizations
-      keff = k_sum(1) / n
-    else
-      n = k_generation_size()
-      keff = k_generation(n)
-    end if
+    call restart_set_keff()
     current_batch = restart_batch
 
     ! Check to make sure source bank is present

src/state_point.cpp

@@ -9,6 +9,7 @@
 
 #include "openmc/capi.h"
 #include "openmc/constants.h"
+#include "openmc/eigenvalue.h"
 #include "openmc/error.h"
 #include "openmc/hdf5_interface.h"
 #include "openmc/message_passing.h"
@@ -275,4 +276,18 @@ void write_tally_results_nr(hid_t file_id)
   }
 }
 
+void restart_set_keff()
+{
+  if (simulation::restart_batch > settings::n_inactive) {
+    for (int i = settings::n_inactive; i < simulation::restart_batch; ++i) {
+      k_sum[0] += simulation::k_generation[i];
+      k_sum[1] += std::pow(simulation::k_generation[i], 2);
+    }
+    int n = settings::gen_per_batch*n_realizations;
+    simulation::keff = k_sum[0] / n;
+  } else {
+    simulation::keff = simulation::k_generation.back();
+  }
+}
+
 } // namespace openmc