diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 4238f3327..4ac508da0 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -59,7 +59,7 @@ module tally real(C_DOUBLE), intent(in), value :: flux end subroutine - subroutine score_general_mg_c(p, i_tally, start_index, filter_index, & + subroutine score_general_mg(p, i_tally, start_index, filter_index, & i_nuclide, atom_density, flux) bind(C) import Particle, C_INT, C_DOUBLE type(Particle), intent(in) :: p @@ -147,559 +147,6 @@ contains ! analog tallies. !=============================================================================== - subroutine score_general_mg(p, i_tally, start_index, filter_index, i_nuclide, & - atom_density, flux) bind(C) - type(Particle), intent(in) :: p - integer(C_INT), intent(in), value :: i_tally - integer(C_INT), intent(in), value :: start_index - integer(C_INT), intent(in), value :: i_nuclide - integer(C_INT), intent(in), value :: filter_index ! for % results - real(C_DOUBLE), intent(in), value :: flux ! flux estimate - real(C_DOUBLE), intent(in), value :: atom_density ! atom/b-cm - - integer :: i ! loop index for scoring bins - integer :: q ! loop index for scoring bins - integer :: score_bin ! scoring bin, e.g. SCORE_FLUX - integer :: score_index ! scoring bin index - integer :: d ! delayed neutron index - integer :: g ! delayed neutron index - integer :: k ! loop index for bank sites - integer :: d_bin ! delayed group bin index - integer :: dg_filter ! index of delayed group filter - real(8) :: score ! analog tally score - real(8) :: p_uvw(3) ! Particle's current uvw - integer :: p_g ! Particle group to use for getting info - ! to tally with. - - associate (t => tallies(i_tally) % obj) - - ! Set the direction and group to use with get_xs - if (t % estimator() == ESTIMATOR_ANALOG .or. & - t % estimator() == ESTIMATOR_COLLISION) then - - if (survival_biasing) then - - ! Then we either are alive and had a scatter (and so g changed), - ! or are dead and g did not change - if (p % alive) then - p_uvw = p % last_uvw - p_g = p % last_g - else - p_uvw = p % coord(p % n_coord) % uvw - p_g = p % g - end if - else if (p % event == EVENT_SCATTER) then - - ! Then the energy group has been changed by the scattering routine - ! meaning gin is now in p % last_g - p_uvw = p % last_uvw - p_g = p % last_g - else - - ! No scatter, no change in g. - p_uvw = p % coord(p % n_coord) % uvw - p_g = p % g - end if - else - - ! No actual collision so g has not changed. - p_uvw = p % coord(p % n_coord) % uvw - p_g = p % g - end if - - ! To significantly reduce de-referencing, point matxs to the - ! macroscopic Mgxs for the material of interest - call set_macro_angle_index_c(p % material, p_uvw) - - ! Do same for nucxs, point it to the microscopic nuclide data of interest - if (i_nuclide >= 0) then - ! And since we haven't calculated this temperature index yet, do so now - call set_nuclide_temperature_index_c(i_nuclide+1, p % sqrtkT) - call set_nuclide_angle_index_c(i_nuclide+1, p_uvw) - end if - - call score_general_mg_c(p, i_tally, start_index, filter_index, i_nuclide, & - atom_density, flux) - - i = 0 - SCORE_LOOP: do q = 1, t % n_score_bins() - i = i + 1 - - ! determine what type of score bin - score_bin = t % score_bins(i) - - ! determine scoring bin index - score_index = start_index + i - - !######################################################################### - ! Determine appropirate scoring value. - - select case(score_bin) - - - case (SCORE_FLUX) - cycle SCORE_LOOP - - - case (SCORE_TOTAL) - cycle SCORE_LOOP - - - case (SCORE_INVERSE_VELOCITY) - cycle SCORE_LOOP - - - case (SCORE_SCATTER) - cycle SCORE_LOOP - - - case (SCORE_NU_SCATTER) - cycle SCORE_LOOP - - - case (SCORE_ABSORPTION) - cycle SCORE_LOOP - - - case (SCORE_FISSION) - cycle SCORE_LOOP - - - case (SCORE_NU_FISSION) - cycle SCORE_LOOP - - - case (SCORE_PROMPT_NU_FISSION) - - if (t % estimator() == ESTIMATOR_ANALOG) then - if (survival_biasing .or. p % fission) then - if (t % energyout_filter() > 0) then - ! Normally, we only need to make contributions to one scoring - ! bin. However, in the case of fission, since multiple fission - ! neutrons were emitted with different energies, multiple - ! outgoing energy bins may have been scored to. The following - ! logic treats this special case and results to multiple bins - call score_fission_eout(p, i_tally, score_index, score_bin) - cycle SCORE_LOOP - end if - end if - if (survival_biasing) then - ! No fission events occur if survival biasing is on -- need to - ! calculate fraction of absorptions that would have resulted in - ! nu-fission - score = p % absorb_wgt * flux - if (i_nuclide >= 0) then - score = score * atom_density * & - get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_PROMPT_NU_FISSION, p_g) / & - get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) - else - score = score * & - get_macro_xs_c(p % material, MG_GET_XS_PROMPT_NU_FISSION, p_g) / & - get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) - end if - else - ! Skip any non-fission events - if (.not. p % fission) cycle SCORE_LOOP - ! If there is no outgoing energy filter, than we only need to - ! score to one bin. For the score to be 'analog', we need to - ! score the number of particles that were banked in the fission - ! bank. Since this was weighted by 1/keff, we multiply by keff - ! to get the proper score. - score = keff * p % wgt_bank * (ONE - sum(p % n_delayed_bank) & - / real(p % n_bank, 8)) * flux - if (i_nuclide >= 0) then - score = score * atom_density * & - get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_FISSION, p_g) / & - get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g) - end if - end if - - else - if (i_nuclide >= 0) then - score = get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_PROMPT_NU_FISSION, p_g) * & - atom_density * flux - else - score = get_macro_xs_c(p % material, MG_GET_XS_PROMPT_NU_FISSION, p_g) * flux - end if - end if - - - case (SCORE_DELAYED_NU_FISSION) - - ! Set the delayedgroup filter index and the number of delayed group bins - dg_filter = t % delayedgroup_filter() - - if (t % estimator() == ESTIMATOR_ANALOG) then - if (survival_biasing .or. p % fission) then - if (t % energyout_filter() > 0) then - ! Normally, we only need to make contributions to one scoring - ! bin. However, in the case of fission, since multiple fission - ! neutrons were emitted with different energies, multiple - ! outgoing energy bins may have been scored to. The following - ! logic treats this special case and results to multiple bins - call score_fission_eout(p, i_tally, score_index, score_bin) - cycle SCORE_LOOP - end if - end if - if (survival_biasing) then - ! No fission events occur if survival biasing is on -- need to - ! calculate fraction of absorptions that would have resulted in - ! nu-fission - if (get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) > ZERO) then - - if (dg_filter > 0) then - select type(filt => filters(t % filter(dg_filter) + 1) % obj) - type is (DelayedGroupFilter) - - ! Loop over all delayed group bins and tally to them - ! individually - do d_bin = 1, filt % n_bins - - ! Get the delayed group for this bin - d = filt % groups(d_bin) - - score = p % absorb_wgt * flux - if (i_nuclide >= 0) then - score = score * & - get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / & - get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) - else - score = score * & - get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / & - get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) - end if - - call score_fission_delayed_dg(i_tally, d_bin, score, score_index) - end do - cycle SCORE_LOOP - end select - else - score = p % absorb_wgt * flux - if (i_nuclide >= 0) then - score = score * & - get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_DELAYED_NU_FISSION, p_g) / & - get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) - else - score = score * & - get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g) / & - get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) - end if - end if - end if - else - ! Skip any non-fission events - if (.not. p % fission) cycle SCORE_LOOP - ! If there is no outgoing energy filter, than we only need to - ! score to one bin. For the score to be 'analog', we need to - ! score the number of particles that were banked in the fission - ! bank. Since this was weighted by 1/keff, we multiply by keff - ! to get the proper score. - - ! Check if the delayed group filter is present - if (dg_filter > 0) then - select type(filt => filters(t % filter(dg_filter) + 1) % obj) - type is (DelayedGroupFilter) - - ! Loop over all delayed group bins and tally to them - ! individually - do d_bin = 1, filt % n_bins - - ! Get the delayed group for this bin - d = filt % groups(d_bin) - - score = keff * p % wgt_bank / p % n_bank * & - p % n_delayed_bank(d) * flux - - if (i_nuclide >= 0) then - score = score * atom_density * & - get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_FISSION, p_g) / & - get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g) - end if - - call score_fission_delayed_dg(i_tally, d_bin, score, score_index) - end do - cycle SCORE_LOOP - end select - else - score = keff * p % wgt_bank / p % n_bank * sum(p % n_delayed_bank) * flux - if (i_nuclide >= 0) then - score = score * atom_density * & - get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_FISSION, p_g) / & - get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g) - end if - end if - end if - else - - ! Check if the delayed group filter is present - if (dg_filter > 0) then - select type(filt => filters(t % filter(dg_filter) + 1) % obj) - type is (DelayedGroupFilter) - - ! Loop over all delayed group bins and tally to them - ! individually - do d_bin = 1, filt % n_bins - - ! Get the delayed group for this bin - d = filt % groups(d_bin) - - if (i_nuclide >= 0) then - score = atom_density * flux * & - get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) - else - score = flux * & - get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) - end if - - call score_fission_delayed_dg(i_tally, d_bin, score, score_index) - end do - cycle SCORE_LOOP - end select - else - if (i_nuclide >= 0) then - score = atom_density * flux * & - get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_DELAYED_NU_FISSION, p_g) - - else - score = flux * & - get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g) - end if - end if - end if - - case (SCORE_DECAY_RATE) - - ! Set the delayedgroup filter index and the number of delayed group bins - dg_filter = t % delayedgroup_filter() - - if (t % estimator() == ESTIMATOR_ANALOG) then - if (survival_biasing) then - ! No fission events occur if survival biasing is on -- need to - ! calculate fraction of absorptions that would have resulted in - ! nu-fission - if (get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) > ZERO) then - - if (dg_filter > 0) then - select type(filt => filters(t % filter(dg_filter) + 1) % obj) - type is (DelayedGroupFilter) - - ! Loop over all delayed group bins and tally to them - ! individually - do d_bin = 1, filt % n_bins - - ! Get the delayed group for this bin - d = filt % groups(d_bin) - - score = p % absorb_wgt * flux - if (i_nuclide >= 0) then - score = score * & - get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_DECAY_RATE, p_g, DG=d) * & - get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / & - get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) - else - score = score * & - get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * & - get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / & - get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) - end if - - call score_fission_delayed_dg(i_tally, d_bin, score, score_index) - end do - cycle SCORE_LOOP - end select - else - - score = ZERO - - ! If the delayed group filter is not present, compute the score - ! by accumulating the absorbed weight times the decay rate times - ! the fraction of the delayed-nu-fission xs to the absorption xs - ! for all delayed groups. - do d = 1, num_delayed_groups - if (i_nuclide >= 0) then - score = score + p % absorb_wgt * flux * & - get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_DECAY_RATE, p_g, DG=d) * & - get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / & - get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) - else - score = score + p % absorb_wgt * flux * & - get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * & - get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / & - get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) - end if - end do - end if - end if - else - ! Skip any non-fission events - if (.not. p % fission) cycle SCORE_LOOP - ! If there is no outgoing energy filter, than we only need to - ! score to one bin. For the score to be 'analog', we need to - ! score the number of particles that were banked in the fission - ! bank. Since this was weighted by 1/keff, we multiply by keff - ! to get the proper score. - - score = ZERO - - ! 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) - - ! Case for tallying delayed emissions - if (g /= 0) then - - ! determine score based on bank site weight and keff. - if (i_nuclide >= 0) then - score = score + keff * atom_density * & - fission_bank_wgt(n_bank - p % n_bank + k) * & - get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_DECAY_RATE, p_g, DG=g) * & - get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_FISSION, p_g) / & - get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g) * flux - else - score = score + keff * & - fission_bank_wgt(n_bank - p % n_bank + k) * & - get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=g) * flux - end if - - ! if the delayed group filter is present, tally to corresponding - ! delayed group bin if it exists - if (dg_filter > 0) then - - ! declare the delayed group filter type - select type(filt => filters(t % filter(dg_filter) + 1) % obj) - type is (DelayedGroupFilter) - - ! loop over delayed group bins until the corresponding bin - ! is found - do d_bin = 1, filt % n_bins - d = 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 - call score_fission_delayed_dg(i_tally, d_bin, score, & - score_index) - end if - end do - end select - - ! Reset the score to zero - score = ZERO - end if - end if - end do - - ! If the delayed group filter is present, cycle because the - ! score_fission_delayed_dg(...) has already tallied the score - if (dg_filter > 0) then - cycle SCORE_LOOP - end if - end if - else - - ! Check if the delayed group filter is present - if (dg_filter > 0) then - select type(filt => filters(t % filter(dg_filter) + 1) % obj) - type is (DelayedGroupFilter) - - ! Loop over all delayed group bins and tally to them - ! individually - do d_bin = 1, filt % n_bins - - ! Get the delayed group for this bin - d = filt % groups(d_bin) - - if (i_nuclide >= 0) then - score = atom_density * flux * & - get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_DECAY_RATE, p_g, DG=d) * & - get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) - else - score = flux * & - get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * & - get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) - end if - - call score_fission_delayed_dg(i_tally, d_bin, score, score_index) - end do - cycle SCORE_LOOP - end select - else - score = ZERO - - ! If the delayed group filter is not present, compute the score - ! by accumulating the absorbed weight times the decay rate times - ! the fraction of the delayed-nu-fission xs to the absorption xs - ! for all delayed groups. - do d = 1, num_delayed_groups - if (i_nuclide >= 0) then - score = score + atom_density * flux * & - get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_DECAY_RATE, p_g, DG=d) * & - get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) - else - score = score + flux * & - get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * & - get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) - end if - end do - end if - end if - - - case (SCORE_KAPPA_FISSION) - - if (t % estimator() == ESTIMATOR_ANALOG) then - if (survival_biasing) then - ! No fission events occur if survival biasing is on -- need to - ! calculate fraction of absorptions that would have resulted in - ! fission - score = p % absorb_wgt * flux - else - ! Skip any non-absorption events - if (p % event == EVENT_SCATTER) cycle SCORE_LOOP - ! All fission events will contribute, so again we can use - ! particle's weight entering the collision as the estimate for the - ! fission reaction rate - score = p % last_wgt * flux - end if - if (i_nuclide >= 0) then - score = score * atom_density * & - get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_KAPPA_FISSION, p_g) / & - get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) - else - score = score * & - get_macro_xs_c(p % material, MG_GET_XS_KAPPA_FISSION, p_g) / & - get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) - end if - else - if (i_nuclide >= 0) then - score = get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_KAPPA_FISSION, p_g) * & - atom_density * flux - else - score = flux * & - get_macro_xs_c(p % material, MG_GET_XS_KAPPA_FISSION, p_g) - - end if - end if - - - case (SCORE_EVENTS) - ! Simply count number of scoring events - score = ONE - - end select - - !######################################################################### - ! Expand score if necessary and add to tally results. -!$omp atomic - t % results(RESULT_VALUE, score_index, filter_index) = & - t % results(RESULT_VALUE, score_index, filter_index) + score - - end do SCORE_LOOP - end associate - end subroutine score_general_mg - !=============================================================================== ! APPLY_DERIVATIVE_TO_SCORE multiply the given score by its relative derivative !=============================================================================== diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 46d5ff60a..d62224d63 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -1761,13 +1761,16 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, } extern "C" void -score_general_mg_c(Particle* p, int i_tally, int start_index, int filter_index, +score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, int i_nuclide, double atom_density, double flux) { //TODO: off-by-one const Tally& tally {*model::tallies[i_tally-1]}; auto results = tally_results(i_tally); + //TODO: off-by-one throughout on score_index in calls to score_fission_eout + //TODO: off-by-one throughout on p->material + // Set the direction and group to use with get_xs double* p_uvw; int p_g; @@ -1804,9 +1807,16 @@ score_general_mg_c(Particle* p, int i_tally, int start_index, int filter_index, p_g = p->g; } - //TODO: set_macro_angle_index - //TODO: set_nuclide_temperature_index - //TODO: est_nuclide_angle_index + // To significantly reduce de-referencing, point matxs to the macroscopic + // Mgxs for the material of interest + data::macro_xs[p->material - 1].set_angle_index(p_uvw); + + // Do same for nucxs, point it to the microscopic nuclide data of interest + if (i_nuclide >= 0) { + // And since we haven't calculated this temperature index yet, do so now + data::nuclides_MG[i_nuclide].set_temperature_index(p->sqrtkT); + data::nuclides_MG[i_nuclide].set_angle_index(p_uvw); + } for (auto i = 0; i < tally.scores_.size(); ++i) { auto score_bin = tally.scores_[i]; @@ -2067,6 +2077,384 @@ score_general_mg_c(Particle* p, int i_tally, int start_index, int filter_index, break; + case SCORE_PROMPT_NU_FISSION: + if (tally.estimator_ == ESTIMATOR_ANALOG) { + if (settings::survival_biasing || p->fission) { + if (tally.energyout_filter_ > 0) { + // Fission has multiple outgoing neutrons so this helper function + // is used to handle scoring the multiple filter bins. + score_fission_eout(p, i_tally, score_index+1, score_bin); + continue; + } + } + if (settings::survival_biasing) { + // No fission events occur if survival biasing is on -- need to + // calculate fraction of absorptions that would have resulted in + // prompt-nu-fission + score = p->absorb_wgt * flux; + if (i_nuclide >= 0) { + score *= atom_density + * get_nuclide_xs(i_nuclide+1, MG_GET_XS_PROMPT_NU_FISSION, p_g) + / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + } else { + score *= + get_macro_xs(p->material, MG_GET_XS_PROMPT_NU_FISSION, p_g) + / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + } + } else { + // Skip any non-fission events + if (!p->fission) continue; + // If there is no outgoing energy filter, than we only need to score + // to one bin. For the score to be 'analog', we need to score the + // number of particles that were banked in the fission bank. Since + // this was weighted by 1/keff, we multiply by keff to get the proper + // score. + auto n_delayed = std::accumulate(p->n_delayed_bank, + p->n_delayed_bank+MAX_DELAYED_GROUPS, 0); + auto prompt_frac = 1. - n_delayed / static_cast(p->n_bank); + score = simulation::keff * p->wgt_bank * prompt_frac * flux; + if (i_nuclide >= 0) { + score *= atom_density + * get_nuclide_xs(i_nuclide+1, MG_GET_XS_FISSION, p_g) + / get_macro_xs(p->material, MG_GET_XS_FISSION, p_g); + } + } + } else { + if (i_nuclide >= 0) { + score = get_nuclide_xs(i_nuclide+1, MG_GET_XS_PROMPT_NU_FISSION, p_g) + * atom_density * flux; + } else { + score = get_macro_xs(p->material, MG_GET_XS_PROMPT_NU_FISSION, p_g) + * flux; + } + } + break; + + + case SCORE_DELAYED_NU_FISSION: + if (tally.estimator_ == ESTIMATOR_ANALOG) { + if (settings::survival_biasing || p->fission) { + if (tally.energyout_filter_ > 0) { + // Fission has multiple outgoing neutrons so this helper function + // is used to handle scoring the multiple filter bins. + score_fission_eout(p, i_tally, score_index+1, score_bin); + continue; + } + } + if (settings::survival_biasing) { + // No fission events occur if survival biasing is on -- need to + // calculate fraction of absorptions that would have resulted in + // delayed-nu-fission + if (get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g) > 0) { + if (tally.delayedgroup_filter_ > 0) { + auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_-1]; + const DelayedGroupFilter& filt + {*dynamic_cast( + model::tally_filters[i_dg_filt].get())}; + // Tally each delayed group bin individually + for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { + auto d = filt.groups_[d_bin]; + score = p->absorb_wgt * flux; + if (i_nuclide >= 0) { + score *= + get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_DELAYED_NU_FISSION, + p_g, nullptr, nullptr, &d) + / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + } else { + score *= + get_macro_xs_c(p->material, MG_GET_XS_DELAYED_NU_FISSION, + p_g, nullptr, nullptr, &d) + / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + } + score_fission_delayed_dg(i_tally, d_bin+1, score, + score_index+1); + } + continue; + } else { + // If the delayed group filter is not present, compute the score + // by multiplying the absorbed weight by the fraction of the + // delayed-nu-fission xs to the absorption xs + score = p->absorb_wgt * flux; + if (i_nuclide >= 0) { + score *= + get_nuclide_xs(i_nuclide+1, MG_GET_XS_DELAYED_NU_FISSION, p_g) + / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + } else { + score *= + get_macro_xs(p->material, MG_GET_XS_DELAYED_NU_FISSION, p_g) + / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + } + } + } + } else { + // Skip any non-fission events + if (!p->fission) continue; + // If there is no outgoing energy filter, than we only need to score + // to one bin. For the score to be 'analog', we need to score the + // number of particles that were banked in the fission bank. Since + // this was weighted by 1/keff, we multiply by keff to get the proper + // score. Loop over the neutrons produced from fission and check which + // ones are delayed. If a delayed neutron is encountered, add its + // contribution to the fission bank to the score. + if (tally.delayedgroup_filter_ > 0) { + auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_-1]; + const DelayedGroupFilter& filt + {*dynamic_cast( + model::tally_filters[i_dg_filt].get())}; + // Tally each delayed group bin individually + for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { + auto d = filt.groups_[d_bin]; + score = simulation::keff * p->wgt_bank / p->n_bank + * p->n_delayed_bank[d-1] * flux; + if (i_nuclide >= 0) { + score *= atom_density + * get_nuclide_xs(i_nuclide+1, MG_GET_XS_FISSION, p_g) + / get_macro_xs(p->material, MG_GET_XS_FISSION, p_g); + } + score_fission_delayed_dg(i_tally, d_bin+1, score, score_index+1); + } + continue; + } else { + // Add the contribution from all delayed groups + auto n_delayed = std::accumulate(p->n_delayed_bank, + p->n_delayed_bank+MAX_DELAYED_GROUPS, 0); + score = simulation::keff * p->wgt_bank / p->n_bank * n_delayed + * flux; + if (i_nuclide >= 0) { + score *= atom_density + * get_nuclide_xs(i_nuclide+1, MG_GET_XS_FISSION, p_g) + / get_macro_xs(p->material, MG_GET_XS_FISSION, p_g); + } + } + } + } else { + if (tally.delayedgroup_filter_ > 0) { + auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_-1]; + const DelayedGroupFilter& filt + {*dynamic_cast( + model::tally_filters[i_dg_filt].get())}; + // Tally each delayed group bin individually + for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { + auto d = filt.groups_[d_bin]; + if (i_nuclide >= 0) { + score = flux * atom_density + * get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_DELAYED_NU_FISSION, + p_g, nullptr, nullptr, &d); + } else { + score = flux + * get_macro_xs_c(p->material, MG_GET_XS_DELAYED_NU_FISSION, + p_g, nullptr, nullptr, &d); + } + score_fission_delayed_dg(i_tally, d_bin+1, score, score_index+1); + } + continue; + } else { + if (i_nuclide >= 0) { + score = flux * atom_density + * get_nuclide_xs(i_nuclide+1, MG_GET_XS_DELAYED_NU_FISSION, p_g); + } else { + score = flux + * get_macro_xs(p->material, MG_GET_XS_DELAYED_NU_FISSION, p_g); + } + } + } + break; + + + case SCORE_DECAY_RATE: + if (tally.estimator_ == ESTIMATOR_ANALOG) { + if (settings::survival_biasing) { + // No fission events occur if survival biasing is on -- need to + // calculate fraction of absorptions that would have resulted in + // delayed-nu-fission + if (get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g) > 0) { + if (tally.delayedgroup_filter_ > 0) { + auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_-1]; + const DelayedGroupFilter& filt + {*dynamic_cast( + model::tally_filters[i_dg_filt].get())}; + // Tally each delayed group bin individually + for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { + auto d = filt.groups_[d_bin]; + score = p->absorb_wgt * flux; + if (i_nuclide >= 0) { + score *= + get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_DECAY_RATE, + p_g, nullptr, nullptr, &d) + *get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_DELAYED_NU_FISSION, + p_g, nullptr, nullptr, &d) + / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + } else { + score *= + get_macro_xs_c(p->material, MG_GET_XS_DECAY_RATE, + p_g, nullptr, nullptr, &d) + * get_macro_xs_c(p->material, MG_GET_XS_DELAYED_NU_FISSION, + p_g, nullptr, nullptr, &d) + / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + } + score_fission_delayed_dg(i_tally, d_bin+1, score, + score_index+1); + } + continue; + } else { + // If the delayed group filter is not present, compute the score + // by multiplying the absorbed weight by the fraction of the + // delayed-nu-fission xs to the absorption xs for all delayed + // groups + score = 0.; + for (auto d = 0; d < data::num_delayed_groups; ++d) { + if (i_nuclide >= 0) { + score += p->absorb_wgt * flux + * get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_DECAY_RATE, + p_g, nullptr, nullptr, &d) + *get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_DELAYED_NU_FISSION, + p_g, nullptr, nullptr, &d) + / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + } else { + score += p->absorb_wgt * flux + * get_macro_xs_c(p->material, MG_GET_XS_DECAY_RATE, + p_g, nullptr, nullptr, &d) + * get_macro_xs_c(p->material, MG_GET_XS_DELAYED_NU_FISSION, + p_g, nullptr, nullptr, &d) + / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + } + } + } + } + } else { + // Skip any non-fission events + if (!p->fission) continue; + // If there is no outgoing energy filter, than we only need to score + // to one bin. For the score to be 'analog', we need to score the + // number of particles that were banked in the fission bank. Since + // this was weighted by 1/keff, we multiply by keff to get the proper + // score. Loop over the neutrons produced from fission and check which + // ones are delayed. If a delayed neutron is encountered, add its + // contribution to the fission bank to the score. + score = 0.; + 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]; + auto g = bank.delayed_group; + if (g != 0) { + if (i_nuclide >= 0) { + score += simulation::keff * atom_density * bank.wgt * flux + * get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_DECAY_RATE, p_g, + nullptr, nullptr, &g) + * get_nuclide_xs(i_nuclide+1, MG_GET_XS_FISSION, p_g) + / get_macro_xs(p->material, MG_GET_XS_FISSION, p_g); + } else { + score += simulation::keff * bank.wgt * flux + * get_macro_xs_c(p->material, MG_GET_XS_DECAY_RATE, p_g, + nullptr, nullptr, &g); + } + if (tally.delayedgroup_filter_ > 0) { + auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_-1]; + const DelayedGroupFilter& filt + {*dynamic_cast( + model::tally_filters[i_dg_filt].get())}; + // Find the corresponding filter bin and then score + for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { + auto d = filt.groups_[d_bin]; + if (d == g) + score_fission_delayed_dg(i_tally, d_bin+1, score, + score_index+1); + } + score = 0.; + } + } + } + if (tally.delayedgroup_filter_ > 0) continue; + } + } else { + if (tally.delayedgroup_filter_ > 0) { + auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_-1]; + const DelayedGroupFilter& filt + {*dynamic_cast( + model::tally_filters[i_dg_filt].get())}; + // Tally each delayed group bin individually + for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { + auto d = filt.groups_[d_bin]; + if (i_nuclide >= 0) { + score += atom_density * flux + * get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_DECAY_RATE, + p_g, nullptr, nullptr, &d) + * get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_DELAYED_NU_FISSION, + p_g, nullptr, nullptr, &d); + } else { + score += flux + * get_macro_xs_c(p->material, MG_GET_XS_DECAY_RATE, + p_g, nullptr, nullptr, &d) + * get_macro_xs_c(p->material, MG_GET_XS_DELAYED_NU_FISSION, + p_g, nullptr, nullptr, &d); + } + score_fission_delayed_dg(i_tally, d_bin+1, score, score_index+1); + } + continue; + } else { + score = 0.; + for (auto d = 0; d < data::num_delayed_groups; ++d) { + if (i_nuclide >= 0) { + score += atom_density * flux + * get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_DECAY_RATE, + p_g, nullptr, nullptr, &d) + * get_nuclide_xs_c(i_nuclide+1, MG_GET_XS_DELAYED_NU_FISSION, + p_g, nullptr, nullptr, &d); + } else { + score += flux + * get_macro_xs_c(p->material, MG_GET_XS_DECAY_RATE, + p_g, nullptr, nullptr, &d) + * get_macro_xs_c(p->material, MG_GET_XS_DELAYED_NU_FISSION, + p_g, nullptr, nullptr, &d); + } + } + } + } + break; + + + case SCORE_KAPPA_FISSION: + if (tally.estimator_ == ESTIMATOR_ANALOG) { + if (settings::survival_biasing) { + // No fission events occur if survival biasing is on -- need to + // calculate fraction of absorptions that would have resulted in + // fission scaled by the Q-value + score = p->absorb_wgt * flux; + } else { + // Skip any non-absorption events + if (p->event == EVENT_SCATTER) continue; + // All fission events will contribute, so again we can use particle's + // weight entering the collision as the estimate for the fission + // reaction rate + score = p->last_wgt * flux; + } + if (i_nuclide >= 0) { + score *= atom_density + * get_nuclide_xs(i_nuclide+1, MG_GET_XS_KAPPA_FISSION, p_g) + / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + } else { + score *= + get_macro_xs(p->material, MG_GET_XS_KAPPA_FISSION, p_g) + / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + } + } else { + if (i_nuclide >= 0) { + score = get_nuclide_xs(i_nuclide+1, MG_GET_XS_KAPPA_FISSION, p_g) + * atom_density * flux; + } else { + score = get_macro_xs(p->material, MG_GET_XS_KAPPA_FISSION, p_g) + * flux; + } + } + break; + + + case SCORE_EVENTS: + // Simply count the number of scoring events + score = 1.; + break; + + default: continue; }