diff --git a/src/simulation.cpp b/src/simulation.cpp index 3591953cb..4565a6391 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -16,6 +16,7 @@ #include "openmc/timer.h" #include "openmc/tallies/filter.h" #include "openmc/tallies/tally.h" +#include "openmc/tallies/trigger.h" #include @@ -27,7 +28,6 @@ namespace openmc { // data/functions from Fortran side extern "C" void accumulate_tallies(); 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_batch_keff(); diff --git a/src/tallies/trigger.F90 b/src/tallies/trigger.F90 index 4aea52326..d65b14c8f 100644 --- a/src/tallies/trigger.F90 +++ b/src/tallies/trigger.F90 @@ -2,106 +2,13 @@ module trigger use, intrinsic :: ISO_C_BINDING - use constants - use eigenvalue, only: openmc_get_keff - use endf, only: reaction_name - use error, only: warning, write_message - use string, only: to_str - use message_passing, only: master - use settings - use simulation_header use trigger_header - use tally, only: TallyObject - use tally_header, only: tallies, n_tallies + use tally_header, only: tallies implicit none - interface - 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 - contains -!=============================================================================== -! CHECK_TRIGGERS checks any user-specified precision triggers' for convergence -! and predicts the number of remainining batches to convergence. -!=============================================================================== - - subroutine check_triggers() bind(C) - - implicit none - - real(8) :: keff_ratio ! uncertainty/threshold ratio for keff - real(8) :: tally_ratio ! max tally uncertainty/threshold ratio - integer :: tally_id ! id for tally with max ratio - integer :: score ! tally score with max ratio - integer :: n_pred_batches ! predicted # batches to satisfy all triggers - - if (.not. master) return - - ! Checks if current_batch is one for which the triggers must be checked - if (current_batch < n_batches .or. (.not. trigger_on)) return - if (mod((current_batch - n_batches), n_batch_interval) /= 0 .and. & - current_batch /= n_max_batches) return - - ! Check the eigenvalue and tally triggers - 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 " // & - trim(to_str(current_batch)), 7) - return - end if - - ! Specify which trigger is unsatisfied - if (keff_ratio >= tally_ratio) then - call write_message("Triggers unsatisfied, max unc./thresh. is " // & - trim(to_str(keff_ratio)) // " for eigenvalue", 7) - else - call write_message("Triggers unsatisfied, max unc./thresh. is " // & - trim(to_str(tally_ratio)) // " for " // trim(reaction_name(score)) & - // " in tally " // trim(to_str(tally_id)), 7) - end if - - ! Estimate batches till triggers are satisfied - if (pred_batches) then - ! Estimate the number of remaining batches to convergence - ! The prediction uses the fact that tally variances are proportional - ! to 1/N where N is the number of the batches/particles - n_batch_interval = int((current_batch-n_inactive) * & - (max(keff_ratio, tally_ratio) ** 2)) + n_inactive-n_batches + 1 - n_pred_batches = n_batch_interval + n_batches - - ! Write the predicted number of batches for the user - if (n_pred_batches > n_max_batches) then - call warning("The estimated number of batches is " // & - trim(to_str(n_pred_batches)) // & - " -- greater than max batches. ") - else - call write_message("The estimated number of batches is " // & - trim(to_str(n_pred_batches)), 7) - end if - end if - end subroutine check_triggers - function n_tally_triggers(i_tally) bind(C) result(n) integer(C_INT), value :: i_tally integer(C_INT) :: n diff --git a/src/tallies/trigger.cpp b/src/tallies/trigger.cpp index 6a6fe231d..79926e934 100644 --- a/src/tallies/trigger.cpp +++ b/src/tallies/trigger.cpp @@ -1,12 +1,15 @@ #include "openmc/tallies/trigger.h" #include +#include #include // for std::pair #include "openmc/capi.h" #include "openmc/constants.h" -#include "openmc/message_passing.h" +#include "openmc/error.h" +#include "openmc/reaction.h" #include "openmc/settings.h" +#include "openmc/simulation.h" #include "openmc/tallies/tally.h" namespace openmc { @@ -40,10 +43,10 @@ extern "C" Trigger* get_tally_trigger(int i_tally, int i_trig); //! 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) +void +check_tally_triggers(double& ratio, int& tally_id, int& score) { - *ratio = 0.; + 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]}; @@ -94,17 +97,17 @@ check_tally_triggers(double* ratio, int* tally_id, int* score) double this_ratio = uncertainty / trigger.threshold; if (trigger.type == VARIANCE) { - this_ratio = std::sqrt(*ratio); + this_ratio = std::sqrt(ratio); } - if (this_ratio > *ratio) { - *ratio = this_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); + score = scores[trigger.score_index-1]; + err = openmc_tally_get_id(i_tally, &tally_id); } } } @@ -114,7 +117,7 @@ check_tally_triggers(double* ratio, int* tally_id, int* score) //! Computes the uncertainty/threshold ratio for the eigenvalue trigger. -extern "C" double +double check_keff_trigger() { if (settings::run_mode != RUN_MODE_EIGENVALUE) return 0.; @@ -141,4 +144,65 @@ check_keff_trigger() return ratio; } +//! See if tally and eigenvalue uncertainties are under trigger thresholds. + +void +check_triggers() +{ + // Make some aliases. + const auto current_batch {simulation::current_batch}; + const auto n_batches {settings::n_batches}; + const auto interval {settings::trigger_batch_interval}; + + // See if the current batch is one for which the triggers must be checked. + if (!settings::trigger_on) return; + if (current_batch < n_batches) return; + if (((current_batch - n_batches) % interval) != 0) return; + + // Check the eigenvalue and tally triggers. + double keff_ratio = check_keff_trigger(); + double tally_ratio; + int tally_id, score; + check_tally_triggers(tally_ratio, tally_id, score); + + // If all the triggers are satisfied, alert the user and return. + if (std::max(keff_ratio, tally_ratio) <= 1.) { + simulation::satisfy_triggers = true; + write_message("Triggers satisfied for batch " + + std::to_string(current_batch), 7); + return; + } + + // At least one trigger is unsatisfied. Let the user know which one. + simulation::satisfy_triggers = false; + std::stringstream msg; + msg << "Triggers unsatisfied, max unc./thresh. is "; + if (keff_ratio >= tally_ratio) { + msg << keff_ratio << " for eigenvalue"; + } else { + msg << tally_ratio << " for " << reaction_name(score) << " in tally " + << tally_id; + } + write_message(msg, 7); + + // Estimate batches til triggers are satisfied. + if (settings::trigger_predict) { + // This calculation assumes tally variance is proportional to 1/N where N is + // the number of batches. + auto max_ratio = std::max(keff_ratio, tally_ratio); + auto n_active = current_batch - settings::n_inactive; + auto n_pred_batches = static_cast(n_active * max_ratio * max_ratio) + + settings::n_inactive + 1; + + std::stringstream msg; + msg << "The estimated number of batches is " << n_pred_batches; + if (n_pred_batches > settings::n_max_batches) { + msg << " --- greater than max batches"; + warning(msg); + } else { + write_message(msg, 7); + } + } +} + } // namespace openmc