Calculation of kinetic parameters relevant to fissile system with an external neutron source

include/openmc/constants.h

@@ -324,8 +324,15 @@ enum TallyScore {
 };
 
 // Global tally parameters
-constexpr int N_GLOBAL_TALLIES {4};
-enum class GlobalTally { K_COLLISION, K_ABSORPTION, K_TRACKLENGTH, LEAKAGE };
+enum class GlobalTally {
+  K_COLLISION,
+  K_ABSORPTION,
+  K_TRACKLENGTH,
+  LEAKAGE,
+  K_TRACKLENGTH_SQ, // for calculation of stddev for generational k
+  SIZE
+};
+constexpr int N_GLOBAL_TALLIES {static_cast<int>(GlobalTally::SIZE)};
 
 // Miscellaneous
 constexpr int C_NONE {-1};

include/openmc/eigenvalue.h

@@ -21,7 +21,8 @@ namespace openmc {
 
 namespace simulation {
 
-extern double keff_generation; //!<  Single-generation k on each processor
+extern array<double, 2>
+  keff_generation;             //!<  Single-generation k on each processor
 extern array<double, 2> k_sum; //!< Used to reduce sum and sum_sq
 extern vector<double> entropy; //!< Shannon entropy at each generation
 extern xt::xtensor<double, 1> source_frac; //!< Source fraction for UFS

include/openmc/simulation.h

@@ -4,6 +4,7 @@
 #ifndef OPENMC_SIMULATION_H
 #define OPENMC_SIMULATION_H
 
+#include "openmc/array.h"
 #include "openmc/mesh.h"
 #include "openmc/particle.h"
 #include "openmc/vector.h"
@@ -45,7 +46,7 @@ extern int64_t work_per_rank;      //!< number of particles per MPI rank
 extern const RegularMesh* entropy_mesh;
 extern const RegularMesh* ufs_mesh;
 
-extern vector<double> k_generation;
+extern vector<array<double, 2>> k_generation;
 extern vector<int64_t> work_index;
 
 } // namespace simulation

include/openmc/state_point.h

@@ -49,6 +49,8 @@ void write_source_bank(hid_t group_id, span<SourceSite> source_bank,
 
 void read_source_bank(
   hid_t group_id, vector<SourceSite>& sites, bool distribute);
+void write_global_tallies(hid_t file_id);
+void read_global_tallies(hid_t file_id);
 void write_tally_results_nr(hid_t file_id);
 void restart_set_keff();
 void write_unstructured_mesh_results();

include/openmc/tallies/tally.h

@@ -225,6 +225,7 @@ extern "C" int32_t n_realizations;
 extern double global_tally_absorption;
 extern double global_tally_collision;
 extern double global_tally_tracklength;
+extern double global_tally_tracklength_sq;
 extern double global_tally_leakage;
 
 //==============================================================================

openmc/lib/tally.py

@@ -127,11 +127,11 @@ def global_tallies():
     """
     ptr = POINTER(c_double)()
     _dll.openmc_global_tallies(ptr)
-    array = as_array(ptr, (4, 3))
+    array = as_array(ptr, (5, 3))
 
     # Get sum, sum-of-squares, and number of realizations
-    sum_ = array[:, 1]
-    sum_sq = array[:, 2]
+    sum_ = array[:4, 1]
+    sum_sq = array[:4, 2]
     n = num_realizations()
 
     # Determine mean

openmc/statepoint.py

@@ -8,6 +8,7 @@
 import numpy as np
 from pathlib import Path
 from uncertainties import ufloat
+from uncertainties.unumpy import uarray
 
 import openmc
 import openmc.checkvalue as cv
@@ -80,6 +81,8 @@ class StatePoint:
         .. versionadded:: 0.13.1
     meshes : dict
         Dictionary whose keys are mesh IDs and whose values are MeshBase objects
+    multiplication : uncertainties.UFloat
+        Estimator for multiplication of initial source neutrons 
     n_batches : int
         Number of batches
     n_inactive : int
@@ -106,6 +109,8 @@ class StatePoint:
         'E', 'wgt', 'delayed_group', 'surf_id', and 'particle', corresponding to
         the position, direction, energy, weight, delayed group, surface ID and
         particle type of the source site, respectively.
+    source_effectiveness : uncertainties.UFloat
+        Source effectiveness factor.        
     source_present : bool
         Indicate whether source sites are present
     sparse : bool
@@ -242,8 +247,8 @@ def global_tallies(self):
                 ('mean', 'f8'), ('std_dev', 'f8')])
             gt['name'] = ['k-collision', 'k-absorption', 'k-tracklength',
                           'leakage']
-            gt['sum'] = data[:,1]
-            gt['sum_sq'] = data[:,2]
+            gt['sum'] = data[:,0]
+            gt['sum_sq'] = data[:,1]
 
             # Calculate mean and sample standard deviation of mean
             n = self.n_realizations
@@ -264,7 +269,13 @@ def k_cmfd(self):
     @property
     def k_generation(self):
         if self.run_mode == 'eigenvalue':
-            return self._f['k_generation'][()]
+            arr = self._f['k_generation'][()]
+            if arr.ndim==1:
+                return arr
+            elif arr.ndim==2:
+                return uarray(arr[:,0],arr[:,1])
+            else:
+                raise ValueError(f'k_generation shape ({arr.shape}) must be either 1d or 2d')
         else:
             return None
 
@@ -317,6 +328,17 @@ def meshes(self):
             self._meshes_read = True
 
         return self._meshes
+
+    @property
+    def multiplication(self):
+        if self.run_mode == 'eigenvalue':
+            k_gen = self.k_generation
+            k_inactive = k_gen[:self.n_inactive]
+            cumprod = np.cumprod(k_inactive)
+            cumprod[-1] *= 1/(1-self.keff)
+            return 1.0+np.sum(cumprod)
+        else:
+            return None        
 
     @property
     def n_batches(self):
@@ -365,6 +387,15 @@ def stride(self):
     @property
     def source(self):
         return self._f['source_bank'][()] if self.source_present else None
+
+    @property
+    def source_efficiency(self):
+        if self.run_mode == 'eigenvalue':
+            k_gen = self.k_generation
+            k_inactive = k_gen[:self.n_inactive]
+            return np.prod(k_inactive/self.keff)
+        else:
+            return None          
 
     @property
     def source_present(self):

src/eigenvalue.cpp

@@ -36,7 +36,7 @@ namespace openmc {
 
 namespace simulation {
 
-double keff_generation;
+array<double, 2> keff_generation;
 array<double, 2> k_sum;
 vector<double> entropy;
 xt::xtensor<double, 1> source_frac;
@@ -52,16 +52,21 @@ void calculate_generation_keff()
   const auto& gt = simulation::global_tallies;
 
   // Get keff for this generation by subtracting off the starting value
-  simulation::keff_generation =
+  simulation::keff_generation[0] =
     gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) -
-    simulation::keff_generation;
+    simulation::keff_generation[0];
+  simulation::keff_generation[1] =
+    gt(GlobalTally::K_TRACKLENGTH_SQ, TallyResult::VALUE) -
+    simulation::keff_generation[1];
 
-  double keff_reduced;
+  array<double, 2> keff_reduced;
 #ifdef OPENMC_MPI
   if (settings::solver_type != SolverType::RANDOM_RAY) {
     // Combine values across all processors
-    MPI_Allreduce(&simulation::keff_generation, &keff_reduced, 1, MPI_DOUBLE,
-      MPI_SUM, mpi::intracomm);
+    MPI_Allreduce(&simulation::keff_generation[0], &keff_reduced[0], 1,
+      MPI_DOUBLE, MPI_SUM, mpi::intracomm);
+    MPI_Allreduce(&simulation::keff_generation[1], &keff_reduced[1], 1,
+      MPI_DOUBLE, MPI_SUM, mpi::intracomm);
   } else {
     // If using random ray, MPI parallelism is provided by domain replication.
     // As such, all fluxes will be reduced at the end of each transport sweep,
@@ -77,10 +82,13 @@ void calculate_generation_keff()
   // Normalize single batch estimate of k
   // TODO: This should be normalized by total_weight, not by n_particles
   if (settings::solver_type != SolverType::RANDOM_RAY) {
-    keff_reduced /= settings::n_particles;
+    keff_reduced[0] /= settings::n_particles;
+    keff_reduced[1] /= settings::n_particles;
   }
-
-  simulation::k_generation.push_back(keff_reduced);
+  double k_mean = keff_reduced[0];
+  double k_std = std::sqrt(
+    (keff_reduced[1] - std::pow(k_mean, 2)) / (settings::n_particles - 1));
+  simulation::k_generation.push_back({k_mean, k_std});
 }
 
 void synchronize_bank()
@@ -379,11 +387,11 @@ void calculate_average_keff()
   if (n <= 0) {
     // For inactive generations, use current generation k as estimate for next
     // generation
-    simulation::keff = simulation::k_generation[i];
+    simulation::keff = simulation::k_generation[i][0];
   } else {
     // Sample mean of keff
-    simulation::k_sum[0] += simulation::k_generation[i];
-    simulation::k_sum[1] += std::pow(simulation::k_generation[i], 2);
+    simulation::k_sum[0] += simulation::k_generation[i][0];
+    simulation::k_sum[1] += std::pow(simulation::k_generation[i][0], 2);
 
     // Determine mean
     simulation::keff = simulation::k_sum[0] / n;
@@ -660,7 +668,13 @@ void write_eigenvalue_hdf5(hid_t group)
 {
   write_dataset(group, "n_inactive", settings::n_inactive);
   write_dataset(group, "generations_per_batch", settings::gen_per_batch);
-  write_dataset(group, "k_generation", simulation::k_generation);
+  auto n = simulation::k_generation.size();
+  xt::xtensor<double, 2> k_generation({n, 2});
+  for (int i = 0; i < n; ++i) {
+    k_generation(i, 0) = simulation::k_generation[i][0];
+    k_generation(i, 1) = simulation::k_generation[i][1];
+  }
+  write_dataset(group, "k_generation", k_generation);
   if (settings::entropy_on) {
     write_dataset(group, "entropy", simulation::entropy);
   }
@@ -675,9 +689,13 @@ void write_eigenvalue_hdf5(hid_t group)
 void read_eigenvalue_hdf5(hid_t group)
 {
   read_dataset(group, "generations_per_batch", settings::gen_per_batch);
-  int n = simulation::restart_batch * settings::gen_per_batch;
+  size_t n = simulation::restart_batch * settings::gen_per_batch;
+  xt::xtensor<double, 2> k_generation({n, 2});
+  read_dataset(group, "k_generation", k_generation);
   simulation::k_generation.resize(n);
-  read_dataset(group, "k_generation", simulation::k_generation);
+  for (int i = 0; i < n; ++i) {
+    simulation::k_generation[i] = {k_generation(i, 0), k_generation(i, 1)};
+  }
   if (settings::entropy_on) {
     read_dataset(group, "entropy", simulation::entropy);
   }

src/output.cpp

@@ -19,6 +19,7 @@
 #endif
 #include "xtensor/xview.hpp"
 
+#include "openmc/array.h"
 #include "openmc/capi.h"
 #include "openmc/cell.h"
 #include "openmc/constants.h"
@@ -370,11 +371,12 @@ void print_build_info()
 void print_columns()
 {
   if (settings::entropy_on) {
-    fmt::print("  Bat./Gen.      k       Entropy         Average k \n"
-               "  =========   ========   ========   ====================\n");
+    fmt::print(
+      "  Bat./Gen.             k             Entropy         Average k \n"
+      "  =========   ====================   ========   ====================\n");
   } else {
-    fmt::print("  Bat./Gen.      k            Average k\n"
-               "  =========   ========   ====================\n");
+    fmt::print("  Bat./Gen.             k                  Average k\n"
+               "  =========   ====================   ====================\n");
   }
 }
 
@@ -392,7 +394,10 @@ void print_generation()
   // write out batch/generation and generation k-effective
   auto batch_and_gen = std::to_string(simulation::current_batch) + "/" +
                        std::to_string(simulation::current_gen);
-  fmt::print("  {:>9}   {:8.5f}", batch_and_gen, simulation::k_generation[idx]);
+  array<double, 2> k_gen;
+  k_gen[0] = simulation::k_generation[idx][0];
+  k_gen[1] = simulation::k_generation[idx][1];
+  fmt::print("  {:>9}   {:8.5f} +/-{:8.5f}", batch_and_gen, k_gen[0], k_gen[1]);
 
   // write out entropy info
   if (settings::entropy_on) {

src/particle.cpp

@@ -478,8 +478,11 @@ void Particle::event_death()
   global_tally_absorption += keff_tally_absorption();
 #pragma omp atomic
   global_tally_collision += keff_tally_collision();
+  double val = keff_tally_tracklength();
 #pragma omp atomic
-  global_tally_tracklength += keff_tally_tracklength();
+  global_tally_tracklength += val;
+#pragma omp atomic
+  global_tally_tracklength_sq += val * val;
 #pragma omp atomic
   global_tally_leakage += keff_tally_leakage();
 

src/random_ray/random_ray_simulation.cpp

@@ -498,6 +498,7 @@ void RandomRaySimulation::simulate()
 
         // Store random ray k-eff into OpenMC's native k-eff variable
         global_tally_tracklength = k_eff_;
+        global_tally_tracklength_sq = std::pow(k_eff_, 2);
       }
 
       // Execute all tallying tasks, if this is an active batch

src/simulation.cpp

@@ -1,5 +1,6 @@
 #include "openmc/simulation.h"
 
+#include "openmc/array.h"
 #include "openmc/bank.h"
 #include "openmc/capi.h"
 #include "openmc/container_util.h"
@@ -318,7 +319,7 @@ int64_t work_per_rank;
 const RegularMesh* entropy_mesh {nullptr};
 const RegularMesh* ufs_mesh {nullptr};
 
-vector<double> k_generation;
+vector<array<double, 2>> k_generation;
 vector<int64_t> work_index;
 
 } // namespace simulation
@@ -506,8 +507,10 @@ void initialize_generation()
       ufs_count_sites();
 
     // Store current value of tracklength k
-    simulation::keff_generation = simulation::global_tallies(
-      GlobalTally::K_TRACKLENGTH, TallyResult::VALUE);
+    auto& gt = simulation::global_tallies;
+    simulation::keff_generation = {
+      gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE),
+      gt(GlobalTally::K_TRACKLENGTH_SQ, TallyResult::VALUE)};
   }
 }
 
@@ -522,6 +525,8 @@ void finalize_generation()
       global_tally_absorption;
     gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) +=
       global_tally_tracklength;
+    gt(GlobalTally::K_TRACKLENGTH_SQ, TallyResult::VALUE) +=
+      global_tally_tracklength_sq;
   }
   gt(GlobalTally::LEAKAGE, TallyResult::VALUE) += global_tally_leakage;
 
@@ -530,6 +535,7 @@ void finalize_generation()
     global_tally_collision = 0.0;
     global_tally_absorption = 0.0;
     global_tally_tracklength = 0.0;
+    global_tally_tracklength_sq = 0.0;
   }
   global_tally_leakage = 0.0;