22500302709

Committed 27 Feb 2026 07:16PM UTC coverage: 81.512% (-0.3%) from 81.826%

Build # 22500302709

Build Type

Pull #3830

github

Committed by

web-flow

Commit Message

Merge 25fbb4266 into b3788f11e

Pull Request Pull Request #3830: Parallelize sampling external sources and threadsafe rejection counters

Run Details

17488 of 25193 branches covered (69.42%)

Branch coverage included in aggregate %.

59 of 66 new or added lines in 6 files covered. (89.39%)

841 existing lines in 44 files now uncovered.

57726 of 67081 relevant lines covered (86.05%)

44920080.48 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

76.41

/src/tallies/trigger.cpp

#include "openmc/tallies/trigger.h"

#include <cmath>
#include <utility> // for std::pair

#include <fmt/core.h>

#include "openmc/capi.h"
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/reaction.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/tallies/tally.h"

namespace openmc {

//==============================================================================
// Global variable definitions
//==============================================================================

namespace settings {
KTrigger keff_trigger;
}

//==============================================================================
// Non-member functions
//==============================================================================

std::pair<double, double> get_tally_uncertainty(
  int i_tally, int score_index, int filter_index)
{
  const auto& tally {model::tallies[i_tally]};

  auto sum = tally->results_(filter_index, score_index, TallyResult::SUM);
  auto sum_sq = tally->results_(filter_index, score_index, TallyResult::SUM_SQ);

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

  return {std_dev, rel_err};
}

//! Find the limiting limiting tally trigger.
//
//! param[out] ratio The uncertainty/threshold ratio for the most limiting
//!   tally trigger
//! param[out] tally_id The ID number of the most limiting tally
//! param[out] score The most limiting tally score bin

void check_tally_triggers(double& ratio, int& tally_id, int& score)
{
  ratio = 0.;
  for (auto i_tally = 0; i_tally < model::tallies.size(); ++i_tally) {
    const Tally& t {*model::tallies[i_tally]};

    // Ignore tallies with less than two realizations.
    if (t.n_realizations_ < 2)
      continue;

    for (const auto& trigger : t.triggers_) {
      // Skip trigger if it is not active
      if (trigger.metric == TriggerMetric::not_active)
        continue;

      const auto& results = t.results_;
      for (auto filter_index = 0; filter_index < results.shape(0);
           ++filter_index) {
        // 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, unless zero scores are ignored for this trigger.
        if (uncert_pair.first == -1 && !trigger.ignore_zeros) {
          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_;
        }
      }
    }
  }
}

//! Compute the uncertainty/threshold ratio for the eigenvalue trigger.

double check_keff_trigger()
{
  if (settings::run_mode != RunMode::EIGENVALUE)
    return 0.0;

  double k_combined[2];
  openmc_get_keff(k_combined);

  double uncertainty = 0.;
  switch (settings::keff_trigger.metric) {
  case TriggerMetric::variance:
    uncertainty = k_combined[1] * k_combined[1];
    break;
  case TriggerMetric::standard_deviation:
    uncertainty = k_combined[1];
    break;
  case TriggerMetric::relative_error:
    uncertainty = k_combined[1] / k_combined[0];
    break;
  default:
    // If it's an unrecognized TriggerMetric or no keff trigger is on,
    // return 0 to stop division by zero where "ratio" is calculated.
    return 0.0;
  }

  double ratio = uncertainty / settings::keff_trigger.threshold;
  if (settings::keff_trigger.metric == TriggerMetric::variance)
    ratio = std::sqrt(ratio);
  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(7, "Triggers satisfied for batch {}", current_batch);
    return;
  }

  // At least one trigger is unsatisfied.  Let the user know which one.
  simulation::satisfy_triggers = false;
  std::string msg;
  if (keff_ratio >= tally_ratio) {
    msg = fmt::format("Triggers unsatisfied, max unc./thresh. is {} for "
                      "eigenvalue",
      keff_ratio);
  } else {
    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);

  // 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<int>(n_active * max_ratio * max_ratio) +
                          settings::n_inactive + 1;

    if (max_ratio == INFINITY) {
      std::string msg =
        fmt::format("One or more tallies with triggers have no scores. Unable "
                    "to estimate the number of remaining batches.");
      write_message(msg, 7);
    } else {
      std::string msg =
        fmt::format("The estimated number of batches is {}", n_pred_batches);
      if (n_pred_batches > settings::n_max_batches) {
        msg.append(" --- greater than max batches");
        warning(msg);
      } else {
        write_message(msg, 7);
      }
    }
  }
}

} // namespace openmc

