From 1c1fa6d753f187ded9209f208e9bb7342e721adc Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 6 Feb 2019 21:03:57 -0500 Subject: [PATCH] Finish moving score_general_ce to C++ --- src/nuclide.cpp | 2 +- src/tallies/tally.F90 | 709 +----------------------------------------- src/tallies/tally.cpp | 455 ++++++++++++++++++++++++++- 3 files changed, 462 insertions(+), 704 deletions(-) diff --git a/src/nuclide.cpp b/src/nuclide.cpp index b64a1681ab..b507b9130e 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -274,7 +274,7 @@ void Nuclide::create_derived() xs_.emplace_back(shape, 0.0); } - reaction_index_.fill(-1); + reaction_index_.fill(C_NONE); for (int i = 0; i < reactions_.size(); ++i) { const auto& rx {reactions_[i]}; diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index cc3b169e42..57560910aa 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -47,16 +47,16 @@ module tally end interface interface - subroutine score_general_ce_c(p, i_tally, start_index, filter_index, & + subroutine score_general_ce(p, i_tally, start_index, filter_index, & i_nuclide, atom_density, flux) bind(C) import Particle, C_INT, C_DOUBLE - type(Particle) :: p - integer(C_INT), value :: i_tally - integer(C_INT), value :: start_index - integer(C_INT), value :: filter_index - integer(C_INT), value :: i_nuclide - real(C_DOUBLE), value :: atom_density - real(C_DOUBLE), value :: flux + 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 :: filter_index + integer(C_INT), intent(in), value :: i_nuclide + real(C_DOUBLE), intent(in), value :: atom_density + real(C_DOUBLE), intent(in), value :: flux end subroutine subroutine score_fission_delayed_dg(i_tally, d_bin, score, score_index) bind(C) @@ -135,699 +135,6 @@ contains ! analog tallies. !=============================================================================== - subroutine score_general_ce(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 :: l ! loop index for nuclides in material - integer :: m ! loop index for reactions - integer :: i_temp ! temperature index - integer :: i_nuc ! index in nuclides array (from material) - integer :: i_energy ! index in nuclide energy grid - 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 - integer :: threshold ! threshold energy index - real(8) :: yield ! delayed neutron yield - real(8) :: atom_density_ ! atom/b-cm - real(8) :: f ! interpolation factor - real(8) :: score ! analog tally score - real(8) :: E ! particle energy - real(8) :: xs ! cross section - - call score_general_ce_c(p, i_tally, start_index, filter_index, & - i_nuclide, atom_density, flux) - - associate (t => tallies(i_tally) % obj) - - ! Pre-collision energy of particle - E = p % last_E - - SCORE_LOOP: do i = 1, t % n_score_bins() - - ! 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) - cycle SCORE_LOOP - - case (SCORE_DELAYED_NU_FISSION) - cycle SCORE_LOOP - - case (SCORE_DECAY_RATE) - - if (material_xs % absorption == ZERO) cycle SCORE_LOOP - - ! Set the delayedgroup filter index - 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 - ! delayed-nu-fission - if (micro_xs(p % event_nuclide) % absorption > ZERO .and. & - nuclides(p % event_nuclide) % fissionable) then - - ! 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) - - ! Compute the yield for this delayed group - yield = nuclides(p % event_nuclide) & - % nu(E, EMISSION_DELAYED, d) - - associate (rxn => nuclides(p % event_nuclide) % & - reactions(nuclides(p % event_nuclide) % & - index_fission(1))) - - ! Compute the score - score = p % absorb_wgt * yield * & - micro_xs(p % event_nuclide) % fission & - / micro_xs(p % event_nuclide) % absorption & - * rxn % product_decay_rate(1 + d) * flux - end associate - - ! Tally to bin - call score_fission_delayed_dg(i_tally, d_bin, score, score_index) - end do - cycle SCORE_LOOP - end select - else - - ! 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. - score = ZERO - - associate (rxn => nuclides(p % event_nuclide) % & - reactions(nuclides(p % event_nuclide) % index_fission(1))) - - ! We need to be careful not to overshoot the number of delayed - ! groups since this could cause the range of the - ! rxn % products array to be exceeded. Hence, we use the size - ! of this array and not the MAX_DELAYED_GROUPS constant for - ! this loop. - do d = 1, rxn % products_size() - 2 - score = score + rxn % product_decay_rate(1 + d) * & - p % absorb_wgt & - * micro_xs(p % event_nuclide) % fission & - * nuclides(p % event_nuclide) % & - nu(E, EMISSION_DELAYED, d) & - / micro_xs(p % event_nuclide) % absorption * flux - end do - end associate - 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. 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 = 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 - - ! Accumulate the decay rate times delayed nu fission score - associate (rxn => nuclides(p % event_nuclide) % & - reactions(nuclides(p % event_nuclide) % index_fission(1))) - - ! determine score based on bank site weight and keff. - score = score + keff * fission_bank_wgt(n_bank - p % n_bank + k) & - * rxn % product_decay_rate(1 + g) * flux - end associate - - ! 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 - end if - else - - ! Check if tally is on a single nuclide - if (i_nuclide >= 0) then - - ! 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) - - ! Compute the yield for this delayed group - yield = nuclides(i_nuclide+1) % nu(E, EMISSION_DELAYED, d) - - associate (rxn => nuclides(i_nuclide+1) % & - reactions(nuclides(i_nuclide+1) % index_fission(1))) - - ! Compute the score and tally to bin - score = micro_xs(i_nuclide+1) % fission * yield * flux * & - atom_density * rxn % product_decay_rate(1 + d) - end associate - - ! Tally to bin - call score_fission_delayed_dg(i_tally, d_bin, score, score_index) - end do - cycle SCORE_LOOP - end select - else - - ! 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. - score = ZERO - - associate (rxn => nuclides(p % event_nuclide) % & - reactions(nuclides(p % event_nuclide) % index_fission(1))) - - ! We need to be careful not to overshoot the number of delayed - ! groups since this could cause the range of the rxn % products - ! array to be exceeded. Hence, we use the size of this array - ! and not the MAX_DELAYED_GROUPS constant for this loop. - do d = 1, rxn % products_size() - 2 - score = score + micro_xs(i_nuclide+1) % fission * flux * & - nuclides(i_nuclide+1) % nu(E, EMISSION_DELAYED) * & - atom_density * rxn % product_decay_rate(1 + d) - end do - end associate - end if - - ! Tally is on total nuclides - 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 nuclides in the current material - if (p % material /= MATERIAL_VOID) then - do l = 1, material_nuclide_size(p % material) - - ! Get atom density - atom_density_ = material_atom_density(p % material, l) - - ! Get index in nuclides array - i_nuc = material_nuclide(p % material, l) - - if (nuclides(i_nuc) % fissionable) then - - ! 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) - - ! Get the yield for the desired nuclide and delayed - ! group - yield = nuclides(i_nuc) % nu(E, EMISSION_DELAYED, d) - - associate (rxn => nuclides(i_nuc) % & - reactions(nuclides(i_nuc) % index_fission(1))) - - ! Compute the score - score = micro_xs(i_nuc) % fission * yield * flux * & - atom_density_ & - * rxn % product_decay_rate(1 + d) - end associate - - ! Tally to bin - call score_fission_delayed_dg(i_tally, d_bin, score, & - score_index) - end do - end if - end do - end if - cycle SCORE_LOOP - end select - else - - score = ZERO - - ! Loop over all nuclides in the current material - if (p % material /= MATERIAL_VOID) then - do l = 1, material_nuclide_size(p % material) - - ! Get atom density - atom_density_ = material_atom_density(p % material, l) - - ! Get index in nuclides array - i_nuc = material_nuclide(p % material, l) - - if (nuclides(i_nuc) % fissionable) then - - associate (rxn => nuclides(i_nuc) % & - reactions(nuclides(i_nuc) % index_fission(1))) - - ! We need to be careful not to overshoot the number of - ! delayed groups since this could cause the range of the - ! rxn % products array to be exceeded. Hence, we use the - ! size of this array and not the MAX_DELAYED_GROUPS - ! constant for this loop. - do d = 1, rxn % products_size() - 2 - - ! Accumulate the contribution from each nuclide - score = score + micro_xs(i_nuc) % fission & - * nuclides(i_nuc) % nu(E, EMISSION_DELAYED) & - * atom_density_ * flux & - * rxn % product_decay_rate(1 + d) - end do - end associate - end if - end do - end if - end if - end if - end if - - case (SCORE_KAPPA_FISSION) - - if (material_xs % absorption == ZERO) cycle SCORE_LOOP - - ! Determine kappa-fission cross section on the fly. The ENDF standard - ! (ENDF-102) states that MT 18 stores the fission energy as the Q_value - ! (fission(1)) - - score = ZERO - - 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 scaled by kappa-fission - associate (nuc => nuclides(p % event_nuclide)) - if (micro_xs(p % event_nuclide) % absorption > ZERO .and. & - nuc % fissionable) then - score = p % absorb_wgt * & - nuc % reactions(nuc % index_fission(1)) % Q_value * & - micro_xs(p % event_nuclide) % fission / & - micro_xs(p % event_nuclide) % absorption * flux - end if - end associate - 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 energy production rate - associate (nuc => nuclides(p % event_nuclide)) - if (nuc % fissionable) then - score = p % last_wgt * & - nuc % reactions(nuc % index_fission(1)) % Q_value * & - micro_xs(p % event_nuclide) % fission / & - micro_xs(p % event_nuclide) % absorption * flux - end if - end associate - end if - - else - if (i_nuclide >= 0) then - associate (nuc => nuclides(i_nuclide+1)) - if (nuc % fissionable) then - score = nuc % reactions(nuc % index_fission(1)) % Q_value * & - micro_xs(i_nuclide+1) % fission * atom_density * flux - end if - end associate - else - if (p % material == MATERIAL_VOID) then - score = ZERO - else - do l = 1, material_nuclide_size(p % material) - ! Determine atom density and index of nuclide - atom_density_ = material_atom_density(p % material, l) - i_nuc = material_nuclide(p % material, l) - - ! If nuclide is fissionable, accumulate kappa fission - associate(nuc => nuclides(i_nuc)) - if (nuc % fissionable) then - score = score + & - nuc % reactions(nuc % index_fission(1)) % Q_value * & - micro_xs(i_nuc) % fission * atom_density_ * flux - end if - end associate - end do - end if - end if - end if - - case (SCORE_EVENTS) - ! Simply count number of scoring events - score = ONE - - case (ELASTIC) - if (t % estimator() == ESTIMATOR_ANALOG) then - ! Check if event MT matches - if (p % event_MT /= ELASTIC) cycle SCORE_LOOP - score = p % last_wgt * flux - - else - if (i_nuclide >= 0) then - if (micro_xs(i_nuclide+1) % elastic == CACHE_INVALID) then - call nuclides(i_nuclide+1) % calculate_elastic_xs() - end if - score = micro_xs(i_nuclide+1) % elastic * atom_density * flux - else - score = ZERO - if (p % material /= MATERIAL_VOID) then - do l = 1, material_nuclide_size(p % material) - ! Get atom density - atom_density_ = material_atom_density(p % material, l) - - ! Get index in nuclides array - i_nuc = material_nuclide(p % material, l) - if (micro_xs(i_nuc) % elastic == CACHE_INVALID) then - call nuclides(i_nuc) % calculate_elastic_xs() - end if - - score = score + micro_xs(i_nuc) % elastic * atom_density_ * flux - end do - end if - end if - end if - - case (SCORE_FISS_Q_PROMPT, SCORE_FISS_Q_RECOV) - - if (material_xs % absorption == ZERO) cycle SCORE_LOOP - - score = ZERO - - 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 scaled by Q-value - associate (nuc => nuclides(p % event_nuclide)) - if (micro_xs(p % event_nuclide) % absorption > ZERO) then - if (score_bin == SCORE_FISS_Q_PROMPT) then - xs = nuclide_fission_q_prompt(nuc % ptr, p % last_E) - else if (score_bin == SCORE_FISS_Q_RECOV) then - xs = nuclide_fission_q_recov(nuc % ptr, p % last_E) - end if - - score = p % absorb_wgt * xs * flux & - * micro_xs(p % event_nuclide) % fission & - / micro_xs(p % event_nuclide) % absorption - end if - end associate - 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 energy production rate - associate (nuc => nuclides(p % event_nuclide)) - if (micro_xs(p % event_nuclide) % absorption > ZERO) then - if (score_bin == SCORE_FISS_Q_PROMPT) then - xs = nuclide_fission_q_prompt(nuc % ptr, p % last_E) - else if (score_bin == SCORE_FISS_Q_RECOV) then - xs = nuclide_fission_q_recov(nuc % ptr, p % last_E) - end if - - score = p % last_wgt * xs * flux & - * micro_xs(p % event_nuclide) % fission & - / micro_xs(p % event_nuclide) % absorption - end if - end associate - end if - - else - if (i_nuclide >= 0) then - associate (nuc => nuclides(i_nuclide+1)) - if (score_bin == SCORE_FISS_Q_PROMPT) then - xs = nuclide_fission_q_prompt(nuc % ptr, E) - else if (score_bin == SCORE_FISS_Q_RECOV) then - xs = nuclide_fission_q_recov(nuc % ptr, E) - end if - - score = micro_xs(i_nuclide+1) % fission * atom_density * flux * xs - end associate - else - if (p % material /= MATERIAL_VOID) then - do l = 1, material_nuclide_size(p % material) - atom_density_ = material_atom_density(p % material, l) - i_nuc = material_nuclide(p % material, l) - - associate (nuc => nuclides(i_nuc)) - if (score_bin == SCORE_FISS_Q_PROMPT) then - xs = nuclide_fission_q_prompt(nuc % ptr, E) - else if (score_bin == SCORE_FISS_Q_RECOV) then - xs = nuclide_fission_q_recov(nuc % ptr, E) - end if - - score = score + micro_xs(i_nuc) % fission * atom_density_ & - * flux * xs - end associate - end do - end if - end if - end if - - case (N_2N, N_3N, N_4N, N_GAMMA, N_P, N_A) - if (t % estimator() == ESTIMATOR_ANALOG) then - ! Check if event MT matches - if (p % event_MT /= score_bin) cycle SCORE_LOOP - score = p % last_wgt * flux - - else - ! Determine index in NuclideMicroXS % reaction array - select case (score_bin) - case (N_GAMMA) - m = 1 - case (N_P) - m = 2 - case (N_A) - m = 3 - case (N_2N) - m = 4 - case (N_3N) - m = 5 - case (N_4N) - m = 6 - end select - - if (i_nuclide >= 0) then - score = micro_xs(i_nuclide+1) % reaction(m) * atom_density * flux - else - score = ZERO - if (p % material /= MATERIAL_VOID) then - do l = 1, material_nuclide_size(p % material) - i_nuc = material_nuclide(p % material, l) - atom_density_ = material_atom_density(p % material, l) - score = score + micro_xs(i_nuc) % reaction(m) * atom_density_ * flux - end do - end if - end if - end if - - case default - if (t % estimator() == ESTIMATOR_ANALOG) then - ! Any other score is assumed to be a MT number. Thus, we just need - ! to check if it matches the MT number of the event - if (p % event_MT /= score_bin) cycle SCORE_LOOP - score = p % last_wgt * flux - - else - ! Any other cross section has to be calculated on-the-fly. For - ! cross sections that are used often (e.g. n2n, ngamma, etc. for - ! depletion), it might make sense to optimize this section or - ! pre-calculate cross sections - if (score_bin > 1) then - ! Set default score - score = ZERO - - if (i_nuclide >= 0) then - m = nuclides(i_nuclide+1) % reaction_index(score_bin) - if (m /= 0) then - ! Retrieve temperature and energy grid index and interpolation - ! factor - i_temp = micro_xs(i_nuclide+1) % index_temp + 1 - if (i_temp > 0) then - i_energy = micro_xs(i_nuclide+1) % index_grid - f = micro_xs(i_nuclide+1) % interp_factor - - associate (rx => nuclides(i_nuclide+1) % reactions(m)) - threshold = rx % xs_threshold(i_temp) - if (i_energy >= threshold) then - score = ((ONE - f) * rx % xs(i_temp, i_energy - & - threshold + 1) + f * rx % xs(i_temp, i_energy - & - threshold + 2)) * atom_density * flux - end if - end associate - else - ! This block is reached if multipole is turned on and we're in - ! the resolved range. Assume xs is zero. - score = ZERO - end if - end if - - else - if (p % material /= MATERIAL_VOID) then - do l = 1, material_nuclide_size(p % material) - ! Get atom density - atom_density_ = material_atom_density(p % material, l) - - ! Get index in nuclides array - i_nuc = material_nuclide(p % material, l) - - m = nuclides(i_nuc) % reaction_index(score_bin) - if (m /= 0) then - ! Retrieve temperature and energy grid index and - ! interpolation factor - i_temp = micro_xs(i_nuc) % index_temp + 1 - if (i_temp > 0) then - i_energy = micro_xs(i_nuc) % index_grid - f = micro_xs(i_nuc) % interp_factor - - associate (rx => nuclides(i_nuc) % reactions(m)) - threshold = rx % xs_threshold(i_temp) - if (i_energy >= threshold) then - score = score + ((ONE - f) * rx % xs(i_temp, i_energy - & - threshold + 1) + f * rx % xs(i_temp, i_energy - & - threshold + 2)) * atom_density_ * flux - end if - end associate - else - ! This block is reached if multipole is turned on and - ! we're in the resolved range. Assume xs is zero. - score = ZERO - end if - end if - end do - end if - end if - - else - call fatal_error("Invalid score type on tally " & - // to_str(t % id()) // ".") - end if - end if - - end select - - !######################################################################### - ! Add derivative information on score for differential tallies. - - if (t % deriv() /= C_NONE) then - call apply_derivative_to_score(p, i_tally, i_nuclide, atom_density, & - score_bin, score) - end if - - !######################################################################### - ! 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_ce - subroutine score_general_mg(p, i_tally, start_index, filter_index, i_nuclide, & atom_density, flux) bind(C) type(Particle), intent(in) :: p diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 6e205daff2..b48fd46216 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -848,7 +848,7 @@ score_fission_eout(Particle* p, int i_tally, int i_score, int score_bin) } extern "C" void -score_general_ce_c(Particle* p, int i_tally, int start_index, int filter_index, +score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, int i_nuclide, double atom_density, double flux) { //TODO: off-by-one @@ -1287,8 +1287,459 @@ score_general_ce_c(Particle* p, int i_tally, int start_index, int filter_index, break; + case SCORE_DECAY_RATE: + if (simulation::material_xs.absorption == 0) continue; + 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 + const auto& nuc {*data::nuclides[p->event_nuclide-1]}; + if (simulation::micro_xs[p->event_nuclide-1].absorption > 0 + && nuc.fissionable_) { + const auto& rxn {*nuc.fission_rx_[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]; + auto yield + = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d); + auto rate = rxn.products_[d].decay_rate_; + score = p->absorb_wgt * yield + * simulation::micro_xs[p->event_nuclide-1].fission + / simulation::micro_xs[p->event_nuclide-1].absorption + * rate * flux; + 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.; + // We need to be careful not to overshoot the number of + // delayed groups since this could cause the range of the + // rxn.products_ array to be exceeded. Hence, we use the size + // of this array and not the MAX_DELAYED_GROUPS constant for + // this loop. + for (auto d = 0; d < rxn.products_.size() - 2; ++d) + { + auto yield + = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d+1); + auto rate = rxn.products_[d+1].decay_rate_; + score += rate * p->absorb_wgt + * simulation::micro_xs[p->event_nuclide-1].fission * yield + / simulation::micro_xs[p->event_nuclide-1].absorption * flux; + } + } + } + } 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) { + const auto& nuc {*data::nuclides[p->event_nuclide-1]}; + const auto& rxn {*nuc.fission_rx_[0]}; + auto rate = rxn.products_[g].decay_rate_; + //score += simulation::keff * bank.wgt * rate * flux; + 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.; + } + } + } + } + } else { + if (i_nuclide >= 0) { + const auto& nuc {*data::nuclides[i_nuclide]}; + const auto& rxn {*nuc.fission_rx_[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]; + auto yield + = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d); + auto rate = rxn.products_[d].decay_rate_; + score = simulation::micro_xs[i_nuclide].fission * yield * flux + * atom_density * rate; + score_fission_delayed_dg(i_tally, d_bin+1, score, score_index+1); + } + continue; + } else { + score = 0.; + // We need to be careful not to overshoot the number of + // delayed groups since this could cause the range of the + // rxn.products_ array to be exceeded. Hence, we use the size + // of this array and not the MAX_DELAYED_GROUPS constant for + // this loop. + for (auto d = 0; d < rxn.products_.size() - 2; ++d) + { + auto yield + = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d+1); + auto rate = rxn.products_[d+1].decay_rate_; + score += simulation::micro_xs[i_nuclide].fission * flux + * yield * atom_density * rate; + } + } + } 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())}; + if (p->material != MATERIAL_VOID) { + const Material& material {*model::materials[p->material-1]}; + for (auto i = 0; i < material.nuclide_.size(); ++i) { + auto j_nuclide = material.nuclide_[i]; + auto atom_density = material.atom_density_(i); + const auto& nuc {*data::nuclides[j_nuclide]}; + if (nuc.fissionable_) { + const auto& rxn {*nuc.fission_rx_[0]}; + // Tally each delayed group bin individually + for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { + auto d = filt.groups_[d_bin]; + auto yield + = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d); + auto rate = rxn.products_[d].decay_rate_; + score = simulation::micro_xs[j_nuclide].fission * yield + * flux * atom_density * rate; + score_fission_delayed_dg(i_tally, d_bin+1, score, + score_index+1); + } + } + } + } + continue; + } else { + score = 0.; + if (p->material != MATERIAL_VOID) { + const Material& material {*model::materials[p->material-1]}; + for (auto i = 0; i < material.nuclide_.size(); ++i) { + auto j_nuclide = material.nuclide_[i]; + auto atom_density = material.atom_density_(i); + const auto& nuc {*data::nuclides[j_nuclide]}; + if (nuc.fissionable_) { + const auto& rxn {*nuc.fission_rx_[0]}; + // We need to be careful not to overshoot the number of + // delayed groups since this could cause the range of the + // rxn.products_ array to be exceeded. Hence, we use the size + // of this array and not the MAX_DELAYED_GROUPS constant for + // this loop. + for (auto d = 0; d < rxn.products_.size() - 2; ++d) + { + auto yield + = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d+1); + auto rate = rxn.products_[d+1].decay_rate_; + score += simulation::micro_xs[j_nuclide].fission + * yield * atom_density * flux * rate; + } + } + } + } + } + } + } + break; + + + case SCORE_KAPPA_FISSION: + if (simulation::material_xs.absorption == 0.) continue; + score = 0.; + // Kappa-fission values are determined from the Q-value listed for the + // fission cross section. + 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 + const auto& nuc {*data::nuclides[p->event_nuclide-1]}; + if (simulation::micro_xs[p->event_nuclide-1].absorption > 0 + && nuc.fissionable_) { + const auto& rxn {*nuc.fission_rx_[0]}; + score = p->absorb_wgt * rxn.q_value_ + * simulation::micro_xs[p->event_nuclide-1].fission + / simulation::micro_xs[p->event_nuclide-1].absorption * 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 + const auto& nuc {*data::nuclides[p->event_nuclide-1]}; + if (simulation::micro_xs[p->event_nuclide-1].absorption > 0 + && nuc.fissionable_) { + const auto& rxn {*nuc.fission_rx_[0]}; + score = p->last_wgt * rxn.q_value_ + * simulation::micro_xs[p->event_nuclide-1].fission + / simulation::micro_xs[p->event_nuclide-1].absorption * flux; + } + } + } else { + if (i_nuclide >= 0) { + const auto& nuc {*data::nuclides[i_nuclide]}; + if (nuc.fissionable_) { + const auto& rxn {*nuc.fission_rx_[0]}; + score = rxn.q_value_ * simulation::micro_xs[i_nuclide].fission + * atom_density * flux; + } + } else { + if (p->material != MATERIAL_VOID) { + const Material& material {*model::materials[p->material-1]}; + for (auto i = 0; i < material.nuclide_.size(); ++i) { + auto j_nuclide = material.nuclide_[i]; + auto atom_density = material.atom_density_(i); + const auto& nuc {*data::nuclides[j_nuclide]}; + if (nuc.fissionable_) { + const auto& rxn {*nuc.fission_rx_[0]}; + score += rxn.q_value_ * simulation::micro_xs[j_nuclide].fission + * atom_density * flux; + } + } + } + } + } + break; + + + case SCORE_EVENTS: + // Simply count the number of scoring events + score = 1.; + break; + + + case ELASTIC: + if (tally.estimator_ == ESTIMATOR_ANALOG) { + // Check if event MT matches + if (p->event_MT != ELASTIC) continue; + score = p->last_wgt * flux; + } else { + if (i_nuclide >= 0) { + if (simulation::micro_xs[i_nuclide].elastic == CACHE_INVALID) + data::nuclides[i_nuclide]->calculate_elastic_xs(); + score = simulation::micro_xs[i_nuclide].elastic * atom_density * flux; + } else { + score = 0.; + if (p->material != MATERIAL_VOID) { + const Material& material {*model::materials[p->material-1]}; + for (auto i = 0; i < material.nuclide_.size(); ++i) { + auto j_nuclide = material.nuclide_[i]; + auto atom_density = material.atom_density_(i); + if (simulation::micro_xs[j_nuclide].elastic == CACHE_INVALID) + data::nuclides[j_nuclide]->calculate_elastic_xs(); + score += simulation::micro_xs[j_nuclide].elastic * atom_density + * flux; + } + } + } + } + break; + + case SCORE_FISS_Q_PROMPT: + case SCORE_FISS_Q_RECOV: + //continue; + if (simulation::material_xs.absorption == 0.) continue; + score = 0.; + 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 + const auto& nuc {*data::nuclides[p->event_nuclide-1]}; + if (simulation::micro_xs[p->event_nuclide-1].absorption > 0) { + double q_value = 0.; + if (score_bin == SCORE_FISS_Q_PROMPT) { + if (nuc.fission_q_prompt_) + q_value = (*nuc.fission_q_prompt_)(p->last_E); + } else if (score_bin == SCORE_FISS_Q_RECOV) { + if (nuc.fission_q_recov_) + q_value = (*nuc.fission_q_recov_)(p->last_E); + } + score = p->absorb_wgt * q_value + * simulation::micro_xs[p->event_nuclide-1].fission + / simulation::micro_xs[p->event_nuclide-1].absorption * 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 + const auto& nuc {*data::nuclides[p->event_nuclide-1]}; + if (simulation::micro_xs[p->event_nuclide-1].absorption > 0) { + double q_value = 0.; + if (score_bin == SCORE_FISS_Q_PROMPT) { + if (nuc.fission_q_prompt_) + q_value = (*nuc.fission_q_prompt_)(p->last_E); + } else if (score_bin == SCORE_FISS_Q_RECOV) { + if (nuc.fission_q_recov_) + q_value = (*nuc.fission_q_recov_)(p->last_E); + } + score = p->last_wgt * q_value + * simulation::micro_xs[p->event_nuclide-1].fission + / simulation::micro_xs[p->event_nuclide-1].absorption * flux; + } + } + } else { + if (i_nuclide >= 0) { + const auto& nuc {*data::nuclides[i_nuclide]}; + double q_value = 0.; + if (score_bin == SCORE_FISS_Q_PROMPT) { + if (nuc.fission_q_prompt_) + q_value = (*nuc.fission_q_prompt_)(p->last_E); + } else if (score_bin == SCORE_FISS_Q_RECOV) { + if (nuc.fission_q_recov_) + q_value = (*nuc.fission_q_recov_)(p->last_E); + } + score = q_value * simulation::micro_xs[i_nuclide].fission + * atom_density * flux; + } else { + if (p->material != MATERIAL_VOID) { + const Material& material {*model::materials[p->material-1]}; + for (auto i = 0; i < material.nuclide_.size(); ++i) { + auto j_nuclide = material.nuclide_[i]; + auto atom_density = material.atom_density_(i); + const auto& nuc {*data::nuclides[j_nuclide]}; + double q_value = 0.; + if (score_bin == SCORE_FISS_Q_PROMPT) { + if (nuc.fission_q_prompt_) + q_value = (*nuc.fission_q_prompt_)(p->last_E); + } else if (score_bin == SCORE_FISS_Q_RECOV) { + if (nuc.fission_q_recov_) + q_value = (*nuc.fission_q_recov_)(p->last_E); + } + score += q_value * simulation::micro_xs[j_nuclide].fission + * atom_density * flux; + } + } + } + } + break; + + + case N_2N: + case N_3N: + case N_4N: + case N_GAMMA: + case N_P: + case N_A: + if (tally.estimator_ == ESTIMATOR_ANALOG) { + // Check if the event MT matches + if (p->event_MT != score_bin) continue; + score = p->last_wgt * flux; + } else { + int m; + switch (score_bin) { + case N_GAMMA: m = 0; break; + case N_P: m = 1; break; + case N_A: m = 2; break; + case N_2N: m = 3; break; + case N_3N: m = 4; break; + case N_4N: m = 5; break; + } + if (i_nuclide >= 0) { + score = simulation::micro_xs[i_nuclide].reaction[m] * atom_density + * flux; + } else { + score = 0.; + if (p->material != MATERIAL_VOID) { + const Material& material {*model::materials[p->material-1]}; + for (auto i = 0; i < material.nuclide_.size(); ++i) { + auto j_nuclide = material.nuclide_[i]; + auto atom_density = material.atom_density_(i); + score += simulation::micro_xs[j_nuclide].reaction[m] + * atom_density * flux; + } + } + } + } + break; + + default: - continue; + if (tally.estimator_ == ESTIMATOR_ANALOG) { + // Any other score is assumed to be a MT number. Thus, we just need + // to check if it matches the MT number of the event + if (p->event_MT != score_bin) continue; + score = p->last_wgt*flux; + } else { + // Any other cross section has to be calculated on-the-fly + if (score_bin < 2) fatal_error("Invalid score type on tally " + + std::to_string(tally.id_)); + score = 0.; + if (i_nuclide >= 0) { + const auto& nuc {*data::nuclides[i_nuclide]}; + auto m = nuc.reaction_index_[score_bin]; + if (m == C_NONE) continue; + const auto& rxn {*nuc.reactions_[m]}; + auto i_temp = simulation::micro_xs[i_nuclide].index_temp; + if (i_temp >= 0) { // Can be false due to multipole + auto i_grid = simulation::micro_xs[i_nuclide].index_grid - 1; + auto f = simulation::micro_xs[i_nuclide].interp_factor; + const auto& xs {rxn.xs_[i_temp]}; + auto threshold = xs.threshold - 1; + if (i_grid >= xs.threshold) { + score = ((1.0 - f) * xs.value[i_grid-threshold] + + f * xs.value[i_grid-threshold+1]) * atom_density * flux; + } + } + } else { + if (p->material != MATERIAL_VOID) { + const Material& material {*model::materials[p->material-1]}; + for (auto i = 0; i < material.nuclide_.size(); ++i) { + auto j_nuclide = material.nuclide_[i]; + auto atom_density = material.atom_density_(i); + const auto& nuc {*data::nuclides[j_nuclide]}; + auto m = nuc.reaction_index_[score_bin]; + if (m == C_NONE) continue; + const auto& rxn {*nuc.reactions_[m]}; + auto i_temp = simulation::micro_xs[j_nuclide].index_temp; + if (i_temp >= 0) { // Can be false due to multipole + auto i_grid = simulation::micro_xs[j_nuclide].index_grid - 1; + auto f = simulation::micro_xs[j_nuclide].interp_factor; + const auto& xs {rxn.xs_[i_temp]}; + auto threshold = xs.threshold - 1; + if (i_grid >= threshold) { + score += ((1.0 - f) * xs.value[i_grid-threshold] + + f * xs.value[i_grid-threshold+1]) * atom_density + * flux; + } + } + } + } + } + } } // Add derivative information on score for differnetial tallies.