Move check_triggers to C++

This commit is contained in:
Sterling Harper 2019-02-01 13:03:00 -05:00
parent b66a4d28b9
commit 5b8806c839
3 changed files with 76 additions and 105 deletions

View file

@ -16,6 +16,7 @@
#include "openmc/timer.h"
#include "openmc/tallies/filter.h"
#include "openmc/tallies/tally.h"
#include "openmc/tallies/trigger.h"
#include <omp.h>
@ -27,7 +28,6 @@ namespace openmc {
// data/functions from Fortran side
extern "C" void accumulate_tallies();
extern "C" void allocate_tally_results();
extern "C" void check_triggers();
extern "C" void init_tally_routines();
extern "C" void load_state_point();
extern "C" void print_batch_keff();

View file

@ -2,106 +2,13 @@ module trigger
use, intrinsic :: ISO_C_BINDING
use constants
use eigenvalue, only: openmc_get_keff
use endf, only: reaction_name
use error, only: warning, write_message
use string, only: to_str
use message_passing, only: master
use settings
use simulation_header
use trigger_header
use tally, only: TallyObject
use tally_header, only: tallies, n_tallies
use tally_header, only: tallies
implicit none
interface
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
contains
!===============================================================================
! CHECK_TRIGGERS checks any user-specified precision triggers' for convergence
! and predicts the number of remainining batches to convergence.
!===============================================================================
subroutine check_triggers() bind(C)
implicit none
real(8) :: keff_ratio ! uncertainty/threshold ratio for keff
real(8) :: tally_ratio ! max tally uncertainty/threshold ratio
integer :: tally_id ! id for tally with max ratio
integer :: score ! tally score with max ratio
integer :: n_pred_batches ! predicted # batches to satisfy all triggers
if (.not. master) return
! Checks if current_batch is one for which the triggers must be checked
if (current_batch < n_batches .or. (.not. trigger_on)) return
if (mod((current_batch - n_batches), n_batch_interval) /= 0 .and. &
current_batch /= n_max_batches) return
! Check the eigenvalue and tally triggers
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 " // &
trim(to_str(current_batch)), 7)
return
end if
! Specify which trigger is unsatisfied
if (keff_ratio >= tally_ratio) then
call write_message("Triggers unsatisfied, max unc./thresh. is " // &
trim(to_str(keff_ratio)) // " for eigenvalue", 7)
else
call write_message("Triggers unsatisfied, max unc./thresh. is " // &
trim(to_str(tally_ratio)) // " for " // trim(reaction_name(score)) &
// " in tally " // trim(to_str(tally_id)), 7)
end if
! Estimate batches till triggers are satisfied
if (pred_batches) then
! Estimate the number of remaining batches to convergence
! The prediction uses the fact that tally variances are proportional
! to 1/N where N is the number of the batches/particles
n_batch_interval = int((current_batch-n_inactive) * &
(max(keff_ratio, tally_ratio) ** 2)) + n_inactive-n_batches + 1
n_pred_batches = n_batch_interval + n_batches
! Write the predicted number of batches for the user
if (n_pred_batches > n_max_batches) then
call warning("The estimated number of batches is " // &
trim(to_str(n_pred_batches)) // &
" -- greater than max batches. ")
else
call write_message("The estimated number of batches is " // &
trim(to_str(n_pred_batches)), 7)
end if
end if
end subroutine check_triggers
function n_tally_triggers(i_tally) bind(C) result(n)
integer(C_INT), value :: i_tally
integer(C_INT) :: n

View file

@ -1,12 +1,15 @@
#include "openmc/tallies/trigger.h"
#include <cmath>
#include <sstream>
#include <utility> // for std::pair
#include "openmc/capi.h"
#include "openmc/constants.h"
#include "openmc/message_passing.h"
#include "openmc/error.h"
#include "openmc/reaction.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/tallies/tally.h"
namespace openmc {
@ -40,10 +43,10 @@ extern "C" Trigger* get_tally_trigger(int i_tally, int i_trig);
//! 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)
void
check_tally_triggers(double& ratio, int& tally_id, int& score)
{
*ratio = 0.;
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]};
@ -94,17 +97,17 @@ check_tally_triggers(double* ratio, int* tally_id, int* score)
double this_ratio = uncertainty / trigger.threshold;
if (trigger.type == VARIANCE) {
this_ratio = std::sqrt(*ratio);
this_ratio = std::sqrt(ratio);
}
if (this_ratio > *ratio) {
*ratio = this_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);
score = scores[trigger.score_index-1];
err = openmc_tally_get_id(i_tally, &tally_id);
}
}
}
@ -114,7 +117,7 @@ check_tally_triggers(double* ratio, int* tally_id, int* score)
//! Computes the uncertainty/threshold ratio for the eigenvalue trigger.
extern "C" double
double
check_keff_trigger()
{
if (settings::run_mode != RUN_MODE_EIGENVALUE) return 0.;
@ -141,4 +144,65 @@ check_keff_trigger()
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("Triggers satisfied for batch "
+ std::to_string(current_batch), 7);
return;
}
// At least one trigger is unsatisfied. Let the user know which one.
simulation::satisfy_triggers = false;
std::stringstream msg;
msg << "Triggers unsatisfied, max unc./thresh. is ";
if (keff_ratio >= tally_ratio) {
msg << keff_ratio << " for eigenvalue";
} else {
msg << tally_ratio << " for " << reaction_name(score) << " in tally "
<< 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;
std::stringstream msg;
msg << "The estimated number of batches is " << n_pred_batches;
if (n_pred_batches > settings::n_max_batches) {
msg << " --- greater than max batches";
warning(msg);
} else {
write_message(msg, 7);
}
}
}
} // namespace openmc