mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 13:45:36 -04:00
Move most tally trigger math to C++
This commit is contained in:
Sterling Harper
parent
a07ea809f2
commit
b66a4d28b9
2 changed files with 148 additions and 156 deletions
|
|
@ -17,14 +17,16 @@ module trigger
|
|||
implicit none
|
||||
|
||||
interface
|
||||
subroutine get_tally_uncertainty(i_tally, score_index, filter_index, &
|
||||
std_dev, rel_err) bind(C)
|
||||
import C_INT, C_DOUBLE
|
||||
integer(C_INT), value :: i_tally
|
||||
integer(C_INT), value :: score_index
|
||||
integer(C_INT), value :: filter_index
|
||||
real(C_DOUBLE) :: std_dev
|
||||
real(C_DOUBLE) :: rel_err
|
||||
function check_keff_trigger() bind(C) result(ratio)
|
||||
import C_DOUBLE
|
||||
real(C_DOUBLE) :: ratio
|
||||
end function
|
||||
|
||||
subroutine check_tally_triggers(ratio, tally_id, score) bind(C)
|
||||
import C_DOUBLE, C_INT
|
||||
real(C_DOUBLE) :: ratio
|
||||
integer(C_INT) :: tally_id
|
||||
integer(C_INT) :: score
|
||||
end subroutine
|
||||
end interface
|
||||
|
||||
|
|
@ -52,13 +54,16 @@ contains
|
|||
if (mod((current_batch - n_batches), n_batch_interval) /= 0 .and. &
|
||||
current_batch /= n_max_batches) return
|
||||
|
||||
! By default, assume all triggers are satisfied
|
||||
satisfy_triggers = .true.
|
||||
|
||||
! Check the eigenvalue and tally triggers
|
||||
call check_keff_trigger(keff_ratio)
|
||||
keff_ratio = check_keff_trigger()
|
||||
call check_tally_triggers(tally_ratio, tally_id, score)
|
||||
|
||||
if (max(keff_ratio, tally_ratio) <= ONE) then
|
||||
satisfy_triggers = .true.
|
||||
else
|
||||
satisfy_triggers = .false.
|
||||
end if
|
||||
|
||||
! Alert the user if the triggers are satisfied
|
||||
if (satisfy_triggers) then
|
||||
call write_message("Triggers satisfied for batch " // &
|
||||
|
|
@ -97,142 +102,17 @@ contains
|
|||
end if
|
||||
end subroutine check_triggers
|
||||
|
||||
!===============================================================================
|
||||
! CHECK_KEFF_TRIGGER computes the uncertainty/threshold ratio for the eigenvalue
|
||||
! trigger and updates the global satisfy_tiggers variable if the trigger is
|
||||
! unsatisfied.
|
||||
!===============================================================================
|
||||
function n_tally_triggers(i_tally) bind(C) result(n)
|
||||
integer(C_INT), value :: i_tally
|
||||
integer(C_INT) :: n
|
||||
n = tallies(i_tally) % obj % n_triggers
|
||||
end function
|
||||
|
||||
subroutine check_keff_trigger(ratio)
|
||||
real(8), intent(out) :: ratio
|
||||
integer(C_INT) :: err
|
||||
real(C_DOUBLE) :: k_combined(2)
|
||||
real(8) :: uncertainty
|
||||
|
||||
ratio = 0
|
||||
if (run_mode == MODE_EIGENVALUE) then
|
||||
if (keff_trigger % trigger_type /= 0) then
|
||||
err = openmc_get_keff(k_combined)
|
||||
select case (keff_trigger % trigger_type)
|
||||
case(VARIANCE)
|
||||
uncertainty = k_combined(2) ** 2
|
||||
case(STANDARD_DEVIATION)
|
||||
uncertainty = k_combined(2)
|
||||
case default
|
||||
uncertainty = k_combined(2) / k_combined(1)
|
||||
end select
|
||||
|
||||
! If uncertainty is above threshold, store uncertainty ratio
|
||||
if (uncertainty > keff_trigger % threshold) then
|
||||
satisfy_triggers = .false.
|
||||
if (keff_trigger % trigger_type == VARIANCE) then
|
||||
ratio = sqrt(uncertainty / keff_trigger % threshold)
|
||||
else
|
||||
ratio = uncertainty / keff_trigger % threshold
|
||||
end if
|
||||
end if
|
||||
end if
|
||||
end if
|
||||
end subroutine check_keff_trigger
|
||||
|
||||
!===============================================================================
|
||||
! CHECK_TALLY_TRIGGERS computes the uncertainty/threshold ratio for all tally
|
||||
! triggers and updates the global satisfy_tiggers variable if any trigger is
|
||||
! unsatisfied.
|
||||
!===============================================================================
|
||||
|
||||
subroutine check_tally_triggers(max_ratio, tally_id, score)
|
||||
|
||||
! Variables to reflect distance to trigger convergence criteria
|
||||
real(8), intent(out) :: max_ratio ! max uncertainty/thresh ratio
|
||||
integer, intent(out) :: tally_id ! id for tally with max ratio
|
||||
integer, intent(out) :: score
|
||||
|
||||
integer :: i ! index in tallies array
|
||||
integer :: j ! index in tally filters
|
||||
integer :: n ! loop index for nuclides
|
||||
integer :: s ! loop index for triggers
|
||||
integer :: filter_index ! index in results array for filters
|
||||
integer :: score_index ! scoring bin index
|
||||
integer(C_INT) :: err
|
||||
real(8) :: uncertainty ! trigger uncertainty
|
||||
real(8) :: std_dev = ZERO ! trigger standard deviation
|
||||
real(8) :: rel_err = ZERO ! trigger relative error
|
||||
real(8) :: ratio ! ratio of the uncertainty/trigger threshold
|
||||
real(C_DOUBLE) :: k_combined(2)
|
||||
|
||||
! Initialize tally trigger maximum uncertainty ratio to zero
|
||||
max_ratio = 0
|
||||
|
||||
! Compute uncertainties for all tallies, scores with triggers
|
||||
TALLY_LOOP: do i = 1, n_tallies
|
||||
associate (t => tallies(i) % obj)
|
||||
|
||||
! Cycle through if only one batch has been simumlate
|
||||
if (t % n_realizations == 1) then
|
||||
cycle TALLY_LOOP
|
||||
end if
|
||||
|
||||
TRIGGER_LOOP: do s = 1, t % n_triggers
|
||||
associate (trigger => t % triggers(s))
|
||||
|
||||
! Initialize trigger uncertainties to zero
|
||||
trigger % std_dev = ZERO
|
||||
trigger % rel_err = ZERO
|
||||
trigger % variance = ZERO
|
||||
|
||||
FILTER_LOOP: do filter_index = 1, t % n_filter_bins()
|
||||
|
||||
! Initialize score index
|
||||
score_index = trigger % score_index
|
||||
|
||||
! Initialize score bin index
|
||||
NUCLIDE_LOOP: do n = 1, t % n_nuclide_bins()
|
||||
|
||||
call get_tally_uncertainty(i, score_index, filter_index, &
|
||||
std_dev, rel_err)
|
||||
|
||||
if (trigger % variance < variance) then
|
||||
trigger % variance = std_dev ** 2
|
||||
end if
|
||||
if (trigger % std_dev < std_dev) then
|
||||
trigger % std_dev = std_dev
|
||||
end if
|
||||
if (trigger % rel_err < rel_err) then
|
||||
trigger % rel_err = rel_err
|
||||
end if
|
||||
|
||||
select case (trigger % type)
|
||||
case(VARIANCE)
|
||||
uncertainty = trigger % variance
|
||||
case(STANDARD_DEVIATION)
|
||||
uncertainty = trigger % std_dev
|
||||
case default
|
||||
uncertainty = trigger % rel_err
|
||||
end select
|
||||
|
||||
if (uncertainty > trigger % threshold) then
|
||||
satisfy_triggers = .false.
|
||||
|
||||
if (trigger % type == VARIANCE) then
|
||||
ratio = sqrt(uncertainty / trigger % threshold)
|
||||
else
|
||||
ratio = uncertainty / trigger % threshold
|
||||
end if
|
||||
|
||||
if (max_ratio < ratio) then
|
||||
max_ratio = ratio
|
||||
score = t % score_bins(trigger % score_index)
|
||||
tally_id = t % id
|
||||
end if
|
||||
end if
|
||||
end do NUCLIDE_LOOP
|
||||
if (t % n_filters() == 0) exit FILTER_LOOP
|
||||
end do FILTER_LOOP
|
||||
end associate
|
||||
end do TRIGGER_LOOP
|
||||
end associate
|
||||
end do TALLY_LOOP
|
||||
end subroutine check_tally_triggers
|
||||
function get_tally_trigger(i_tally, i_trig) bind(C) result(trigger)
|
||||
integer(C_INT), value :: i_tally
|
||||
integer(C_INT), value :: i_trig
|
||||
type(C_PTR) :: trigger
|
||||
trigger = C_LOC(tallies(i_tally) % obj % triggers(i_trig))
|
||||
end function
|
||||
|
||||
end module trigger
|
||||
|
|
|
|||
|
|
@ -1,15 +1,18 @@
|
|||
#include "openmc/tallies/trigger.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <utility> // for std::pair
|
||||
|
||||
#include "openmc/capi.h"
|
||||
#include "openmc/constants.h"
|
||||
#include "openmc/message_passing.h"
|
||||
#include "openmc/settings.h"
|
||||
#include "openmc/tallies/tally.h"
|
||||
|
||||
namespace openmc {
|
||||
|
||||
extern "C" void
|
||||
get_tally_uncertainty(int i_tally, int score_index, int filter_index,
|
||||
double* std_dev, double* rel_err)
|
||||
static std::pair<double, double>
|
||||
get_tally_uncertainty(int i_tally, int score_index, int filter_index)
|
||||
{
|
||||
int n;
|
||||
int err = openmc_tally_get_n_realizations(i_tally, &n);
|
||||
|
|
@ -20,13 +23,122 @@ get_tally_uncertainty(int i_tally, int score_index, int filter_index,
|
|||
auto sum_sq = results(filter_index-1, score_index-1, RESULT_SUM_SQ);
|
||||
|
||||
auto mean = sum / n;
|
||||
*std_dev = std::sqrt((sum_sq/n - mean*mean) / (n-1));
|
||||
double std_dev = std::sqrt((sum_sq/n - mean*mean) / (n-1));
|
||||
double rel_err = (mean != 0.) ? std_dev / std::abs(mean) : 0.;
|
||||
|
||||
if (mean > 0) {
|
||||
*rel_err = *std_dev / mean;
|
||||
} else {
|
||||
*rel_err = 0.;
|
||||
return {std_dev, rel_err};
|
||||
}
|
||||
|
||||
// Functions defined in F90.
|
||||
extern "C" int n_tally_triggers(int i_tally);
|
||||
extern "C" Trigger* get_tally_trigger(int i_tally, int i_trig);
|
||||
|
||||
//! Finds 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
|
||||
|
||||
extern "C" void
|
||||
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) {
|
||||
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;
|
||||
|
||||
//TODO: off-by-one
|
||||
for (auto i_trig = 1; i_trig < n_tally_triggers(i_tally)+1; ++i_trig) {
|
||||
auto& trigger {*get_tally_trigger(i_tally, i_trig)};
|
||||
|
||||
trigger.std_dev = 0.;
|
||||
trigger.rel_err = 0.;
|
||||
trigger.variance = 0.;
|
||||
|
||||
const auto& results = tally_results(i_tally);
|
||||
//TODO: off-by-one
|
||||
for (auto filter_index = 1; filter_index < results.shape()[0]+1;
|
||||
++filter_index) {
|
||||
//TODO: off-by-one
|
||||
for (auto score_index = 1; score_index < results.shape()[1]+1;
|
||||
++score_index) {
|
||||
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;
|
||||
}
|
||||
|
||||
double uncertainty;
|
||||
switch (trigger.type) {
|
||||
case VARIANCE:
|
||||
uncertainty = trigger.variance;
|
||||
break;
|
||||
case STANDARD_DEVIATION:
|
||||
uncertainty = trigger.std_dev;
|
||||
break;
|
||||
case RELATIVE_ERROR:
|
||||
uncertainty = trigger.rel_err;
|
||||
}
|
||||
|
||||
double this_ratio = uncertainty / trigger.threshold;
|
||||
if (trigger.type == VARIANCE) {
|
||||
this_ratio = std::sqrt(*ratio);
|
||||
}
|
||||
|
||||
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];
|
||||
err = openmc_tally_get_id(i_tally, tally_id);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//! Computes the uncertainty/threshold ratio for the eigenvalue trigger.
|
||||
|
||||
extern "C" double
|
||||
check_keff_trigger()
|
||||
{
|
||||
if (settings::run_mode != RUN_MODE_EIGENVALUE) return 0.;
|
||||
if (settings::keff_trigger.type == 0) return 0.;
|
||||
|
||||
double k_combined[2];
|
||||
int err = openmc_get_keff(k_combined);
|
||||
|
||||
double uncertainty = 0.;
|
||||
switch (settings::keff_trigger.type) {
|
||||
case VARIANCE:
|
||||
uncertainty = k_combined[1] * k_combined[1];
|
||||
break;
|
||||
case STANDARD_DEVIATION:
|
||||
uncertainty = k_combined[1];
|
||||
break;
|
||||
case RELATIVE_ERROR:
|
||||
uncertainty = k_combined[1] / k_combined[0];
|
||||
}
|
||||
|
||||
double ratio = uncertainty / settings::keff_trigger.threshold;
|
||||
if (settings::keff_trigger.type == VARIANCE)
|
||||
ratio = std::sqrt(ratio);
|
||||
return ratio;
|
||||
}
|
||||
|
||||
} // namespace openmc
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue