diff --git a/src/tallies/trigger.cpp b/src/tallies/trigger.cpp index 2c155980e7..79af658793 100644 --- a/src/tallies/trigger.cpp +++ b/src/tallies/trigger.cpp @@ -37,6 +37,10 @@ std::pair get_tally_uncertainty( int n = tally->n_realizations_; auto mean = sum / n; + + // if the result has no contributions, return an invalid pair + if (mean == 0) return {-1 , -1}; + double std_dev = std::sqrt((sum_sq / n - mean * mean) / (n - 1)); double rel_err = (mean != 0.) ? std_dev / std::abs(mean) : 0.; @@ -68,42 +72,49 @@ void check_tally_triggers(double& ratio, int& tally_id, int& score) const auto& results = t.results_; for (auto filter_index = 0; filter_index < results.shape()[0]; ++filter_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; + // Compute the tally uncertainty metrics. + auto uncert_pair = + get_tally_uncertainty(i_tally, trigger.score_index, filter_index); - // Pick out the relevant uncertainty metric for this trigger. - double uncertainty; - switch (trigger.metric) { - case TriggerMetric::variance: - uncertainty = std_dev * std_dev; - break; - case TriggerMetric::standard_deviation: - uncertainty = std_dev; - break; - case TriggerMetric::relative_error: - uncertainty = rel_err; - break; - case TriggerMetric::not_active: - UNREACHABLE(); - } + // if there is a score without contributions, set ratio to inf and + // exit early + if (uncert_pair.first == -1) { + ratio = INFINITY; + score = t.scores_[trigger.score_index]; + tally_id = t.id_; + return; + } - // Compute the uncertainty / threshold ratio. - double this_ratio = uncertainty / trigger.threshold; - if (trigger.metric == TriggerMetric::variance) { - this_ratio = std::sqrt(ratio); - } + double std_dev = uncert_pair.first; + double rel_err = uncert_pair.second; - // If this is the most uncertain value, set the output variables. - if (this_ratio > ratio) { - ratio = this_ratio; - score = t.scores_[trigger.score_index]; - tally_id = t.id_; - } + // Pick out the relevant uncertainty metric for this trigger. + double uncertainty; + switch (trigger.metric) { + case TriggerMetric::variance: + uncertainty = std_dev * std_dev; + break; + case TriggerMetric::standard_deviation: + uncertainty = std_dev; + break; + case TriggerMetric::relative_error: + uncertainty = rel_err; + break; + case TriggerMetric::not_active: + UNREACHABLE(); + } + + // Compute the uncertainty / threshold ratio. + double this_ratio = uncertainty / trigger.threshold; + 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; + score = t.scores_[trigger.score_index]; + tally_id = t.id_; } } } @@ -181,9 +192,13 @@ void check_triggers() "eigenvalue", keff_ratio); } else { - msg = fmt::format( - "Triggers unsatisfied, max unc./thresh. is {} for {} in tally {}", - tally_ratio, reaction_name(score), tally_id); + if (tally_ratio == INFINITY) { + msg = fmt::format("Triggers unsatisfied, no result tallied for score {} in tally {}", reaction_name(score), tally_id); + } else{ + msg = fmt::format( + "Triggers unsatisfied, max unc./thresh. is {} for {} in tally {}", + tally_ratio, reaction_name(score), tally_id); + } } write_message(msg, 7); diff --git a/tests/unit_tests/test_triggers.py b/tests/unit_tests/test_triggers.py new file mode 100644 index 0000000000..4fe6e044ab --- /dev/null +++ b/tests/unit_tests/test_triggers.py @@ -0,0 +1,75 @@ + +import openmc + +def test_tally_trigger(run_in_tmpdir): + pincell = openmc.examples.pwr_pin_cell() + + # create a tally filter on the materials + mat_filter = openmc.MaterialFilter(pincell.materials) + + # create a tally with triggers applied + tally = openmc.Tally() + tally.filters = [mat_filter] + tally.scores = ['scatter'] + + trigger = openmc.Trigger('rel_err', 0.05) + trigger.scores = ['scatter'] + + tally.triggers = [trigger] + + pincell.tallies = [tally] + + pincell.settings.trigger_active = True + pincell.settings.trigger_max_batches = 100 + pincell.settings.trigger_batch_interval = 5 + + sp_file = pincell.run() + with openmc.StatePoint(sp_file) as sp: + expected_realizations = sp.n_realizations + + # adding other scores to the tally should not change the + # number of batches required to satisfy the trigger + tally.scores = ['total', 'absorption', 'scatter'] + + sp_file = pincell.run() + + with openmc.StatePoint(sp_file) as sp: + realizations = sp.n_realizations + + assert realizations == expected_realizations + + +def test_tally_trigger_null_score(run_in_tmpdir): + pincell = openmc.examples.pwr_pin_cell() + + # create a tally filter on the materials + mat_filter = openmc.MaterialFilter(pincell.materials) + + # apply a tally with a score that be tallied in this model + tally = openmc.Tally() + tally.filters = [mat_filter] + tally.scores = ['pair-production'] + + trigger = openmc.Trigger('rel_err', 0.05) + trigger.scores = ['pair-production'] + + tally.triggers = [trigger] + + pincell.tallies = [tally] + + pincell.settings.trigger_active = True + pincell.settings.trigger_max_batches = 50 + pincell.settings.trigger_batch_interval = 5 + + sp_file = pincell.run() + + with openmc.StatePoint(sp_file) as sp: + # verify that the tally mean is zero + tally_out = sp.get_tally(id=tally.id) + assert all(tally_out.mean == 0.0) + + # we expect that this simulation will run + # up to the max allowed batches + total_batches = sp.n_realizations + sp.n_inactive + assert total_batches == pincell.settings.trigger_max_batches +