22920521397

Committed 10 Mar 2026 07:28PM UTC coverage: 81.221% (-0.5%) from 81.721%

Build # 22920521397

Build Type

Pull #3810

github

Committed by

web-flow

Commit Message

Merge c81063b52 into 1dc4aa988

Pull Request Pull Request #3810: Implementation of migration area score

Coverage Stats

17031 of 24387 branches covered (69.84%)

Branch coverage included in aggregate %.

37 of 52 new or added lines in 5 files covered. (71.15%)

2537 existing lines in 86 files now uncovered.

57019 of 66784 relevant lines covered (85.38%)

35780397.68 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(	308✔
31	int i_tally, int score_index, int filter_index)
32	{
33	const auto& tally {model::tallies[i_tally]};	308✔
34
35	auto sum = tally->results_(filter_index, score_index, TallyResult::SUM);	308✔
36	auto sum_sq = tally->results_(filter_index, score_index, TallyResult::SUM_SQ);	308✔
37
38	int n = tally->n_realizations_;	308✔
39	auto mean = sum / n;	308✔
40
41	// if the result has no contributions, return an invalid pair
42	if (mean == 0)	308✔
43	return {-1, -1};	91✔
44
45	double std_dev = std::sqrt((sum_sq / n - mean * mean) / (n - 1));	217✔
46	double rel_err = (mean != 0.) ? std_dev / std::abs(mean) : 0.;	217✔
47
48	return {std_dev, rel_err};	217✔
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)	385✔
59	{
60	ratio = 0.;	385✔
61	for (auto i_tally = 0; i_tally < model::tallies.size(); ++i_tally) {	805✔
62	const Tally& t {*model::tallies[i_tally]};	504!
63
64	// Ignore tallies with less than two realizations.
65	if (t.n_realizations_ < 2)	504!
66	continue;	×
67
68	for (const auto& trigger : t.triggers_) {	693✔
69	// Skip trigger if it is not active
70	if (trigger.metric == TriggerMetric::not_active)	273!
71	continue;	×
72
73	const auto& results = t.results_;
74	for (auto filter_index = 0; filter_index < results.shape(0);	994!
75	++filter_index) {
76	// Compute the tally uncertainty metrics.
77	auto uncert_pair =	308✔
78	get_tally_uncertainty(i_tally, trigger.score_index, filter_index);	308✔
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) {	308✔
83	ratio = INFINITY;	84✔
84	score = t.scores_[trigger.score_index];	84✔
85	tally_id = t.id_;	84✔
86	return;	84✔
87	}
88
89	double std_dev = uncert_pair.first;	224✔
90	double rel_err = uncert_pair.second;	224✔
91
92	// Pick out the relevant uncertainty metric for this trigger.
93	double uncertainty;	224✔
94	switch (trigger.metric) {	224!
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:	224✔
102	uncertainty = rel_err;	224✔
103	break;	224✔
104	case TriggerMetric::not_active:
105	UNREACHABLE();
106	}
107
108	// Compute the uncertainty / threshold ratio.
109	double this_ratio = uncertainty / trigger.threshold;	224✔
110	if (trigger.metric == TriggerMetric::variance) {	224!
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) {	224✔
116	ratio = this_ratio;	175✔
117	score = t.scores_[trigger.score_index];	175✔
118	tally_id = t.id_;	175✔
119	}
120	}
121	}
122	}
123	}
124
125	//! Compute the uncertainty/threshold ratio for the eigenvalue trigger.
126
127	double check_keff_trigger()	385✔
128	{
129	if (settings::run_mode != RunMode::EIGENVALUE)	385✔
130	return 0.0;
131
132	double k_combined[2];	364✔
133	openmc_get_keff(k_combined);	364✔
134
135	double uncertainty = 0.;	364✔
136	switch (settings::keff_trigger.metric) {	364!
137	case TriggerMetric::variance:	×
138	uncertainty = k_combined[1] * k_combined[1];	×
139	break;	×
140	case TriggerMetric::standard_deviation:	280✔
141	uncertainty = k_combined[1];	280✔
142	break;	280✔
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;	280✔
153	if (settings::keff_trigger.metric == TriggerMetric::variance)	280!
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()	87,048✔
161	{
162	// Make some aliases.
163	const auto current_batch {simulation::current_batch};	87,048✔
164	const auto n_batches {settings::n_batches};	87,048✔
165	const auto interval {settings::trigger_batch_interval};	87,048✔
166
167	// See if the current batch is one for which the triggers must be checked.
168	if (!settings::trigger_on)	87,048✔
169	return;	86,712✔
170	if (current_batch < n_batches)	1,484✔
171	return;
172	if (((current_batch - n_batches) % interval) != 0)	735✔
173	return;
174
175	// Check the eigenvalue and tally triggers.
176	double keff_ratio = check_keff_trigger();	385✔
177	double tally_ratio;	385✔
178	int tally_id, score;	385✔
179	check_tally_triggers(tally_ratio, tally_id, score);	385✔
180
181	// If all the triggers are satisfied, alert the user and return.
182	if (std::max(keff_ratio, tally_ratio) <= 1.) {	490✔
183	simulation::satisfy_triggers = true;	49✔
184	write_message(7, "Triggers satisfied for batch {}", current_batch);	49✔
185	return;	49✔
186	}
187
188	// At least one trigger is unsatisfied. Let the user know which one.
189	simulation::satisfy_triggers = false;	336✔
190	std::string msg;	336✔
191	if (keff_ratio >= tally_ratio) {	336✔
192	msg = fmt::format("Triggers unsatisfied, max unc./thresh. is {} for "	504✔
193	"eigenvalue",
194	keff_ratio);	252✔
195	} else {
196	if (tally_ratio == INFINITY) {	84!
197	msg = fmt::format(	168✔
198	"Triggers unsatisfied, no result tallied for score {} in tally {}",
199	reaction_name(score), tally_id);	252✔
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);	336✔
207
208	// Estimate batches til triggers are satisfied.
209	if (settings::trigger_predict) {	336✔
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);	63!
213	auto n_active = current_batch - settings::n_inactive;	63✔
214	auto n_pred_batches = static_cast<int>(n_active * max_ratio * max_ratio) +	63✔
215	settings::n_inactive + 1;	63✔
216
217	if (max_ratio == INFINITY) {	63!
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 =	63✔
224	fmt::format("The estimated number of batches is {}", n_pred_batches);	63✔
225	if (n_pred_batches > settings::n_max_batches) {	63✔
226	msg.append(" --- greater than max batches");	21✔
227	warning(msg);	21✔
228	} else {
229	write_message(msg, 7);	42✔
230	}
231	}	63✔
232	}
233	}	336✔
234
235	} // namespace openmc

openmc-dev / openmc / 22920521397

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