Translate print_results to C++

include/openmc/math_functions.h

@@ -34,7 +34,7 @@ extern "C" double normal_percentile(double p);
 //! \return The requested percentile
 //==============================================================================
 
-extern "C" double t_percentile_c(double p, int df);
+extern "C" double t_percentile(double p, int df);
 
 //==============================================================================
 //! Calculate the n-th order Legendre polynomials at the value of x.

include/openmc/output.h

@@ -5,6 +5,7 @@
 #define OPENMC_OUTPUT_H
 
 #include <string>
+#include <utility> // for pair
 
 #include "openmc/particle.h"
 
@@ -65,5 +66,13 @@ void print_batch_keff();
 
 void print_runtime();
 
+//! Display results for global tallies including k-effective estimators
+
+void print_results();
+
+//! Calculate the mean and standard deviation for a tally result
+
+std::pair<double, double> mean_stdev(const double* x, int n);
+
 } // namespace openmc
 #endif // OPENMC_OUTPUT_H

src/eigenvalue.cpp

@@ -331,7 +331,7 @@ void calculate_average_keff()
       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);
+        t_value = t_percentile(1.0 - alpha/2.0, n - 1);
       } else {
         t_value = 1.0;
       }

src/math.F90

@@ -11,8 +11,7 @@ module math
 
   interface
 
-    pure function t_percentile(p, df) bind(C, name='t_percentile_c') &
-         result(t)
+    pure function t_percentile(p, df) bind(C) result(t)
       use ISO_C_BINDING
       implicit none
       real(C_DOUBLE), value, intent(in) :: p

src/math_functions.cpp

@@ -60,7 +60,7 @@ double normal_percentile(double p) {
 }
 
 
-double t_percentile_c(double p, int df){
+double t_percentile(double p, int df){
   double t;
 
   if (df == 1) {

src/output.F90

@@ -93,73 +93,6 @@ contains
 
   end subroutine header
 
-!===============================================================================
-! PRINT_RESULTS displays various estimates of k-effective as well as the global
-! leakage rate.
-!===============================================================================
-
-  subroutine print_results() bind(C)
-
-    integer :: n       ! number of realizations
-    real(8) :: alpha   ! significance level for CI
-    real(8) :: t_n1    ! t-value with N-1 degrees of freedom
-    real(8) :: t_n3    ! t-value with N-3 degrees of freedom
-    real(8) :: x(2)    ! mean and standard deviation
-    real(C_DOUBLE) :: k_combined(2)
-    integer(C_INT) :: err
-
-    ! display header block for results
-    call header("Results", 4)
-
-    n = n_realizations
-
-    if (confidence_intervals) then
-      ! Calculate t-value for confidence intervals
-      alpha = ONE - CONFIDENCE_LEVEL
-      t_n1 = t_percentile(ONE - alpha/TWO, n - 1)
-      t_n3 = t_percentile(ONE - alpha/TWO, n - 3)
-    else
-      t_n1 = ONE
-      t_n3 = ONE
-    end if
-
-    ! write global tallies
-    if (n > 1) then
-      associate (r => global_tallies(RESULT_SUM:RESULT_SUM_SQ, :))
-        if (run_mode == MODE_EIGENVALUE) then
-          x(:) = mean_stdev(r(:, K_COLLISION), n)
-          write(ou,102) "k-effective (Collision)", x(1), t_n1 * x(2)
-          x(:) = mean_stdev(r(:, K_TRACKLENGTH), n)
-          write(ou,102) "k-effective (Track-length)", x(1), t_n1 * x(2)
-          x(:) = mean_stdev(r(:, K_ABSORPTION), n)
-          write(ou,102) "k-effective (Absorption)", x(1), t_n1 * x(2)
-          if (n > 3) then
-            err = openmc_get_keff(k_combined)
-            write(ou,102) "Combined k-effective", k_combined(1), &
-                 t_n3 * k_combined(2)
-          end if
-        end if
-        x(:) = mean_stdev(r(:, LEAKAGE), n)
-        write(ou,102) "Leakage Fraction", x(1), t_n1 * x(2)
-      end associate
-    else
-      if (master) call warning("Could not compute uncertainties -- only one &
-           &active batch simulated!")
-
-      if (run_mode == MODE_EIGENVALUE) then
-        write(ou,103) "k-effective (Collision)", global_tallies(RESULT_SUM, K_COLLISION) / n
-        write(ou,103) "k-effective (Track-length)", global_tallies(RESULT_SUM, K_TRACKLENGTH) / n
-        write(ou,103) "k-effective (Absorption)", global_tallies(RESULT_SUM, K_ABSORPTION) / n
-      end if
-      write(ou,103) "Leakage Fraction", global_tallies(RESULT_SUM, LEAKAGE) / n
-    end if
-    write(ou,*)
-
-102 format (1X,A,T30,"= ",F8.5," +/- ",F8.5)
-103 format (1X,A,T30,"= ",F8.5)
-
-  end subroutine print_results
-
 !===============================================================================
 ! WRITE_TALLIES creates an output file and writes out the mean values of all
 ! tallies and their standard deviations

src/output.cpp

@@ -9,19 +9,23 @@
 #include <ctime>
 
 #include <omp.h>
+#include "xtensor/xview.hpp"
 
 #include "openmc/capi.h"
 #include "openmc/cell.h"
 #include "openmc/constants.h"
 #include "openmc/eigenvalue.h"
+#include "openmc/error.h"
 #include "openmc/geometry.h"
 #include "openmc/lattice.h"
+#include "openmc/math_functions.h"
 #include "openmc/message_passing.h"
 #include "openmc/plot.h"
 #include "openmc/settings.h"
 #include "openmc/simulation.h"
 #include "openmc/surface.h"
 #include "openmc/timer.h"
+#include "openmc/tallies/tally.h"
 
 namespace openmc {
 
@@ -102,7 +106,7 @@ header(const char* msg, int level) {
 
   // Print header based on verbosity level.
   if (settings::verbosity >= level) {
-    std::cout << out.str() << "\n\n";
+    std::cout << '\n' << out.str() << "\n\n";
   }
 }
 
@@ -493,4 +497,84 @@ void print_runtime()
   std::cout.flags(f);
 }
 
+//==============================================================================
+
+void print_results()
+{
+  // Save state of cout
+  auto f {std::cout.flags()};
+
+  // display header block for results
+  header("Results", 4);
+
+  // Calculate t-value for confidence intervals
+  int n = n_realizations;
+  double alpha, t_n1, t_n3;
+  if (settings::confidence_intervals) {
+    alpha = 1.0 - CONFIDENCE_LEVEL;
+    t_n1 = t_percentile(1.0 - alpha/2.0, n - 1);
+    t_n3 = t_percentile(1.0 - alpha/2.0, n - 3);
+  } else {
+    t_n1 = 1.0;
+    t_n3 = 1.0;
+  }
+
+  // Set formatting for floats
+  std::cout << std::fixed << std::setprecision(5);
+
+  // write global tallies
+  auto gt = global_tallies();
+  double mean, stdev;
+  if (n > 1) {
+    if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+      std::tie(mean, stdev) = mean_stdev(&gt(K_COLLISION, 0), n);
+      std::cout << " k-effective (Collision)     = "
+        << mean << " +/- " << t_n1 * stdev << '\n';
+      std::tie(mean, stdev) = mean_stdev(&gt(K_TRACKLENGTH, 0), n);
+      std::cout << " k-effective (Track-length)  = "
+        << mean << " +/- " << t_n1 * stdev << '\n';
+      std::tie(mean, stdev) = mean_stdev(&gt(K_ABSORPTION, 0), n);
+      std::cout << " k-effective (Absorption)    = "
+        << mean << " +/- " << t_n1 * stdev << '\n';
+      if (n > 3) {
+        double k_combined[2];
+        openmc_get_keff(k_combined);
+        std::cout << " Combined k-effective        = "
+          << k_combined[0] << " +/- " << t_n3 * k_combined[1] << '\n';
+      }
+    }
+    std::tie(mean, stdev) = mean_stdev(&gt(LEAKAGE, 0), n);
+    std::cout << " Leakage Fraction            = "
+      << mean << " +/- " << t_n1 * stdev << '\n';
+  } else {
+    if (mpi::master) warning("Could not compute uncertainties -- only one "
+      "active batch simulated!");
+
+    if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+      std::cout << " k-effective (Collision)    = "
+        << gt(K_COLLISION, RESULT_SUM) / n << '\n';
+      std::cout << " k-effective (Track-length) = "
+        << gt(K_TRACKLENGTH, RESULT_SUM) / n << '\n';
+      std::cout << " k-effective (Absorption)   = "
+        << gt(K_ABSORPTION, RESULT_SUM) / n << '\n';
+    }
+    std::cout << " Leakage Fraction           = "
+      << gt(LEAKAGE, RESULT_SUM) / n << '\n';
+  }
+  std::cout << '\n';
+
+  // Restore state of cout
+  std::cout.flags(f);
+}
+
+//==============================================================================
+
+std::pair<double, double> mean_stdev(const double* x, int n)
+{
+  double mean = x[RESULT_SUM] / n;
+  double stdev = n > 1 ? std::sqrt((x[RESULT_SUM_SQ]/n
+    - mean*mean)/(n - 1)) : 0.0;
+  return {mean, stdev};
+}
+
 } // namespace openmc

src/simulation.cpp

@@ -31,7 +31,6 @@ extern "C" void allocate_tally_results();
 extern "C" void check_triggers();
 extern "C" void init_tally_routines();
 extern "C" void load_state_point();
-extern "C" void print_results();
 extern "C" void setup_active_tallies();
 extern "C" void simulation_init_f();
 extern "C" void simulation_finalize_f();