diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 360ac0ce8..4dd2d9ae2 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -55,6 +55,14 @@ module tally integer(C_INT), value :: score_index end subroutine + subroutine score_fission_eout(p, i_tally, i_score, score_bin) bind(C) + import Particle, C_INT + type(Particle) :: p + integer(C_INT), value :: i_tally + integer(C_INT), value :: i_score + integer(C_INT), value :: score_bin + end subroutine + subroutine score_analog_tally_ce(p) bind(C) import Particle type(Particle), intent(in) :: p @@ -1229,7 +1237,7 @@ contains ! Add derivative information on score for differential tallies. if (t % deriv() /= C_NONE) then - call apply_derivative_to_score(p, t, i_nuclide, atom_density, & + call apply_derivative_to_score(p, i_tally, i_nuclide, atom_density, & score_bin, score) end if @@ -2037,205 +2045,18 @@ contains end associate end subroutine score_general_mg -!=============================================================================== -! SCORE_FISSION_EOUT handles a special case where we need to store neutron -! production rate with an outgoing energy filter (think of a fission matrix). In -! this case, we may need to score to multiple bins if there were multiple -! neutrons produced with different energies. -!=============================================================================== - - subroutine score_fission_eout(p, i_tally, i_score, score_bin) - - type(Particle), intent(in) :: p - integer, intent(in) :: i_tally - integer, intent(in) :: i_score ! index for score - integer, intent(in) :: score_bin - - integer :: i ! index of outgoing energy filter - integer :: j ! index of delayedgroup filter - integer :: d ! delayed group - integer :: g ! another delayed group - integer :: d_bin ! delayed group bin index - integer :: k ! loop index for bank sites - integer :: l ! loop index for tally filters - integer :: f ! index in filters array - integer :: b ! index of filter bin - integer :: i_match ! matching bin index on energyout filter - integer :: i_bin ! index of matching filter bin - integer :: bin_energyout ! original outgoing energy bin - integer :: i_filter ! index for matching filter bin combination - real(8) :: filter_weight ! combined weight of all filters - real(8) :: score ! actual score - real(8) :: E_out ! energy of fission bank site - integer :: g_out ! energy group of fission bank site - - associate (t => tallies(i_tally) % obj) - - ! save original outgoing energy bin and score index - i = t % filter(t % energyout_filter()) + 1 - i_bin = filter_matches(i) % i_bin() - bin_energyout = filter_matches(i) % bins_data(i_bin) - - ! declare the energyout filter type - select type(eo_filt => filters(i) % obj) - type is (EnergyoutFilter) - - ! Since the creation of fission sites is weighted such that it is - ! expected to create n_particles sites, we need to multiply the - ! score by keff to get the true nu-fission rate. Otherwise, the sum - ! of all nu-fission rates would be ~1.0. - - ! loop over number of particles banked - do k = 1, p % n_bank - - ! get the delayed group - g = fission_bank_delayed_group(n_bank - p % n_bank + k) - - ! determine score based on bank site weight and keff - score = keff * fission_bank_wgt(n_bank - p % n_bank + k) - - ! Add derivative information for differential tallies. Note that the - ! i_nuclide and atom_density arguments do not matter since this is an - ! analog estimator. - if (t % deriv() /= C_NONE) then - call apply_derivative_to_score(p, t, 0, ZERO, SCORE_NU_FISSION, score) - end if - - if (.not. run_CE .and. eo_filt % matches_transport_groups) then - - ! determine outgoing energy group from fission bank - g_out = int(fission_bank_E(n_bank - p % n_bank + k)) - - ! modify the value so that g_out = 1 corresponds to the highest - ! energy bin - g_out = eo_filt % n_bins - g_out + 1 - - ! change outgoing energy bin - call filter_matches(i) % bins_set_data(i_bin, g_out) - - else - - ! determine outgoing energy from fission bank - if (run_CE) then - E_out = fission_bank_E(n_bank - p % n_bank + k) - else - E_out = energy_bin_avg(int(fission_bank_E(n_bank - p % n_bank + k))) - end if - - ! If this outgoing energy falls within the energyout filter's range, - ! set the appropriate filter_matches bin. - i_match = eo_filt % search(E_out) - if (i_match == -1) cycle - call filter_matches(i) % bins_set_data(i_bin, i_match) - - end if - - ! Case for tallying prompt neutrons - if (score_bin == SCORE_NU_FISSION .or. & - (score_bin == SCORE_PROMPT_NU_FISSION .and. g == 0)) then - - ! determine scoring index and weight for this filter combination - i_filter = 1 - do l = 1, t % n_filters() - i_filter = i_filter + (filter_matches(t % filter(l) + 1) & - % bins_data(filter_matches(t % filter(l) + 1) % i_bin()) - 1) * & - t % stride(l) - end do - - ! Add score to tally -!$omp atomic - t % results(RESULT_VALUE, i_score, i_filter) = & - t % results(RESULT_VALUE, i_score, i_filter) + score - - ! Case for tallying delayed emissions - else if (score_bin == SCORE_DELAYED_NU_FISSION .and. g /= 0) then - - ! Get the index of delayed group filter - j = t % delayedgroup_filter() - - ! if the delayed group filter is present, tally to corresponding - ! delayed group bin if it exists - if (j > 0) then - - ! declare the delayed group filter type - select type(dg_filt => filters(t % filter(j) + 1) % obj) - type is (DelayedGroupFilter) - - ! loop over delayed group bins until the corresponding bin is - ! found - do d_bin = 1, dg_filt % n_bins - d = dg_filt % groups(d_bin) - - ! check whether the delayed group of the particle is equal to - ! the delayed group of this bin - if (d == g) then - - ! Reset scoring index and filter weight - i_filter = 1 - filter_weight = ONE - - ! determine scoring index and weight for this filter - ! combination - do l = 1, t % n_filters() - f = t % filter(l) + 1 - b = filter_matches(f) % i_bin() - i_filter = i_filter + (filter_matches(f) & - % bins_data(b) - 1) * t % stride(l) - filter_weight = filter_weight * filter_matches(f) & - % weights_data(b) - end do - - call score_fission_delayed_dg(i_tally, d_bin, & - score * filter_weight, i_score) - end if - end do - end select - - ! if the delayed group filter is not present, add score to tally - else - - ! Reset scoring index and filter weight - i_filter = 1 - filter_weight = ONE - - ! determine scoring index and weight for this filter combination - do l = 1, t % n_filters() - f = t % filter(l) + 1 - b = filter_matches(f) % i_bin() - i_filter = i_filter + (filter_matches(f) % bins_data(b) - 1) & - * t % stride(l) - filter_weight = filter_weight * filter_matches(f) & - % weights_data(b) - end do - - ! Add score to tally -!$omp atomic - t % results(RESULT_VALUE, i_score, i_filter) = & - t % results(RESULT_VALUE, i_score, i_filter) + score * filter_weight - end if - end if - end do - end select - - ! reset outgoing energy bin and score index - call filter_matches(i) % bins_set_data(i_bin, bin_energyout) - - end associate - - end subroutine score_fission_eout - !=============================================================================== ! APPLY_DERIVATIVE_TO_SCORE multiply the given score by its relative derivative !=============================================================================== - subroutine apply_derivative_to_score(p, t, i_nuclide, atom_density, & - score_bin, score) - type(Particle), intent(in) :: p - type(TallyObject), intent(in) :: t - integer, intent(in) :: i_nuclide - real(8), intent(in) :: atom_density ! atom/b-cm - integer, intent(in) :: score_bin - real(8), intent(inout) :: score + subroutine apply_derivative_to_score(p, i_tally, i_nuclide, atom_density, & + score_bin, score) bind(C) + type(Particle), intent(in) :: p + integer(C_INT), value, intent(in) :: i_tally + integer(C_INT), value, intent(in) :: i_nuclide + real(C_DOUBLE), value, intent(in) :: atom_density ! atom/b-cm + integer(C_INT), value, intent(in) :: score_bin + real(C_DOUBLE), intent(inout) :: score type(TallyDerivative), pointer :: deriv integer :: l @@ -2244,6 +2065,8 @@ contains real(8) :: flux_deriv real(8) :: dsig_s, dsig_a, dsig_f, cum_dsig + associate (t => tallies(i_tally) % obj) + if (score == ZERO) return ! If our score was previously c then the new score is @@ -2778,7 +2601,7 @@ contains &analog and collision estimators.") end select end select - !end associate + end associate end subroutine apply_derivative_to_score !=============================================================================== diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 79552db1f..ea948ae40 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -44,6 +44,10 @@ extern "C" void score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, int i_nuclide, double atom_density, double flux); +extern "C" void +apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide, + double atom_density, int score_bin, double* score); + extern "C" int energy_filter_search(const EnergyFilter* filt, double val); @@ -699,6 +703,147 @@ score_fission_delayed_dg(int i_tally, int d_bin, double score, int score_index) dg_match.bins_[i_bin-1] = original_bin; } +//! Helper function for nu-fission tallies with energyout filters. +// +//! In this case, we may need to score to multiple bins if there were multiple +//! neutrons produced with different energies. + +extern "C" void +score_fission_eout(Particle* p, int i_tally, int i_score, int score_bin) +{ + //TODO: off-by-one + const Tally& tally {*model::tallies[i_tally-1]}; + auto results = tally_results(i_tally); + auto i_eout_filt = tally.filters()[tally.energyout_filter_-1]; + auto i_bin = simulation::filter_matches[i_eout_filt].i_bin_; + auto bin_energyout = simulation::filter_matches[i_eout_filt].bins_[i_bin-1]; + + const EnergyoutFilter& eo_filt + {*dynamic_cast(model::tally_filters[i_eout_filt].get())}; + + // Note that the score below is weighted by keff. Since the creation of + // fission sites is weighted such that it is expected to create n_particles + // sites, we need to multiply the score by keff to get the true nu-fission + // rate. Otherwise, the sum of all nu-fission rates would be ~1.0. + + // loop over number of particles banked + for (auto i = 0; i < p->n_bank; ++i) { + auto i_bank = simulation::n_bank - p->n_bank + i; + const auto& bank = simulation::fission_bank[i_bank]; + + // get the delayed group + auto g = bank.delayed_group; + + // determine score based on bank site weight and keff + double score = simulation::keff * bank.wgt; + + // Add derivative information for differential tallies. Note that the + // i_nuclide and atom_density arguments do not matter since this is an + // analog estimator. + if (tally.deriv_ != C_NONE) + apply_derivative_to_score(p, i_tally, 0, 0., SCORE_NU_FISSION, &score); + + if (!settings::run_CE && eo_filt.matches_transport_groups_) { + + // determine outgoing energy group from fission bank + auto g_out = static_cast(bank.E); + + // modify the value so that g_out = 1 corresponds to the highest energy + // bin + g_out = eo_filt.n_bins_ - g_out + 1; + + // change outgoing energy bin + simulation::filter_matches[i_eout_filt].bins_[i_bin-1] = g_out; + + } else { + + double E_out; + if (settings::run_CE) { + E_out = bank.E; + } else { + E_out = data::energy_bin_avg[static_cast(bank.E)]; + } + + //TODO: do this without the extern "C" function + auto i_match = energy_filter_search(&eo_filt, E_out); + if (i_match == -1) continue; + simulation::filter_matches[i_eout_filt].bins_[i_bin-1] = i_match; + + } + + // Case for tallying prompt neutrons + if (score_bin == SCORE_NU_FISSION + || (score_bin == SCORE_PROMPT_NU_FISSION && g == 0)) { + + // Find the filter scoring index for this filter combination + //TODO: should this include a weight? + int filter_index = 1; + for (auto j = 0; j < tally.filters().size(); ++j) { + auto i_filt = tally.filters(j); + auto& match {simulation::filter_matches[i_filt]}; + auto i_bin = match.i_bin_; + filter_index += (match.bins_[i_bin-1] - 1) * tally.strides(j); + } + + // Update tally results + #pragma omp atomic + results(filter_index-1, i_score-1, RESULT_VALUE) += score; + + } else if (score_bin == SCORE_DELAYED_NU_FISSION && g != 0) { + + // Get the index of the delayed group filter + auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_-1]; + + // If the delayed group filter is present, tally to corresponding delayed + // group bin if it exists + if (i_dg_filt >= 0) { + const DelayedGroupFilter& dg_filt {*dynamic_cast( + model::tally_filters[i_dg_filt].get())}; + + // Loop over delayed group bins until the corresponding bin is found + for (auto d_bin = 0; d_bin < dg_filt.n_bins_; ++d_bin) { + if (dg_filt.groups_[d_bin] == g) { + // Find the filter index and weight for this filter combination + int filter_index = 1; + double filter_weight = 1.; + for (auto j = 0; j < tally.filters().size(); ++j) { + auto i_filt = tally.filters(j); + auto& match {simulation::filter_matches[i_filt]}; + auto i_bin = match.i_bin_; + filter_index += (match.bins_[i_bin-1] - 1) * tally.strides(j); + filter_weight *= match.weights_[i_bin-1]; + } + + score_fission_delayed_dg(i_tally, d_bin+1, score*filter_weight, + i_score); + } + } + + // If the delayed group filter is not present, add score to tally + } else { + + // Find the filter index and weight for this filter combination + int filter_index = 1; + double filter_weight = 1.; + for (auto j = 0; j < tally.filters().size(); ++j) { + auto i_filt = tally.filters(j); + auto& match {simulation::filter_matches[i_filt]}; + auto i_bin = match.i_bin_; + filter_index += (match.bins_[i_bin-1] - 1) * tally.strides(j); + filter_weight *= match.weights_[i_bin-1]; + } + + // Update tally results + #pragma omp atomic + results(filter_index-1, i_score-1, RESULT_VALUE) += score*filter_weight; + } + } + } + + // Reset outgoing energy bin and score index + simulation::filter_matches[i_eout_filt].bins_[i_bin-1] = bin_energyout; +} + //! Tally rates for when the user requests a tally on all nuclides. void