Merge pull request #2344 from pshriwise/fix_triggers

Fix triggers

src/tallies/trigger.cpp

@@ -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.;
 
@@ -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);
 

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
+