1	#include "openmc/tallies/trigger.h"
2
3	#include <cmath>
4	#include <utility> // for std::pair
5
6	#include <fmt/core.h>
7
8	#include "openmc/capi.h"
9	#include "openmc/constants.h"
10	#include "openmc/error.h"
11	#include "openmc/reaction.h"
12	#include "openmc/settings.h"
13	#include "openmc/simulation.h"
14	#include "openmc/tallies/tally.h"
15
16	namespace openmc {
17
18	//==============================================================================
19	// Global variable definitions
20	//==============================================================================
21
22	namespace settings {
23	KTrigger keff_trigger;
24	}
25
26	//==============================================================================
27	// Non-member functions
28	//==============================================================================
29
30	std::pair<double, double> get_tally_uncertainty(	484✔
31	int i_tally, int score_index, int filter_index)
32	{
33	const auto& tally {model::tallies[i_tally]};	484✔
34
35	auto sum = tally->results_(filter_index, score_index, TallyResult::SUM);	484✔
36	auto sum_sq = tally->results_(filter_index, score_index, TallyResult::SUM_SQ);	484✔
37
38	int n = tally->n_realizations_;	484✔
39	auto mean = sum / n;	484✔
40
41	// if the result has no contributions, return an invalid pair
42	if (mean == 0)	484✔
43	return {-1, -1};	143✔
44
45	double std_dev = std::sqrt((sum_sq / n - mean * mean) / (n - 1));	341✔
46	double rel_err = (mean != 0.) ? std_dev / std::abs(mean) : 0.;	341✔
47
48	return {std_dev, rel_err};	341✔
49	}
50
51	//! Find the limiting limiting tally trigger.
52	//
53	//! param[out] ratio The uncertainty/threshold ratio for the most limiting
54	//! tally trigger
55	//! param[out] tally_id The ID number of the most limiting tally
56	//! param[out] score The most limiting tally score bin
57
58	void check_tally_triggers(double& ratio, int& tally_id, int& score)	605✔
59	{
60	ratio = 0.;	605✔
61	for (auto i_tally = 0; i_tally < model::tallies.size(); ++i_tally) {	1,265✔
62	const Tally& t {*model::tallies[i_tally]};	792!
63
64	// Ignore tallies with less than two realizations.
65	if (t.n_realizations_ < 2)	792!
66	continue;	×
67
68	for (const auto& trigger : t.triggers_) {	1,089✔
69	// Skip trigger if it is not active
70	if (trigger.metric == TriggerMetric::not_active)	429!
71	continue;	×
72
73	const auto& results = t.results_;
74	for (auto filter_index = 0; filter_index < results.shape(0);	1,562!
75	++filter_index) {
76	// Compute the tally uncertainty metrics.
77	auto uncert_pair =	484✔
78	get_tally_uncertainty(i_tally, trigger.score_index, filter_index);	484✔
79
80	// If there is a score without contributions, set ratio to inf and
81	// exit early, unless zero scores are ignored for this trigger.
82	if (uncert_pair.first == -1 && !trigger.ignore_zeros) {	484✔
83	ratio = INFINITY;	132✔
84	score = t.scores_[trigger.score_index];	132✔
85	tally_id = t.id_;	132✔
86	return;	132✔
87	}
88
89	double std_dev = uncert_pair.first;	352✔
90	double rel_err = uncert_pair.second;	352✔
91
92	// Pick out the relevant uncertainty metric for this trigger.
93	double uncertainty;	352✔
94	switch (trigger.metric) {	352!
95	case TriggerMetric::variance:	×
96	uncertainty = std_dev * std_dev;	×
97	break;	×
98	case TriggerMetric::standard_deviation:	×
99	uncertainty = std_dev;	×
100	break;	×
101	case TriggerMetric::relative_error:	352✔
102	uncertainty = rel_err;	352✔
103	break;	352✔
104	case TriggerMetric::not_active:
105	UNREACHABLE();
106	}
107
108	// Compute the uncertainty / threshold ratio.
109	double this_ratio = uncertainty / trigger.threshold;	352✔
110	if (trigger.metric == TriggerMetric::variance) {	352!
111	this_ratio = std::sqrt(ratio);	×
112	}
113
114	// If this is the most uncertain value, set the output variables.
115	if (this_ratio > ratio) {	352✔
116	ratio = this_ratio;	275✔
117	score = t.scores_[trigger.score_index];	275✔
118	tally_id = t.id_;	275✔
119	}
120	}
121	}
122	}
123	}
124
125	//! Compute the uncertainty/threshold ratio for the eigenvalue trigger.
126
127	double check_keff_trigger()	605✔
128	{
129	if (settings::run_mode != RunMode::EIGENVALUE)	605✔
130	return 0.0;
131
132	double k_combined[2];	572✔
133	openmc_get_keff(k_combined);	572✔
134
135	double uncertainty = 0.;	572✔
136	switch (settings::keff_trigger.metric) {	572!
137	case TriggerMetric::variance:	×
138	uncertainty = k_combined[1] * k_combined[1];	×
139	break;	×
140	case TriggerMetric::standard_deviation:	440✔
141	uncertainty = k_combined[1];	440✔
142	break;	440✔
143	case TriggerMetric::relative_error:	×
144	uncertainty = k_combined[1] / k_combined[0];	×
145	break;	×
146	default:
147	// If it's an unrecognized TriggerMetric or no keff trigger is on,
148	// return 0 to stop division by zero where "ratio" is calculated.
149	return 0.0;
150	}
151
152	double ratio = uncertainty / settings::keff_trigger.threshold;	440✔
153	if (settings::keff_trigger.metric == TriggerMetric::variance)	440!
154	ratio = std::sqrt(ratio);	×
155	return ratio;
156	}
157
158	//! See if tally and eigenvalue uncertainties are under trigger thresholds.
159
160	void check_triggers()	136,865✔
161	{
162	// Make some aliases.
163	const auto current_batch {simulation::current_batch};	136,865✔
164	const auto n_batches {settings::n_batches};	136,865✔
165	const auto interval {settings::trigger_batch_interval};	136,865✔
166
167	// See if the current batch is one for which the triggers must be checked.
168	if (!settings::trigger_on)	136,865✔
169	return;	136,337✔
170	if (current_batch < n_batches)	2,332✔
171	return;
172	if (((current_batch - n_batches) % interval) != 0)	1,155✔
173	return;
174
175	// Check the eigenvalue and tally triggers.
176	double keff_ratio = check_keff_trigger();	605✔
177	double tally_ratio;	605✔
178	int tally_id, score;	605✔
179	check_tally_triggers(tally_ratio, tally_id, score);	605✔
180
181	// If all the triggers are satisfied, alert the user and return.
182	if (std::max(keff_ratio, tally_ratio) <= 1.) {	770✔
183	simulation::satisfy_triggers = true;	77✔
184	write_message(7, "Triggers satisfied for batch {}", current_batch);	77✔
185	return;	77✔
186	}
187
188	// At least one trigger is unsatisfied. Let the user know which one.
189	simulation::satisfy_triggers = false;	528✔
190	std::string msg;	528✔
191	if (keff_ratio >= tally_ratio) {	528✔
192	msg = fmt::format("Triggers unsatisfied, max unc./thresh. is {} for "	792✔
193	"eigenvalue",
194	keff_ratio);	396✔
195	} else {
196	if (tally_ratio == INFINITY) {	132!
197	msg = fmt::format(	264✔
198	"Triggers unsatisfied, no result tallied for score {} in tally {}",
199	reaction_name(score), tally_id);	396✔
200	} else {
201	msg = fmt::format(	×
202	"Triggers unsatisfied, max unc./thresh. is {} for {} in tally {}",
203	tally_ratio, reaction_name(score), tally_id);	×
204	}
205	}
206	write_message(msg, 7);	528✔
207
208	// Estimate batches til triggers are satisfied.
209	if (settings::trigger_predict) {	528✔
210	// This calculation assumes tally variance is proportional to 1/N where N is
211	// the number of batches.
212	auto max_ratio = std::max(keff_ratio, tally_ratio);	99!
213	auto n_active = current_batch - settings::n_inactive;	99✔
214	auto n_pred_batches = static_cast<int>(n_active * max_ratio * max_ratio) +	99✔
215	settings::n_inactive + 1;	99✔
216
217	if (max_ratio == INFINITY) {	99!
UNCOV 218	std::string msg =	×
219	fmt::format("One or more tallies with triggers have no scores. Unable "
220	"to estimate the number of remaining batches.");	×
221	write_message(msg, 7);	×
222	} else {	×
223	std::string msg =	99✔
224	fmt::format("The estimated number of batches is {}", n_pred_batches);	99✔
225	if (n_pred_batches > settings::n_max_batches) {	99✔
226	msg.append(" --- greater than max batches");	33✔
227	warning(msg);	33✔
228	} else {
229	write_message(msg, 7);	66✔
230	}
231	}	99✔
232	}
233	}	528✔
234
235	} // namespace openmc

openmc-dev / openmc / 22500302709

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous