From a732efbe51d29d87a8ed7f655f4122aba3a36dd7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 5 Feb 2019 12:06:26 -0600 Subject: [PATCH] Translate print_results to C++ --- include/openmc/math_functions.h | 2 +- include/openmc/output.h | 9 ++++ src/eigenvalue.cpp | 2 +- src/math.F90 | 3 +- src/math_functions.cpp | 2 +- src/output.F90 | 67 ------------------------- src/output.cpp | 86 ++++++++++++++++++++++++++++++++- src/simulation.cpp | 1 - 8 files changed, 98 insertions(+), 74 deletions(-) diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h index 87710fb6b..e07374df3 100644 --- a/include/openmc/math_functions.h +++ b/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. diff --git a/include/openmc/output.h b/include/openmc/output.h index c80ca1167..0f26caa4f 100644 --- a/include/openmc/output.h +++ b/include/openmc/output.h @@ -5,6 +5,7 @@ #define OPENMC_OUTPUT_H #include +#include // 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 mean_stdev(const double* x, int n); + } // namespace openmc #endif // OPENMC_OUTPUT_H diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index b3b8b850c..2b0f731cb 100644 --- a/src/eigenvalue.cpp +++ b/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; } diff --git a/src/math.F90 b/src/math.F90 index 5607a9556..5a0a582a4 100644 --- a/src/math.F90 +++ b/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 diff --git a/src/math_functions.cpp b/src/math_functions.cpp index fe2c53230..919c305c4 100644 --- a/src/math_functions.cpp +++ b/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) { diff --git a/src/output.F90 b/src/output.F90 index 78e784220..209d588b2 100644 --- a/src/output.F90 +++ b/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 diff --git a/src/output.cpp b/src/output.cpp index 39ebf0dce..567a588fa 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -9,19 +9,23 @@ #include #include +#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(>(K_COLLISION, 0), n); + std::cout << " k-effective (Collision) = " + << mean << " +/- " << t_n1 * stdev << '\n'; + std::tie(mean, stdev) = mean_stdev(>(K_TRACKLENGTH, 0), n); + std::cout << " k-effective (Track-length) = " + << mean << " +/- " << t_n1 * stdev << '\n'; + std::tie(mean, stdev) = mean_stdev(>(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(>(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 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 diff --git a/src/simulation.cpp b/src/simulation.cpp index 3d362890c..cf6000398 100644 --- a/src/simulation.cpp +++ b/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();