Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
93 changes: 54 additions & 39 deletions src/tallies/trigger.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,10 @@ std::pair<double, double> 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.;

Expand Down Expand Up @@ -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;

// 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_;
}
// Compute the tally uncertainty metrics.
auto uncert_pair =
get_tally_uncertainty(i_tally, trigger.score_index, filter_index);

// 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;
}

double std_dev = uncert_pair.first;
double rel_err = uncert_pair.second;

// 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_;
}
}
}
Expand Down Expand Up @@ -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);

Expand Down
75 changes: 75 additions & 0 deletions tests/unit_tests/test_triggers.py
Original file line number Diff line number Diff line change
@@ -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