diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 70d20041b0..c657b9722b 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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)) diff --git a/src/tallies/trigger.F90 b/src/tallies/trigger.F90 index 56ed00aba6..785e7b97b0 100644 --- a/src/tallies/trigger.F90 +++ b/src/tallies/trigger.F90 @@ -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). diff --git a/src/tallies/trigger_header.F90 b/src/tallies/trigger_header.F90 index 83d2fa95c0..ccf497d1da 100644 --- a/src/tallies/trigger_header.F90 +++ b/src/tallies/trigger_header.F90 @@ -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