From b66a4d28b9ced77b1e1b88c7bbeb2225627841cd Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Thu, 31 Jan 2019 16:18:46 -0500 Subject: [PATCH] Move most tally trigger math to C++ --- src/tallies/trigger.F90 | 176 +++++++--------------------------------- src/tallies/trigger.cpp | 128 +++++++++++++++++++++++++++-- 2 files changed, 148 insertions(+), 156 deletions(-) diff --git a/src/tallies/trigger.F90 b/src/tallies/trigger.F90 index 93068f4f5..4aea52326 100644 --- a/src/tallies/trigger.F90 +++ b/src/tallies/trigger.F90 @@ -17,14 +17,16 @@ module trigger implicit none interface - subroutine get_tally_uncertainty(i_tally, score_index, filter_index, & - std_dev, rel_err) bind(C) - import C_INT, C_DOUBLE - integer(C_INT), value :: i_tally - integer(C_INT), value :: score_index - integer(C_INT), value :: filter_index - real(C_DOUBLE) :: std_dev - real(C_DOUBLE) :: rel_err + function check_keff_trigger() bind(C) result(ratio) + import C_DOUBLE + real(C_DOUBLE) :: ratio + end function + + subroutine check_tally_triggers(ratio, tally_id, score) bind(C) + import C_DOUBLE, C_INT + real(C_DOUBLE) :: ratio + integer(C_INT) :: tally_id + integer(C_INT) :: score end subroutine end interface @@ -52,13 +54,16 @@ contains if (mod((current_batch - n_batches), n_batch_interval) /= 0 .and. & current_batch /= n_max_batches) return - ! By default, assume all triggers are satisfied - satisfy_triggers = .true. - ! Check the eigenvalue and tally triggers - call check_keff_trigger(keff_ratio) + keff_ratio = check_keff_trigger() call check_tally_triggers(tally_ratio, tally_id, score) + if (max(keff_ratio, tally_ratio) <= ONE) then + satisfy_triggers = .true. + else + satisfy_triggers = .false. + end if + ! Alert the user if the triggers are satisfied if (satisfy_triggers) then call write_message("Triggers satisfied for batch " // & @@ -97,142 +102,17 @@ contains end if end subroutine check_triggers -!=============================================================================== -! CHECK_KEFF_TRIGGER computes the uncertainty/threshold ratio for the eigenvalue -! trigger and updates the global satisfy_tiggers variable if the trigger is -! unsatisfied. -!=============================================================================== + function n_tally_triggers(i_tally) bind(C) result(n) + integer(C_INT), value :: i_tally + integer(C_INT) :: n + n = tallies(i_tally) % obj % n_triggers + end function - subroutine check_keff_trigger(ratio) - real(8), intent(out) :: ratio - integer(C_INT) :: err - real(C_DOUBLE) :: k_combined(2) - real(8) :: uncertainty - - ratio = 0 - if (run_mode == MODE_EIGENVALUE) then - if (keff_trigger % trigger_type /= 0) then - err = openmc_get_keff(k_combined) - select case (keff_trigger % trigger_type) - case(VARIANCE) - uncertainty = k_combined(2) ** 2 - case(STANDARD_DEVIATION) - uncertainty = k_combined(2) - case default - uncertainty = k_combined(2) / k_combined(1) - end select - - ! If uncertainty is above threshold, store uncertainty ratio - if (uncertainty > keff_trigger % threshold) then - satisfy_triggers = .false. - if (keff_trigger % trigger_type == VARIANCE) then - ratio = sqrt(uncertainty / keff_trigger % threshold) - else - ratio = uncertainty / keff_trigger % threshold - end if - end if - end if - end if - end subroutine check_keff_trigger - -!=============================================================================== -! CHECK_TALLY_TRIGGERS computes the uncertainty/threshold ratio for all tally -! triggers and updates the global satisfy_tiggers variable if any trigger is -! unsatisfied. -!=============================================================================== - - subroutine check_tally_triggers(max_ratio, tally_id, score) - - ! Variables to reflect distance to trigger convergence criteria - real(8), intent(out) :: max_ratio ! max uncertainty/thresh ratio - integer, intent(out) :: tally_id ! id for tally with max ratio - integer, intent(out) :: score - - integer :: i ! index in tallies array - integer :: j ! index in tally filters - integer :: n ! loop index for nuclides - integer :: s ! loop index for triggers - integer :: filter_index ! index in results array for filters - integer :: score_index ! scoring bin index - integer(C_INT) :: err - real(8) :: uncertainty ! trigger uncertainty - real(8) :: std_dev = ZERO ! trigger standard deviation - real(8) :: rel_err = ZERO ! trigger relative error - real(8) :: ratio ! ratio of the uncertainty/trigger threshold - real(C_DOUBLE) :: k_combined(2) - - ! Initialize tally trigger maximum uncertainty ratio to zero - max_ratio = 0 - - ! Compute uncertainties for all tallies, scores with triggers - TALLY_LOOP: do i = 1, n_tallies - associate (t => tallies(i) % obj) - - ! Cycle through if only one batch has been simumlate - if (t % n_realizations == 1) then - cycle TALLY_LOOP - end if - - TRIGGER_LOOP: do s = 1, t % n_triggers - associate (trigger => t % triggers(s)) - - ! Initialize trigger uncertainties to zero - trigger % std_dev = ZERO - trigger % rel_err = ZERO - trigger % variance = ZERO - - FILTER_LOOP: do filter_index = 1, t % n_filter_bins() - - ! Initialize score index - score_index = trigger % score_index - - ! Initialize score bin index - NUCLIDE_LOOP: do n = 1, t % n_nuclide_bins() - - call get_tally_uncertainty(i, score_index, filter_index, & - std_dev, rel_err) - - if (trigger % variance < variance) then - trigger % variance = std_dev ** 2 - end if - if (trigger % std_dev < std_dev) then - trigger % std_dev = std_dev - end if - if (trigger % rel_err < rel_err) then - trigger % rel_err = rel_err - end if - - select case (trigger % type) - case(VARIANCE) - uncertainty = trigger % variance - case(STANDARD_DEVIATION) - uncertainty = trigger % std_dev - case default - uncertainty = trigger % rel_err - end select - - if (uncertainty > trigger % threshold) then - satisfy_triggers = .false. - - if (trigger % type == VARIANCE) then - ratio = sqrt(uncertainty / trigger % threshold) - else - ratio = uncertainty / trigger % threshold - end if - - if (max_ratio < ratio) then - max_ratio = ratio - score = t % score_bins(trigger % score_index) - tally_id = t % id - end if - end if - end do NUCLIDE_LOOP - if (t % n_filters() == 0) exit FILTER_LOOP - end do FILTER_LOOP - end associate - end do TRIGGER_LOOP - end associate - end do TALLY_LOOP - end subroutine check_tally_triggers + function get_tally_trigger(i_tally, i_trig) bind(C) result(trigger) + integer(C_INT), value :: i_tally + integer(C_INT), value :: i_trig + type(C_PTR) :: trigger + trigger = C_LOC(tallies(i_tally) % obj % triggers(i_trig)) + end function end module trigger diff --git a/src/tallies/trigger.cpp b/src/tallies/trigger.cpp index d3069d9a0..6a6fe231d 100644 --- a/src/tallies/trigger.cpp +++ b/src/tallies/trigger.cpp @@ -1,15 +1,18 @@ #include "openmc/tallies/trigger.h" #include +#include // for std::pair #include "openmc/capi.h" +#include "openmc/constants.h" +#include "openmc/message_passing.h" +#include "openmc/settings.h" #include "openmc/tallies/tally.h" namespace openmc { -extern "C" void -get_tally_uncertainty(int i_tally, int score_index, int filter_index, - double* std_dev, double* rel_err) +static std::pair +get_tally_uncertainty(int i_tally, int score_index, int filter_index) { int n; int err = openmc_tally_get_n_realizations(i_tally, &n); @@ -20,13 +23,122 @@ get_tally_uncertainty(int i_tally, int score_index, int filter_index, auto sum_sq = results(filter_index-1, score_index-1, RESULT_SUM_SQ); auto mean = sum / n; - *std_dev = std::sqrt((sum_sq/n - mean*mean) / (n-1)); + double std_dev = std::sqrt((sum_sq/n - mean*mean) / (n-1)); + double rel_err = (mean != 0.) ? std_dev / std::abs(mean) : 0.; - if (mean > 0) { - *rel_err = *std_dev / mean; - } else { - *rel_err = 0.; + return {std_dev, rel_err}; +} + +// Functions defined in F90. +extern "C" int n_tally_triggers(int i_tally); +extern "C" Trigger* get_tally_trigger(int i_tally, int i_trig); + +//! Finds the limiting limiting tally trigger. +// +//! param[out] ratio The uncertainty/threshold ratio for the most limiting +//! tally trigger +//! param[out] tally_id The ID number of the most limiting tally +//! param[out] score The most limiting tally score bin + +extern "C" void +check_tally_triggers(double* ratio, int* tally_id, int* score) +{ + *ratio = 0.; + //TODO: off-by-one + for (auto i_tally = 1; i_tally < model::tallies.size()+1; ++i_tally) { + const Tally& t {*model::tallies[i_tally-1]}; + + // Ignore tallies with less than two realizations. + int n_reals; + int err = openmc_tally_get_n_realizations(i_tally, &n_reals); + if (n_reals < 2) continue; + + //TODO: off-by-one + for (auto i_trig = 1; i_trig < n_tally_triggers(i_tally)+1; ++i_trig) { + auto& trigger {*get_tally_trigger(i_tally, i_trig)}; + + trigger.std_dev = 0.; + trigger.rel_err = 0.; + trigger.variance = 0.; + + const auto& results = tally_results(i_tally); + //TODO: off-by-one + for (auto filter_index = 1; filter_index < results.shape()[0]+1; + ++filter_index) { + //TODO: off-by-one + for (auto score_index = 1; score_index < results.shape()[1]+1; + ++score_index) { + auto uncert_pair = get_tally_uncertainty(i_tally, score_index, + filter_index); + double std_dev = uncert_pair.first; + double rel_err = uncert_pair.second; + if (trigger.std_dev < std_dev) { + trigger.std_dev = std_dev; + trigger.variance = std_dev * std_dev; + } + if (trigger.rel_err < rel_err) { + trigger.rel_err = rel_err; + } + + double uncertainty; + switch (trigger.type) { + case VARIANCE: + uncertainty = trigger.variance; + break; + case STANDARD_DEVIATION: + uncertainty = trigger.std_dev; + break; + case RELATIVE_ERROR: + uncertainty = trigger.rel_err; + } + + double this_ratio = uncertainty / trigger.threshold; + if (trigger.type == VARIANCE) { + this_ratio = std::sqrt(*ratio); + } + + if (this_ratio > *ratio) { + *ratio = this_ratio; + int* scores; + int junk; + err = openmc_tally_get_scores(i_tally, &scores, &junk); + //TODO: off-by-one + *score = scores[trigger.score_index-1]; + err = openmc_tally_get_id(i_tally, tally_id); + } + } + } + } } } +//! Computes the uncertainty/threshold ratio for the eigenvalue trigger. + +extern "C" double +check_keff_trigger() +{ + if (settings::run_mode != RUN_MODE_EIGENVALUE) return 0.; + if (settings::keff_trigger.type == 0) return 0.; + + double k_combined[2]; + int err = openmc_get_keff(k_combined); + + double uncertainty = 0.; + switch (settings::keff_trigger.type) { + case VARIANCE: + uncertainty = k_combined[1] * k_combined[1]; + break; + case STANDARD_DEVIATION: + uncertainty = k_combined[1]; + break; + case RELATIVE_ERROR: + uncertainty = k_combined[1] / k_combined[0]; + } + + double ratio = uncertainty / settings::keff_trigger.threshold; + if (settings::keff_trigger.type == VARIANCE) + ratio = std::sqrt(ratio); + return ratio; +} + } // namespace openmc