crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Clean up tally triggers

Browse source
This commit is contained in:
Sterling Harper 2019-02-02 12:46:43 -05:00
parent 3ea876235e
commit 9fe05ad9c8
5 changed files with 42 additions and 64 deletions
Download patch file Download diff file Expand all files Collapse all files

5
include/openmc/constants.h
View file

@ -385,11 +385,6 @@ constexpr int IN_BOTTOM {10}; // z min
constexpr int OUT_TOP {11}; // z max
constexpr int IN_TOP {12}; // z max
// Tally trigger types and threshold
constexpr int VARIANCE {1};
constexpr int RELATIVE_ERROR {2};
constexpr int STANDARD_DEVIATION {3};
// Global tally parameters
constexpr int N_GLOBAL_TALLIES {4};
constexpr int K_COLLISION {0};

16
include/openmc/tallies/trigger.h
View file

@ -8,29 +8,27 @@
namespace openmc {
//==============================================================================
// Structs
// Type definitions
//==============================================================================
enum class TriggerMetric {
variance, relative_error, standard_deviation, not_active
};
//! Stops the simulation early if a desired tally uncertainty is reached.
struct Trigger
{
Trigger(int type_, double threshold_, int score_index_)
: type(type_), threshold(threshold_), score_index(score_index_) {}
int type; //!< variance, std_dev, or rel_err
TriggerMetric metric;
double threshold; //!< uncertainty value below which trigger is satisfied
int score_index; //!< index of the relevant score in the tally's arrays
double variance {0.};
double std_dev {0.};
double rel_err {0.};
};
//! Stops the simulation early if a desired k-effective uncertainty is reached.
struct KTrigger
{
int type {0};
TriggerMetric metric {TriggerMetric::not_active};
double threshold {0.};
};

6
src/settings.cpp
View file

@ -155,11 +155,11 @@ void get_run_parameters(pugi::xml_node node_base)
if (check_for_node(node_keff_trigger, "type")) {
auto temp = get_node_value(node_keff_trigger, "type", true, true);
if (temp == "std_dev") {
keff_trigger.type = STANDARD_DEVIATION;
keff_trigger.metric = TriggerMetric::standard_deviation;
} else if (temp == "variance") {
keff_trigger.type = VARIANCE;
keff_trigger.metric = TriggerMetric::variance;
} else if (temp == "rel_err") {
keff_trigger.type = RELATIVE_ERROR;
keff_trigger.metric = TriggerMetric::relative_error;
} else {
fatal_error("Unrecognized keff trigger type " + temp);
}

15
src/tallies/tally.cpp
View file

@ -119,15 +119,15 @@ Tally::init_triggers(pugi::xml_node node, int i_tally)
for (auto trigger_node: node.children("trigger")) {
// Read the trigger type.
int trigger_type;
TriggerMetric metric;
if (check_for_node(trigger_node, "type")) {
auto type_str = get_node_value(trigger_node, "type");
if (type_str == "std_dev") {
trigger_type = STANDARD_DEVIATION;
metric = TriggerMetric::standard_deviation;
} else if (type_str == "variance") {
trigger_type = VARIANCE;
metric = TriggerMetric::variance;
} else if (type_str == "rel_err") {
trigger_type = RELATIVE_ERROR;
metric = TriggerMetric::relative_error;
} else {
std::stringstream msg;
msg << "Unknown trigger type \"" << type_str << "\" in tally " << id_;
@ -170,8 +170,7 @@ Tally::init_triggers(pugi::xml_node node, int i_tally)
if (score_str == "all") {
triggers_.reserve(triggers_.size() + n_tally_scores);
for (auto i_score = 0; i_score < n_tally_scores; ++i_score) {
//TODO: off-by-one
triggers_.emplace_back(trigger_type, threshold, i_score+1);
triggers_.push_back({metric, threshold, i_score});
}
} else {
int i_score = 0;
@ -184,9 +183,7 @@ Tally::init_triggers(pugi::xml_node node, int i_tally)
<< id_ << " but it was listed in a trigger on that tally";
fatal_error(msg);
}
//TODO: off-by-one
triggers_.emplace_back(trigger_type, threshold, i_score+1);
std::cout << i_score+1 << "\n";
triggers_.push_back({metric, threshold, i_score});
}
}
}

64
src/tallies/trigger.cpp
View file

@ -33,9 +33,8 @@ get_tally_uncertainty(int i_tally, int score_index, int filter_index)
int err = openmc_tally_get_n_realizations(i_tally, &n);
auto results = tally_results(i_tally);
//TODO: off-by-one
auto sum = results(filter_index-1, score_index-1, RESULT_SUM);
auto sum_sq = results(filter_index-1, score_index-1, RESULT_SUM_SQ);
auto sum = results(filter_index, score_index, RESULT_SUM);
auto sum_sq = results(filter_index, score_index, RESULT_SUM_SQ);
auto mean = sum / n;
double std_dev = std::sqrt((sum_sq/n - mean*mean) / (n-1));
@ -57,62 +56,51 @@ check_tally_triggers(double& ratio, int& tally_id, int& score)
ratio = 0.;
//TODO: off-by-one
for (auto i_tally = 1; i_tally < model::tallies.size()+1; ++i_tally) {
//TODO: can I mike trigger and t const?
Tally& t {*model::tallies[i_tally-1]};
const Tally& t {*model::tallies[i_tally-1]};
// Ignore tallies with less than two realizations.
int n_reals;
int err = openmc_tally_get_n_realizations(i_tally, &n_reals);
if (n_reals < 2) continue;
for (auto& trigger : t.triggers_) {
trigger.std_dev = 0.;
trigger.rel_err = 0.;
trigger.variance = 0.;
for (const auto& trigger : t.triggers_) {
const auto& results = tally_results(i_tally);
//TODO: off-by-one
for (auto filter_index = 1; filter_index < results.shape()[0]+1;
for (auto filter_index = 0; filter_index < results.shape()[0];
++filter_index) {
//TODO: off-by-one
for (auto score_index = 1; score_index < results.shape()[1]+1;
++score_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;
if (trigger.std_dev < std_dev) {
trigger.std_dev = std_dev;
trigger.variance = std_dev * std_dev;
}
if (trigger.rel_err < rel_err) {
trigger.rel_err = rel_err;
}
// Pick out the relevant uncertainty metric for this trigger.
double uncertainty;
switch (trigger.type) {
case VARIANCE:
uncertainty = trigger.variance;
switch (trigger.metric) {
case TriggerMetric::variance:
uncertainty = std_dev * std_dev;
break;
case STANDARD_DEVIATION:
uncertainty = trigger.std_dev;
case TriggerMetric::standard_deviation:
uncertainty = std_dev;
break;
case RELATIVE_ERROR:
uncertainty = trigger.rel_err;
case TriggerMetric::relative_error:
uncertainty = rel_err;
}
// Compute the uncertainty / threshold ratio.
double this_ratio = uncertainty / trigger.threshold;
if (trigger.type == VARIANCE) {
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;
int* scores;
int junk;
err = openmc_tally_get_scores(i_tally, &scores, &junk);
//TODO: off-by-one
score = scores[trigger.score_index-1];
score = scores[trigger.score_index];
err = openmc_tally_get_id(i_tally, &tally_id);
}
}
@ -127,25 +115,25 @@ double
check_keff_trigger()
{
if (settings::run_mode != RUN_MODE_EIGENVALUE) return 0.;
if (settings::keff_trigger.type == 0) return 0.;
if (settings::keff_trigger.metric == TriggerMetric::not_active) return 0.;
double k_combined[2];
int err = openmc_get_keff(k_combined);
double uncertainty = 0.;
switch (settings::keff_trigger.type) {
case VARIANCE:
switch (settings::keff_trigger.metric) {
case TriggerMetric::variance:
uncertainty = k_combined[1] * k_combined[1];
break;
case STANDARD_DEVIATION:
case TriggerMetric::standard_deviation:
uncertainty = k_combined[1];
break;
case RELATIVE_ERROR:
case TriggerMetric::relative_error:
uncertainty = k_combined[1] / k_combined[0];
}
double ratio = uncertainty / settings::keff_trigger.threshold;
if (settings::keff_trigger.type == VARIANCE)
if (settings::keff_trigger.metric == TriggerMetric::variance)
ratio = std::sqrt(ratio);
return ratio;
}