From 355ce75f0ac38456c3e63b9d65f4b4cac307aff3 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 6 Feb 2019 13:39:57 -0500 Subject: [PATCH] Start moving score_general_ce to C++ --- src/tallies/tally.F90 | 473 +++--------------------------------------- src/tallies/tally.cpp | 460 +++++++++++++++++++++++++++++++++++++++- 2 files changed, 484 insertions(+), 449 deletions(-) diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 4dd2d9ae24..cc3b169e42 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -47,6 +47,18 @@ module tally end interface interface + subroutine score_general_ce_c(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 + end subroutine + subroutine score_fission_delayed_dg(i_tally, d_bin, score, score_index) bind(C) import C_INT, C_DOUBLE integer(C_INT), value :: i_tally @@ -154,6 +166,9 @@ contains 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 @@ -174,480 +189,42 @@ contains case (SCORE_FLUX) - if (t % estimator() == ESTIMATOR_ANALOG) then - ! All events score to a flux bin. We actually use a collision - ! estimator in place of an analog one since there is no way to count - ! 'events' exactly for the flux - if (survival_biasing) then - ! We need to account for the fact that some weight was already - ! absorbed - score = p % last_wgt + p % absorb_wgt - else - score = p % last_wgt - end if - - if (p % type == NEUTRON .or. p % type == PHOTON) then - score = score / material_xs % total * flux - else - score = ZERO - end if - - else - ! For flux, we need no cross section - score = flux - end if + cycle SCORE_LOOP case (SCORE_TOTAL) - if (t % estimator() == ESTIMATOR_ANALOG) then - ! All events will score to the total reaction rate. We can just - ! use the weight of the particle entering the collision as the - ! score - if (survival_biasing) then - ! We need to account for the fact that some weight was already - ! absorbed - score = (p % last_wgt + p % absorb_wgt) * flux - else - score = p % last_wgt * flux - end if - - else - if (i_nuclide >= 0) then - score = micro_xs(i_nuclide+1) % total * atom_density * flux - else - score = material_xs % total * flux - end if - end if + cycle SCORE_LOOP case (SCORE_INVERSE_VELOCITY) - if (t % estimator() == ESTIMATOR_ANALOG) then - ! All events score to an inverse velocity bin. We actually use a - ! collision estimator in place of an analog one since there is no way - ! to count 'events' exactly for the inverse velocity - if (survival_biasing) then - ! We need to account for the fact that some weight was already - ! absorbed - score = p % last_wgt + p % absorb_wgt - else - score = p % last_wgt - end if - - ! Score the flux weighted inverse velocity with velocity in units of - ! cm/s - score = score / material_xs % total & - / (sqrt(TWO * E / MASS_NEUTRON_EV) * C_LIGHT * 100.0_8) * flux - - else - ! For inverse velocity, we don't need a cross section. The velocity is - ! in units of cm/s. - score = flux / (sqrt(TWO * E / MASS_NEUTRON_EV) * C_LIGHT * 100.0_8) - end if + cycle SCORE_LOOP case (SCORE_SCATTER) - if (t % estimator() == ESTIMATOR_ANALOG) then - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP - ! Since only scattering events make it here, again we can use - ! the weight entering the collision as the estimator for the - ! reaction rate - score = p % last_wgt * flux - - else - if (i_nuclide >= 0) then - score = (micro_xs(i_nuclide+1) % total & - - micro_xs(i_nuclide+1) % absorption) * atom_density * flux - else - score = (material_xs % total - material_xs % absorption) * flux - end if - end if + cycle SCORE_LOOP case (SCORE_NU_SCATTER) - ! Only analog estimators are available. - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP - ! For scattering production, we need to use the pre-collision weight - ! times the yield as the estimate for the number of neutrons exiting a - ! reaction with neutrons in the exit channel - if (p % event_MT == ELASTIC .or. p % event_MT == N_LEVEL .or. & - (p % event_MT >= N_N1 .and. p % event_MT <= N_NC)) then - ! Don't waste time on very common reactions we know have - ! multiplicities of one. - score = p % last_wgt * flux - else - m = nuclides(p % event_nuclide) % reaction_index(p % event_MT) - - ! Get yield and apply to score - associate (rxn => nuclides(p % event_nuclide) % reactions(m)) - score = p % last_wgt * flux & - * rxn % product_yield(1, E) - end associate - end if + cycle SCORE_LOOP case (SCORE_ABSORPTION) - if (t % estimator() == ESTIMATOR_ANALOG) then - if (survival_biasing) then - ! No absorption events actually occur if survival biasing is on -- - ! just use weight absorbed in survival biasing - score = p % absorb_wgt * flux - else - ! Skip any event where the particle wasn't absorbed - if (p % event == EVENT_SCATTER) cycle SCORE_LOOP - ! All fission and absorption events will contribute here, so we - ! can just use the particle's weight entering the collision - score = p % last_wgt * flux - end if - - else - if (i_nuclide >= 0) then - score = micro_xs(i_nuclide+1) % absorption * atom_density * flux - else - score = material_xs % absorption * flux - end if - end if + cycle SCORE_LOOP case (SCORE_FISSION) - - if (material_xs % absorption == ZERO) cycle SCORE_LOOP - - 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 - if (micro_xs(p % event_nuclide) % absorption > ZERO) then - score = p % absorb_wgt * micro_xs(p % event_nuclide) % fission & - / micro_xs(p % event_nuclide) % absorption * flux - else - score = ZERO - end if - 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 * micro_xs(p % event_nuclide) % fission & - / micro_xs(p % event_nuclide) % absorption * flux - end if - - else - if (i_nuclide >= 0) then - score = micro_xs(i_nuclide+1) % fission * atom_density * flux - else - score = material_xs % fission * flux - end if - end if + cycle SCORE_LOOP case (SCORE_NU_FISSION) - - if (material_xs % absorption == ZERO) cycle SCORE_LOOP - - 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 (micro_xs(p % event_nuclide) % absorption > ZERO) then - score = p % absorb_wgt * micro_xs(p % event_nuclide) % & - nu_fission / micro_xs(p % event_nuclide) % absorption * flux - else - score = ZERO - 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 * flux - end if - - else - if (i_nuclide >= 0) then - score = micro_xs(i_nuclide+1) % nu_fission * atom_density * flux - else - score = material_xs % nu_fission * flux - end if - end if + cycle SCORE_LOOP case (SCORE_PROMPT_NU_FISSION) - - if (material_xs % absorption == ZERO) cycle SCORE_LOOP - - 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 - ! prompt-nu-fission - if (micro_xs(p % event_nuclide) % absorption > ZERO) then - score = p % absorb_wgt * micro_xs(p % event_nuclide) % fission & - * nuclides(p % event_nuclide) % nu(E, EMISSION_PROMPT) & - / micro_xs(p % event_nuclide) % absorption * flux - else - score = ZERO - 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 as prompt neutrons. 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 - end if - - else - if (i_nuclide >= 0) then - score = micro_xs(i_nuclide+1) % fission * nuclides(i_nuclide+1) % & - nu(E, EMISSION_PROMPT) * atom_density * flux - 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) - - ! Accumulate the contribution from each nuclide - score = score + micro_xs(i_nuc) % fission * nuclides(i_nuc) % & - nu(E, EMISSION_PROMPT) * atom_density_ * flux - end do - end if - end if - end if + cycle SCORE_LOOP case (SCORE_DELAYED_NU_FISSION) - - if (material_xs % absorption == ZERO) cycle SCORE_LOOP - - ! 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 - ! 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) - - ! Compute the score and tally to bin - score = p % absorb_wgt * yield & - * micro_xs(p % event_nuclide) % fission & - / micro_xs(p % event_nuclide) % absorption * flux - 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 multiplying the absorbed weight by the fraction of the - ! delayed-nu-fission xs to the absorption xs - score = p % absorb_wgt * micro_xs(p % event_nuclide) % fission & - * nuclides(p % event_nuclide) % nu(E, EMISSION_DELAYED) & - / micro_xs(p % event_nuclide) % absorption * flux - 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. - - ! 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 score and tally to bin - score = keff * p % wgt_bank / p % n_bank * & - p % n_delayed_bank(d) * flux - - call score_fission_delayed_dg(i_tally, d_bin, score, score_index) - end do - cycle SCORE_LOOP - end select - else - - ! Add the contribution from all delayed groups - score = keff * p % wgt_bank / p % n_bank * & - sum(p % n_delayed_bank) * flux - end if - 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) - - ! Compute the score and tally to bin - score = micro_xs(i_nuclide+1) % fission * yield * & - atom_density * flux - 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 multiplying the delayed-nu-fission macro xs by the flux - score = micro_xs(i_nuclide+1) % fission * nuclides(i_nuclide+1) % & - nu(E, EMISSION_DELAYED) * atom_density * flux - 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) - - ! 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) - - ! Compute the score and tally to bin - score = micro_xs(i_nuc) % fission * yield & - * atom_density_ * flux - call score_fission_delayed_dg(i_tally, d_bin, score, & - score_index) - end do - 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) - - ! Accumulate the contribution from each nuclide - score = score + micro_xs(i_nuc) % fission * nuclides(i_nuc) %& - nu(E, EMISSION_DELAYED) * atom_density_ * flux - end do - end if - end if - end if - end if + cycle SCORE_LOOP case (SCORE_DECAY_RATE) diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index ea948ae401..6e205daff2 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -7,6 +7,7 @@ #include "openmc/message_passing.h" #include "openmc/mgxs_interface.h" #include "openmc/nuclide.h" +#include "openmc/reaction_product.h" #include "openmc/settings.h" #include "openmc/simulation.h" #include "openmc/tallies/derivative.h" @@ -51,6 +52,8 @@ apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide, extern "C" int energy_filter_search(const EnergyFilter* filt, double val); +extern "C" double get_precursor_decay_rate(int i_nuclide, int d); + //============================================================================== // Global variable definitions //============================================================================== @@ -814,7 +817,7 @@ score_fission_eout(Particle* p, int i_tally, int i_score, int score_bin) filter_weight *= match.weights_[i_bin-1]; } - score_fission_delayed_dg(i_tally, d_bin+1, score*filter_weight, + score_fission_delayed_dg(i_tally, d_bin+1, score*filter_weight, i_score); } } @@ -844,6 +847,461 @@ score_fission_eout(Particle* p, int i_tally, int i_score, int score_bin) simulation::filter_matches[i_eout_filt].bins_[i_bin-1] = bin_energyout; } +extern "C" void +score_general_ce_c(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); + + // Get the pre-collision energy of the particle. + auto E = p->last_E; + + for (auto i = 0; i < tally.scores_.size(); ++i) { + auto score_bin = tally.scores_[i]; + auto score_index = start_index + i; + + double score; + + //TODO: off-by-one throughout on p->event_nuclide + //TODO: off-by-one throughout on score_index in calls to score_fission_eout + //TODO: off-by-one throughout on p->material + + switch (score_bin) { + + + case SCORE_FLUX: + if (tally.estimator_ == ESTIMATOR_ANALOG) { + // All events score to a flux bin. We actually use a collision estimator + // in place of an analog one since there is no way to count 'events' + // exactly for the flux + if (settings::survival_biasing) { + // We need to account for the fact that some weight was already + // absorbed + score = p->last_wgt + p->absorb_wgt; + } else { + score = p->last_wgt; + } + + if (p->type == static_cast(ParticleType::neutron) || + p->type == static_cast(ParticleType::photon)) { + score *= flux / simulation::material_xs.total; + } else { + score = 0.; + } + } else { + score = flux; + } + break; + + + case SCORE_TOTAL: + if (tally.estimator_ == ESTIMATOR_ANALOG) { + // All events will score to the total reaction rate. We can just use + // use the weight of the particle entering the collision as the score + if (settings::survival_biasing) { + // We need to account for the fact that some weight was already + // absorbed + score = (p->last_wgt + p->absorb_wgt) * flux; + } else { + score = p->last_wgt * flux; + } + + } else { + if (i_nuclide >= 0) { + score = simulation::micro_xs[i_nuclide].total * atom_density * flux; + } else { + score = simulation::material_xs.total * flux; + } + } + break; + + + case SCORE_INVERSE_VELOCITY: + if (tally.estimator_ == ESTIMATOR_ANALOG) { + // All events score to an inverse velocity bin. We actually use a + // collision estimator in place of an analog one since there is no way + // to count 'events' exactly for the inverse velocity + if (settings::survival_biasing) { + // We need to account for the fact that some weight was already + // absorbed + score = (p->last_wgt + p->absorb_wgt); + } else { + score = p->last_wgt; + } + score *= flux / simulation::material_xs.total; + } else { + score = flux; + } + // Score inverse velocity in units of s/cm. + score /= std::sqrt(2. * E / MASS_NEUTRON_EV) * C_LIGHT * 100.; + break; + + + case SCORE_SCATTER: + if (tally.estimator_ == ESTIMATOR_ANALOG) { + // Skip any event where the particle didn't scatter + if (p->event != EVENT_SCATTER) continue; + // Since only scattering events make it here, again we can use the + // weight entering the collision as the estimator for the reaction rate + score = p->last_wgt * flux; + } else { + if (i_nuclide >= 0) { + score = (simulation::micro_xs[i_nuclide].total + - simulation::micro_xs[i_nuclide].absorption) * atom_density * flux; + } else { + score = (simulation::material_xs.total + - simulation::material_xs.absorption) * flux; + } + } + break; + + + case SCORE_NU_SCATTER: + // Only analog estimators are available. + // Skip any event where the particle didn't scatter + if (p->event != EVENT_SCATTER) continue; + // For scattering production, we need to use the pre-collision weight + // times the yield as the estimate for the number of neutrons exiting a + // reaction with neutrons in the exit channel + if (p->event_MT == ELASTIC || p->event_MT == N_LEVEL || + (p->event_MT >= N_N1 && p->event_MT <= N_NC)) { + // Don't waste time on very common reactions we know have + // multiplicities of one. + score = p->last_wgt * flux; + } else { + // Get yield and apply to score + auto m = + data::nuclides[p->event_nuclide-1]->reaction_index_[p->event_MT]; + const auto& rxn {*data::nuclides[p->event_nuclide-1]->reactions_[m]}; + score = p->last_wgt * flux * (*rxn.products_[0].yield_)(E); + } + break; + + + case SCORE_ABSORPTION: + if (tally.estimator_ == ESTIMATOR_ANALOG) { + if (settings::survival_biasing) { + // No absorption events actually occur if survival biasing is on -- + // just use weight absorbed in survival biasing + score = p->absorb_wgt * flux; + } else { + // Skip any event where the particle wasn't absorbed + if (p->event == EVENT_SCATTER) continue; + // All fission and absorption events will contribute here, so we + // can just use the particle's weight entering the collision + score = p->last_wgt * flux; + } + } else { + if (i_nuclide >= 0) { + score = simulation::micro_xs[i_nuclide].absorption * atom_density + * flux; + } else { + score = simulation::material_xs.absorption * flux; + } + } + break; + + + case SCORE_FISSION: + 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 + // fission + if (simulation::micro_xs[p->event_nuclide-1].absorption > 0) { + score = p->absorb_wgt + * simulation::micro_xs[p->event_nuclide-1].fission + / simulation::micro_xs[p->event_nuclide-1].absorption * flux; + } else { + score = 0.; + } + } 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 + * simulation::micro_xs[p->event_nuclide-1].fission + / simulation::micro_xs[p->event_nuclide-1].absorption * flux; + } + } else { + if (i_nuclide >= 0) { + score = simulation::micro_xs[i_nuclide].fission * atom_density * flux; + } else { + score = simulation::material_xs.fission * flux; + } + } + break; + + + case SCORE_NU_FISSION: + if (simulation::material_xs.absorption == 0) continue; + 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 + // nu-fission + if (simulation::micro_xs[p->event_nuclide-1].absorption > 0) { + score = p->absorb_wgt + * simulation::micro_xs[p->event_nuclide-1].nu_fission + / simulation::micro_xs[p->event_nuclide-1].absorption * flux; + } else { + score = 0.; + } + } 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. + score = simulation::keff * p->wgt_bank * flux; + } + } else { + if (i_nuclide >= 0) { + score = simulation::micro_xs[i_nuclide].nu_fission * atom_density + * flux; + } else { + score = simulation::material_xs.nu_fission * flux; + } + } + break; + + + case SCORE_PROMPT_NU_FISSION: + if (simulation::material_xs.absorption == 0) continue; + 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 + if (simulation::micro_xs[p->event_nuclide-1].absorption > 0) { + score = p->absorb_wgt + * simulation::micro_xs[p->event_nuclide-1].fission + * data::nuclides[p->event_nuclide-1] + ->nu(E, ReactionProduct::EmissionMode::prompt) + / simulation::micro_xs[p->event_nuclide-1].absorption * flux; + } else { + score = 0.; + } + } 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; + } + } else { + if (i_nuclide >= 0) { + score = simulation::micro_xs[i_nuclide].fission + * data::nuclides[i_nuclide] + ->nu(E, ReactionProduct::EmissionMode::prompt) + * atom_density * flux; + } else { + score = 0.; + // Add up contributions from each nuclide in the material. + 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].fission + * data::nuclides[j_nuclide] + ->nu(E, ReactionProduct::EmissionMode::prompt) + * atom_density * flux; + } + } + } + } + break; + + + case SCORE_DELAYED_NU_FISSION: + if (simulation::material_xs.absorption == 0) continue; + 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 (simulation::micro_xs[p->event_nuclide-1].absorption > 0 + && data::nuclides[p->event_nuclide-1]->fissionable_) { + 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 = data::nuclides[p->event_nuclide-1] + ->nu(E, ReactionProduct::EmissionMode::delayed, d); + score = p->absorb_wgt * yield + * simulation::micro_xs[p->event_nuclide-1].fission + / simulation::micro_xs[p->event_nuclide-1].absorption * 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 + score = p->absorb_wgt + * simulation::micro_xs[p->event_nuclide-1].fission + * data::nuclides[p->event_nuclide-1] + ->nu(E, ReactionProduct::EmissionMode::delayed) + / 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. + 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; + 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; + } + } + } else { + if (i_nuclide >= 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 = data::nuclides[i_nuclide] + ->nu(E, ReactionProduct::EmissionMode::delayed, d); + score = simulation::micro_xs[i_nuclide].fission * yield + * atom_density * 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 delayed-nu-fission macro xs by the flux + score = simulation::micro_xs[i_nuclide].fission + * data::nuclides[i_nuclide] + ->nu(E, ReactionProduct::EmissionMode::delayed) + * atom_density * flux; + } + } 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); + // 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 = data::nuclides[j_nuclide] + ->nu(E, ReactionProduct::EmissionMode::delayed, d); + score = simulation::micro_xs[j_nuclide].fission * yield + * atom_density * flux; + 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); + score += simulation::micro_xs[j_nuclide].fission + * data::nuclides[j_nuclide] + ->nu(E, ReactionProduct::EmissionMode::delayed) + * atom_density * flux; + } + } + } + } + } + break; + + + default: + continue; + } + + // Add derivative information on score for differnetial tallies. + if (tally.deriv_ != C_NONE) + apply_derivative_to_score(p, i_tally, i_nuclide, atom_density, score_bin, + &score); + + // Update tally results + #pragma omp atomic + results(filter_index-1, score_index, RESULT_VALUE) += score; + } +} + //! Tally rates for when the user requests a tally on all nuclides. void