Simplify tally triggers

This commit is contained in:
Sterling Harper 2019-01-31 12:09:09 -05:00
parent abe9bf4a6f
commit a90d8f3e86
3 changed files with 148 additions and 144 deletions

View file

@ -1728,10 +1728,7 @@ contains
call trigger_scores % next_entry(elem, i_elem)
if (i_elem == 0) exit
score_name = trim(elem % key)
! Store the score name and index in the trigger
t % triggers(trig_ind) % score_name = score_name
t % triggers(trig_ind) % score_index = elem % value
! Set the trigger convergence threshold type
@ -1758,7 +1755,6 @@ contains
else
! Store the score name and index
t % triggers(trig_ind) % score_name = trim(score_name)
t % triggers(trig_ind) % score_index = &
trigger_scores % get(trim(score_name))

View file

@ -8,6 +8,7 @@ module trigger
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 mesh_header, only: RegularMesh, meshes
@ -33,44 +34,50 @@ contains
implicit none
! Variables to reflect distance to trigger convergence criteria
real(8) :: max_ratio ! max uncertainty/thresh ratio
integer :: tally_id ! id for tally with max ratio
character(len=52) :: name ! "eigenvalue" or tally score
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
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 trigger and output the result
call check_tally_triggers(max_ratio, tally_id, name)
! By default, assume all triggers are satisfied
satisfy_triggers = .true.
! When trigger threshold is reached, write information
! Check the eigenvalue and tally triggers
call check_keff_trigger(keff_ratio)
call check_tally_triggers(tally_ratio, tally_id, score)
! 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)
! When trigger is not reached write convergence info for user
elseif (name == "eigenvalue") then
call write_message("Triggers unsatisfied, max unc./thresh. is " // &
trim(to_str(max_ratio)) // " for " // trim(name), 7)
else
call write_message("Triggers unsatisfied, max unc./thresh. is " // &
trim(to_str(max_ratio)) // " for " // trim(name) // &
" in tally " // trim(to_str(tally_id)), 7)
return
end if
! If batch_interval is not set, estimate batches till triggers are satisfied
if (pred_batches .and. .not. satisfy_triggers) then
! 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_ratio ** 2)) + n_inactive-n_batches + 1
(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
@ -85,18 +92,56 @@ contains
end if
end subroutine check_triggers
!===============================================================================
! CHECK_TALLY_TRIGGERS checks whether uncertainties are below the threshold,
! and finds the maximum uncertainty/threshold ratio for all 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.
!===============================================================================
subroutine check_tally_triggers(max_ratio, tally_id, name)
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(inout) :: max_ratio ! max uncertainty/thresh ratio
integer, intent(inout) :: tally_id ! id for tally with max ratio
character(len=52), intent(inout) :: name ! "eigenvalue" or tally score
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
@ -114,126 +159,90 @@ contains
! Initialize tally trigger maximum uncertainty ratio to zero
max_ratio = 0
if (master) then
! Compute uncertainties for all tallies, scores with triggers
TALLY_LOOP: do i = 1, n_tallies
associate (t => tallies(i) % obj)
! By default, assume all triggers are satisfied
satisfy_triggers = .true.
! Check eigenvalue trigger
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
if (max_ratio < ratio) then
max_ratio = ratio
name = "eigenvalue"
end if
end if
end if
! Cycle through if only one batch has been simumlate
if (t % n_realizations == 1) then
cycle TALLY_LOOP
end if
! Compute uncertainties for all tallies, scores with triggers
TALLY_LOOP: do i = 1, n_tallies
associate (t => tallies(i) % obj)
TRIGGER_LOOP: do s = 1, t % n_triggers
associate (trigger => t % triggers(s))
! Cycle through if only one batch has been simumlate
if (t % n_realizations == 1) then
cycle TALLY_LOOP
! Initialize trigger uncertainties to zero
trigger % std_dev = ZERO
trigger % rel_err = ZERO
trigger % variance = ZERO
! Mesh current tally triggers require special treatment
if (t % type() == TALLY_MESH_SURFACE) then
call compute_tally_current(t, trigger)
else
! Initialize bins, filter level
do j = 1, t % n_filters()
call filter_matches(t % filter(j)) % bins_clear()
call filter_matches(t % filter(j)) % bins_push_back(0)
end do
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_trigger_uncertainty(std_dev, rel_err, &
score_index, filter_index, t)
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 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
! Mesh current tally triggers require special treatment
if (t % type() == TALLY_MESH_SURFACE) then
call compute_tally_current(t, trigger)
else
! Initialize bins, filter level
do j = 1, t % n_filters()
call filter_matches(t % filter(j)) % bins_clear()
call filter_matches(t % filter(j)) % bins_push_back(0)
end do
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_trigger_uncertainty(std_dev, rel_err, &
score_index, filter_index, t)
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 (t % triggers(s) % type)
case(VARIANCE)
uncertainty = trigger % variance
case(STANDARD_DEVIATION)
uncertainty = trigger % std_dev
case default
uncertainty = trigger % rel_err
end select
if (uncertainty > t % triggers(s) % threshold) then
satisfy_triggers = .false.
if (t % triggers(s) % type == VARIANCE) then
ratio = sqrt(uncertainty / t % triggers(s) % threshold)
else
ratio = uncertainty / t % triggers(s) % threshold
end if
if (max_ratio < ratio) then
max_ratio = ratio
name = t % triggers(s) % score_name
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 if
end associate
end do TRIGGER_LOOP
end associate
end do TALLY_LOOP
end if
end do TRIGGER_LOOP
end associate
end do TALLY_LOOP
end subroutine check_tally_triggers
!===============================================================================
! COMPUTE_TALLY_CURRENT computes the current for a mesh current tally with
! precision trigger(s).

View file

@ -14,7 +14,6 @@ module trigger_header
type, public :: TriggerObject
integer :: type ! "variance", "std_dev" or "rel_err"
real(8) :: threshold ! a convergence threshold
character(len=52) :: score_name ! the name of the score
integer :: score_index ! the index of the score
real(8) :: variance = ZERO ! temp variance container
real(8) :: std_dev = ZERO ! temp std. dev. container