diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 7fc7dcedec..b55ea061f5 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -385,11 +385,6 @@ constexpr int IN_BOTTOM {10}; // z min constexpr int OUT_TOP {11}; // z max constexpr int IN_TOP {12}; // z max -// Tally trigger types and threshold -constexpr int VARIANCE {1}; -constexpr int RELATIVE_ERROR {2}; -constexpr int STANDARD_DEVIATION {3}; - // Global tally parameters constexpr int N_GLOBAL_TALLIES {4}; constexpr int K_COLLISION {0}; diff --git a/include/openmc/tallies/trigger.h b/include/openmc/tallies/trigger.h index 458c2925a0..f4f44e20ad 100644 --- a/include/openmc/tallies/trigger.h +++ b/include/openmc/tallies/trigger.h @@ -8,29 +8,27 @@ namespace openmc { //============================================================================== -// Structs +// Type definitions //============================================================================== +enum class TriggerMetric { + variance, relative_error, standard_deviation, not_active +}; + //! Stops the simulation early if a desired tally uncertainty is reached. struct Trigger { - Trigger(int type_, double threshold_, int score_index_) - : type(type_), threshold(threshold_), score_index(score_index_) {} - - int type; //!< variance, std_dev, or rel_err + TriggerMetric metric; double threshold; //!< uncertainty value below which trigger is satisfied int score_index; //!< index of the relevant score in the tally's arrays - double variance {0.}; - double std_dev {0.}; - double rel_err {0.}; }; //! Stops the simulation early if a desired k-effective uncertainty is reached. struct KTrigger { - int type {0}; + TriggerMetric metric {TriggerMetric::not_active}; double threshold {0.}; }; diff --git a/src/settings.cpp b/src/settings.cpp index a7543d5e61..0457ea13ff 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -155,11 +155,11 @@ void get_run_parameters(pugi::xml_node node_base) if (check_for_node(node_keff_trigger, "type")) { auto temp = get_node_value(node_keff_trigger, "type", true, true); if (temp == "std_dev") { - keff_trigger.type = STANDARD_DEVIATION; + keff_trigger.metric = TriggerMetric::standard_deviation; } else if (temp == "variance") { - keff_trigger.type = VARIANCE; + keff_trigger.metric = TriggerMetric::variance; } else if (temp == "rel_err") { - keff_trigger.type = RELATIVE_ERROR; + keff_trigger.metric = TriggerMetric::relative_error; } else { fatal_error("Unrecognized keff trigger type " + temp); } diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index f217efd110..81d84456a4 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -119,15 +119,15 @@ Tally::init_triggers(pugi::xml_node node, int i_tally) for (auto trigger_node: node.children("trigger")) { // Read the trigger type. - int trigger_type; + TriggerMetric metric; if (check_for_node(trigger_node, "type")) { auto type_str = get_node_value(trigger_node, "type"); if (type_str == "std_dev") { - trigger_type = STANDARD_DEVIATION; + metric = TriggerMetric::standard_deviation; } else if (type_str == "variance") { - trigger_type = VARIANCE; + metric = TriggerMetric::variance; } else if (type_str == "rel_err") { - trigger_type = RELATIVE_ERROR; + metric = TriggerMetric::relative_error; } else { std::stringstream msg; msg << "Unknown trigger type \"" << type_str << "\" in tally " << id_; @@ -170,8 +170,7 @@ Tally::init_triggers(pugi::xml_node node, int i_tally) if (score_str == "all") { triggers_.reserve(triggers_.size() + n_tally_scores); for (auto i_score = 0; i_score < n_tally_scores; ++i_score) { - //TODO: off-by-one - triggers_.emplace_back(trigger_type, threshold, i_score+1); + triggers_.push_back({metric, threshold, i_score}); } } else { int i_score = 0; @@ -184,9 +183,7 @@ Tally::init_triggers(pugi::xml_node node, int i_tally) << id_ << " but it was listed in a trigger on that tally"; fatal_error(msg); } - //TODO: off-by-one - triggers_.emplace_back(trigger_type, threshold, i_score+1); - std::cout << i_score+1 << "\n"; + triggers_.push_back({metric, threshold, i_score}); } } } diff --git a/src/tallies/trigger.cpp b/src/tallies/trigger.cpp index 0bc783b963..d8f3ed4b7a 100644 --- a/src/tallies/trigger.cpp +++ b/src/tallies/trigger.cpp @@ -33,9 +33,8 @@ get_tally_uncertainty(int i_tally, int score_index, int filter_index) int err = openmc_tally_get_n_realizations(i_tally, &n); auto results = tally_results(i_tally); - //TODO: off-by-one - auto sum = results(filter_index-1, score_index-1, RESULT_SUM); - auto sum_sq = results(filter_index-1, score_index-1, RESULT_SUM_SQ); + auto sum = results(filter_index, score_index, RESULT_SUM); + auto sum_sq = results(filter_index, score_index, RESULT_SUM_SQ); auto mean = sum / n; double std_dev = std::sqrt((sum_sq/n - mean*mean) / (n-1)); @@ -57,62 +56,51 @@ 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) { - //TODO: can I mike trigger and t const? - Tally& t {*model::tallies[i_tally-1]}; + 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; - for (auto& trigger : t.triggers_) { - trigger.std_dev = 0.; - trigger.rel_err = 0.; - trigger.variance = 0.; - + for (const auto& trigger : t.triggers_) { const auto& results = tally_results(i_tally); - //TODO: off-by-one - for (auto filter_index = 1; filter_index < results.shape()[0]+1; + for (auto filter_index = 0; filter_index < results.shape()[0]; ++filter_index) { - //TODO: off-by-one - for (auto score_index = 1; score_index < results.shape()[1]+1; - ++score_index) { + for (auto score_index = 0; score_index < results.shape()[1]; + ++score_index) { + // Compute the tally uncertainty metrics. 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; - } + // Pick out the relevant uncertainty metric for this trigger. double uncertainty; - switch (trigger.type) { - case VARIANCE: - uncertainty = trigger.variance; + switch (trigger.metric) { + case TriggerMetric::variance: + uncertainty = std_dev * std_dev; break; - case STANDARD_DEVIATION: - uncertainty = trigger.std_dev; + case TriggerMetric::standard_deviation: + uncertainty = std_dev; break; - case RELATIVE_ERROR: - uncertainty = trigger.rel_err; + case TriggerMetric::relative_error: + uncertainty = rel_err; } + // Compute the uncertainty / threshold ratio. double this_ratio = uncertainty / trigger.threshold; - if (trigger.type == VARIANCE) { + if (trigger.metric == TriggerMetric::variance) { this_ratio = std::sqrt(ratio); } + // If this is the most uncertain value, set the output variables. 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]; + score = scores[trigger.score_index]; err = openmc_tally_get_id(i_tally, &tally_id); } } @@ -127,25 +115,25 @@ double check_keff_trigger() { if (settings::run_mode != RUN_MODE_EIGENVALUE) return 0.; - if (settings::keff_trigger.type == 0) return 0.; + if (settings::keff_trigger.metric == TriggerMetric::not_active) return 0.; double k_combined[2]; int err = openmc_get_keff(k_combined); double uncertainty = 0.; - switch (settings::keff_trigger.type) { - case VARIANCE: + switch (settings::keff_trigger.metric) { + case TriggerMetric::variance: uncertainty = k_combined[1] * k_combined[1]; break; - case STANDARD_DEVIATION: + case TriggerMetric::standard_deviation: uncertainty = k_combined[1]; break; - case RELATIVE_ERROR: + case TriggerMetric::relative_error: uncertainty = k_combined[1] / k_combined[0]; } double ratio = uncertainty / settings::keff_trigger.threshold; - if (settings::keff_trigger.type == VARIANCE) + if (settings::keff_trigger.metric == TriggerMetric::variance) ratio = std::sqrt(ratio); return ratio; }