diff --git a/src/constants.F90 b/src/constants.F90 index b87c461701..43ca861901 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -231,9 +231,17 @@ module constants ! Tally estimator types integer, parameter :: & - ESTIMATOR_COLLISION = 1, & + ESTIMATOR_ANALOG = 1, & ESTIMATOR_TRACKLENGTH = 2 + ! Event types for tallies + integer, parameter :: & + EVENT_SURFACE = -2, & + EVENT_LATTICE = -1, & + EVENT_SCATTER = 1, & + EVENT_ABSORB = 2, & + EVENT_FISSION = 3 + ! Tally score type integer, parameter :: N_SCORE_TYPES = 15 integer, parameter :: & diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 1d3af286ff..20503d0dc9 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -738,6 +738,14 @@ contains ! Set tally type to volume by default t % type = TALLY_VOLUME + ! It's desirable to use a track-length esimator for tallies since + ! generally more events will score to the tally, reducing the + ! variance. However, for tallies that require information on + ! post-collision parameters (e.g. tally with an energyout filter) the + ! analog esimator must be used. + + t % estimator = ESTIMATOR_TRACKLENGTH + ! Copy material id t % id = tally_(i) % id @@ -836,6 +844,9 @@ contains t % energy_out(j) = str_to_real(words(j)) end do t % n_filter_bins(FILTER_ENERGYOUT) = n_words - 1 + + ! Set tally estimator to analog + t % estimator = ESTIMATOR_ANALOG end if ! Read macro reactions @@ -863,20 +874,44 @@ contains t % score_bins(j) % scalar = SCORE_SCATTER case ('nu-scatter') t % score_bins(j) % scalar = SCORE_NU_SCATTER + + ! Set tally estimator to analog + t % estimator = ESTIMATOR_ANALOG case ('scatter-1') t % score_bins(j) % scalar = SCORE_SCATTER_1 + + ! Set tally estimator to analog + t % estimator = ESTIMATOR_ANALOG case ('scatter-2') t % score_bins(j) % scalar = SCORE_SCATTER_2 + + ! Set tally estimator to analog + t % estimator = ESTIMATOR_ANALOG case ('scatter-3') t % score_bins(j) % scalar = SCORE_SCATTER_3 + + ! Set tally estimator to analog + t % estimator = ESTIMATOR_ANALOG case ('n1n') t % score_bins(j) % scalar = SCORE_N_1N + + ! Set tally estimator to analog + t % estimator = ESTIMATOR_ANALOG case ('n2n') t % score_bins(j) % scalar = SCORE_N_2N + + ! Set tally estimator to analog + t % estimator = ESTIMATOR_ANALOG case ('n3n') t % score_bins(j) % scalar = SCORE_N_3N + + ! Set tally estimator to analog + t % estimator = ESTIMATOR_ANALOG case ('n4n') t % score_bins(j) % scalar = SCORE_N_4N + + ! Set tally estimator to analog + t % estimator = ESTIMATOR_ANALOG case ('absorption') t % score_bins(j) % scalar = SCORE_ABSORPTION if (t % n_filter_bins(FILTER_ENERGYOUT) > 0) then diff --git a/src/particle_header.F90 b/src/particle_header.F90 index 8f81fbeb1f..692f7de357 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -1,6 +1,6 @@ module particle_header - use constants, only: NEUTRON, ONE, NONE + use constants, only: NEUTRON, ONE, NONE, ZERO use geometry_header, only: BASE_UNIVERSE implicit none @@ -53,6 +53,9 @@ module particle_header real(8) :: last_wgt ! last particle weight real(8) :: last_E ! last energy + ! What event last took place + integer :: event + ! Post-collision physical data integer :: n_bank ! number of fission sites banked diff --git a/src/physics.F90 b/src/physics.F90 index 68c701ae6a..9e6bffef31 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -16,7 +16,8 @@ module physics use random_lcg, only: prn use search, only: binary_search use string, only: to_str - use tally, only: score_tally, score_surface_current + use tally, only: score_analog_tally, score_tracklength_tally, & + score_surface_current implicit none @@ -66,7 +67,6 @@ contains ! Initialize number of events to zero n_event = 0 - ! find energy index, interpolation factor do while (p % alive) ! Calculate microscopic and macroscopic cross sections -- note: if the @@ -95,15 +95,22 @@ contains coord => coord % next end do + ! Score track-length tallies + if (tallies_on) call score_tracklength_tally(p, distance) + if (d_collision > d_boundary) then last_cell = p % coord % cell p % coord % cell = NONE if (lattice_crossed /= NONE) then + ! Particle crosses lattice boundary p % surface = NONE call cross_lattice(p, lattice_crossed) + p % event = EVENT_LATTICE else + ! Particle crosses surface p % surface = surface_crossed call cross_surface(p, last_cell) + p % event = EVENT_SURFACE end if else ! collision @@ -148,10 +155,6 @@ contains type(Particle), pointer :: p - integer :: MT ! ENDF reaction number - logical :: scattered ! was this a scattering reaction? - logical :: fissioned ! was this a fission reaction? - ! Store pre-collision particle properties p % last_wgt = p % wgt p % last_E = p % E @@ -166,13 +169,13 @@ contains if (tallies_on) call score_surface_current(p) ! Sample nuclide/reaction for the material the particle is in - call sample_reaction(p, MT) + call sample_reaction(p) ! Display information about collision if (verbosity >= 10 .or. trace) then - message = " " // trim(reaction_name(MT)) // ". Energy = " // & - trim(to_str(p % E * 1e6_8)) // " eV." - call write_message() +!!$ message = " " // trim(reaction_name(MT)) // ". Energy = " // & +!!$ trim(to_str(p % E * 1e6_8)) // " eV." +!!$ call write_message() end if ! check for very low energy @@ -182,21 +185,12 @@ contains call warning() end if - ! Check if particle scattered or fissioned - if (survival_biasing) then - fissioned = .false. - scattered = .true. - else - fissioned = is_fission(MT) - scattered = is_scatter(MT) - end if - ! Score collision estimator tallies for any macro tallies -- this is done ! after a collision has occurred rather than before because we need ! information on the outgoing energy for any tallies with an outgoing energy ! filter - if (tallies_on) call score_tally(p, scattered, fissioned) + if (tallies_on) call score_analog_tally(p) ! Reset number of particles banked during collision p % n_bank = 0 @@ -211,10 +205,9 @@ contains ! disappearance are treated implicitly. !=============================================================================== - subroutine sample_reaction(p, MT) + subroutine sample_reaction(p) type(Particle), pointer :: p - integer, intent(out) :: MT ! ENDF MT number of reaction that occured integer :: i ! index over nuclides in a material integer :: index_nuclide ! index in nuclides array @@ -281,7 +274,7 @@ contains ! See if disappearance reaction happens if (prob > cutoff) then p % alive = .false. - MT = N_DISAPPEAR + p % event = EVENT_ABSORB return end if end if @@ -337,16 +330,24 @@ contains if (micro_xs(index_nuclide) % use_ptable) then + ! In the case that the particle is in the unresolved resonance + ! region and probability tables are being used, we should only + ! check the first fission reaction + prob = prob + micro_xs(index_nuclide) % fission if (prob > cutoff) then rxn => nuc % reactions(nuc % index_fission(1)) call create_fission_sites(p, index_nuclide, rxn, .true.) p % alive = .false. - MT = rxn % MT + p % event = EVENT_FISSION return end if + else + ! With no probability tables, we need to loop through each fission + ! reaction type + do i = 1, nuc % n_fission rxn => nuc % reactions(nuc % index_fission(i)) @@ -365,7 +366,7 @@ contains if (prob > cutoff) then call create_fission_sites(p, index_nuclide, rxn, .true.) p % alive = .false. - MT = rxn % MT + p % event = EVENT_FISSION return end if end do @@ -418,8 +419,6 @@ contains end if - ! Set MT to be returned - MT = 2 else ! ======================================================================= ! INELASTIC SCATTERING @@ -460,10 +459,12 @@ contains ! Perform collision physics for inelastics scattering call inelastic_scatter(p, nuc, rxn) - ! Set MT to be returned - MT = rxn % MT end if + ! If we made it this far, it means that a scattering reaction took place + ! since the absorption and fission blocks had return statements. + p % event = EVENT_SCATTER + end subroutine sample_reaction !=============================================================================== diff --git a/src/tally.F90 b/src/tally.F90 index 79e31fecdb..085e63024b 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -232,34 +232,26 @@ contains end subroutine add_map_element !=============================================================================== -! SCORE_TALLY contains the main logic for scoring user-specified tallies +! SCORE_ANALOG_TALLY keeps track of how many events occur in a specified cell, +! energy range, etc. Note that since these are "analog" tallies, they are only +! triggered at every collision, not every event !=============================================================================== - subroutine score_tally(p, scattered, fissioned) + subroutine score_analog_tally(p) type(Particle), pointer :: p - logical, intent(in) :: scattered - logical, intent(in) :: fissioned integer :: i integer :: j - integer :: k - integer :: n integer :: bins(N_FILTER_TYPES) - integer :: bin_energyout integer :: score_index - integer :: score_index0 integer :: score_bin - integer :: mesh_bin real(8) :: score real(8) :: last_wgt real(8) :: wgt real(8) :: mu - real(8) :: E_out - logical :: has_energyout_bin - logical :: analog - type(TallyObject), pointer :: t - type(StructuredMesh), pointer :: m + logical :: found_bin + type(TallyObject), pointer :: t => null() ! Copy particle's pre- and post-collision weight and angle last_wgt = p % last_wgt @@ -275,87 +267,14 @@ contains ! Surface current tallies are treated separately if (t % type == TALLY_SURFACE_CURRENT) cycle + ! Skip tallies that are handled by track-length estimators + if (t % estimator == ESTIMATOR_TRACKLENGTH) cycle + ! ======================================================================= ! DETERMINE SCORING BIN COMBINATION - ! determine mesh bin - if (t % n_filter_bins(FILTER_MESH) > 0) then - m => meshes(t % mesh) - - ! Determine if we're in the mesh first - call get_mesh_bin(m, p % coord0 % xyz, mesh_bin) - if (mesh_bin == NO_BIN_FOUND) cycle - bins(FILTER_MESH) = mesh_bin - else - bins(FILTER_MESH) = 1 - end if - - ! determine next universe bin - if (t % n_filter_bins(FILTER_UNIVERSE) > 0) then - bins(FILTER_UNIVERSE) = get_next_bin(FILTER_UNIVERSE, p % coord % universe, i) - if (bins(FILTER_UNIVERSE) == NO_BIN_FOUND) cycle - else - bins(FILTER_UNIVERSE) = 1 - end if - - ! determine next material bin - if (t % n_filter_bins(FILTER_MATERIAL) > 0) then - bins(FILTER_MATERIAL) = get_next_bin(FILTER_MATERIAL, p % material, i) - if (bins(FILTER_MATERIAL) == NO_BIN_FOUND) cycle - else - bins(FILTER_MATERIAL) = 1 - end if - - ! determine next cell bin - if (t % n_filter_bins(FILTER_CELL) > 0) then - bins(FILTER_CELL) = get_next_bin(FILTER_CELL, p % coord % cell, i) - if (bins(FILTER_CELL) == NO_BIN_FOUND) cycle - else - bins(FILTER_CELL) = 1 - end if - - ! determine next cellborn bin - if (t % n_filter_bins(FILTER_CELLBORN) > 0) then - bins(FILTER_CELLBORN) = get_next_bin(FILTER_CELLBORN, p % cell_born, i) - if (bins(FILTER_CELLBORN) == NO_BIN_FOUND) cycle - else - bins(FILTER_CELLBORN) = 1 - end if - - ! determine next surface bin - if (t % n_filter_bins(FILTER_SURFACE) > 0) then - bins(FILTER_SURFACE) = get_next_bin(FILTER_SURFACE, p % surface, i) - if (bins(FILTER_SURFACE) == NO_BIN_FOUND) cycle - else - bins(FILTER_SURFACE) = 1 - end if - - ! determine incoming energy bin - n = t % n_filter_bins(FILTER_ENERGYIN) - if (n > 0) then - ! check if energy of the particle is within energy bins - if (p % last_E < t % energy_in(1) .or. & - p % last_E > t % energy_in(n + 1)) cycle - - ! search to find incoming energy bin - bins(FILTER_ENERGYIN) = binary_search(t % energy_in, n + 1, p % last_E) - else - bins(FILTER_ENERGYIN) = 1 - end if - - ! determine outgoing energy bin - n = t % n_filter_bins(FILTER_ENERGYOUT) - if (n > 0) then - ! check if energy of the particle is within energy bins - if (p % E < t % energy_out(1) .or. p % E > t % energy_out(n + 1)) cycle - - ! search to find incoming energy bin - bins(FILTER_ENERGYOUT) = binary_search(t % energy_out, n + 1, p % E) - has_energyout_bin = .true. - else - bins(FILTER_ENERGYOUT) = 1 - has_energyout_bin = .false. - end if + call get_scoring_bins(p, i, bins, found_bin) + if (.not. found_bin) cycle ! ======================================================================= ! CALCULATE SCORES AND ACCUMULATE TALLY @@ -372,130 +291,152 @@ contains ! determine what type of score bin score_bin = t % score_bins(j) % scalar - ! determine if we need outgoing angle - analog = (score_bin == SCORE_NU_SCATTER .or. & - score_bin == SCORE_SCATTER_1 .or. score_bin == SCORE_SCATTER_2 .or. & - score_bin == SCORE_SCATTER_3 .or. score_bin == SCORE_N_1N .or. & - score_bin == SCORE_N_2N .or. score_bin == SCORE_N_3N .or. & - score_bin == SCORE_N_4N) + select case (score_bin) + case (SCORE_FLUX) + ! All events score to a flux bin. We actually use a collision + ! estimator since there is no way to count 'events' exactly for the + ! flux - if (has_energyout_bin .or. analog) then - ! If this tally has an outgoing energy filter, the only supported - ! reaction is scattering or nu-fission. Note that some scoring - ! quantities can only be scored in analog such as scattering - ! moments and (n,xn) + score = last_wgt / material_xs % total - if ((.not. scattered) .and. (.not. fissioned)) cycle + case (SCORE_TOTAL) + ! 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 (scattered) then - ! since scattering has already occured, we do not need to - ! multiply by the scattering cross section + score = last_wgt - select case (score_bin) - case (SCORE_SCATTER) - score = last_wgt - case (SCORE_NU_SCATTER) - score = wgt - case (SCORE_SCATTER_1) - score = last_wgt * mu - case (SCORE_SCATTER_2) - score = last_wgt * 0.5*(3.0*mu*mu - ONE) - case (SCORE_SCATTER_3) - score = last_wgt * 0.5*(5.0*mu*mu*mu - 3.0*mu) - case (SCORE_N_1N) - if (wgt == last_wgt) then - score = last_wgt - else - cycle - end if - case (SCORE_N_2N) - if (int(wgt/last_wgt) == 2) then - score = last_wgt - else - cycle - end if - case (SCORE_N_3N) - if (int(wgt/last_wgt) == 3) then - score = last_wgt - else - cycle - end if - case (SCORE_N_4N) - if (int(wgt/last_wgt) == 4) then - score = last_wgt - else - cycle - end if - case (SCORE_NU_FISSION) - cycle - case default - message = "Invalid score type on analog tally." - call fatal_error() - end select - end if + case (SCORE_SCATTER) + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle - if (fissioned) then + ! Since only scattering events make it here, again we can use the + ! weight entering the collision as the estimator for the reaction + ! rate - if (score_bin /= SCORE_NU_FISSION) cycle + score = last_wgt + case (SCORE_NU_SCATTER) + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle + + ! For scattering production, we need to use the post-collision + ! weight as the estimate for the number of neutrons exiting a + ! reaction with neutrons in the exit channel + + score = wgt + + case (SCORE_SCATTER_1) + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle + + ! The first scattering moment can be determined by using the rate + ! of scattering reactions multiplied by the cosine of the change in + ! neutron's angle due to the collision + + score = last_wgt * mu + + case (SCORE_SCATTER_2) + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle + + ! The second scattering moment can be determined in a similar + ! manner to the first scattering moment + + score = last_wgt * 0.5*(3.0*mu*mu - ONE) + + case (SCORE_SCATTER_3) + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle + + ! The first scattering moment can be determined by using the rate + ! of scattering reactions multiplied by the cosine of the change in + ! neutron's angle due to the collision + + score = last_wgt * 0.5*(5.0*mu*mu*mu - 3.0*mu) + + case (SCORE_N_1N) + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle + + ! Skip any events where weight of particle changed + if (wgt /= last_wgt) cycle + + ! All events that reach this point are (n,1n) reactions + score = last_wgt + + case (SCORE_N_2N) + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle + + ! Skip any scattering events where the weight didn't double + if (nint(wgt/last_wgt) /= 2) cycle + + ! All events that reach this point are (n,2n) reactions + score = last_wgt + + case (SCORE_N_3N) + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle + + ! Skip any scattering events where the weight didn't double + if (nint(wgt/last_wgt) /= 3) cycle + + ! All events that reach this point are (n,3n) reactions + score = last_wgt + + case (SCORE_N_4N) + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle + + ! Skip any scattering events where the weight didn't double + if (nint(wgt/last_wgt) /= 4) cycle + + ! All events that reach this point are (n,3n) reactions + score = last_wgt + + case (SCORE_ABSORPTION) + ! Skip any event where the particle wasn't absorbed + if (p % event == EVENT_SCATTER) cycle + + ! All fission and absorption events will contribute here, so we can + ! just use the particle's weight entering the collision + + score = last_wgt + + case (SCORE_FISSION) + ! Skip any non-fission events + if (p % event /= EVENT_FISSION) cycle + + ! 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 = last_wgt + + case (SCORE_NU_FISSION) + ! Skip any non-fission events + if (p % event /= EVENT_FISSION) cycle + + if (t % n_filter_bins(FILTER_ENERGYOUT) > 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 scores to multiple bins + ! neutrons were emitted with different energies, multiple outgoing + ! energy bins may have been scored to. The following logic treats + ! this special case and scores to multiple bins - ! save original outgoing energy bin and score index - bin_energyout = bins(FILTER_ENERGYOUT) - score_index0 = score_index - - ! Since the creation of fission sites is weighted such that it - ! is expected to create n_particles sites, we need to multiply - ! the score by keff to get the true nu-fission rate. Otherwise, - ! the sum of all nu-fission rates would be ~1.0. - - score = keff - - ! loop over number of particles banked - do k = 1, p % n_bank - ! determine outgoing energy from fission bank - E_out = fission_bank(n_bank - p % n_bank + k) % E - - ! change outgoing energy bin - bins(FILTER_ENERGYOUT) = binary_search(t % energy_out, n + 1, E_out) - - ! determine scoring index - score_index = sum((bins - 1) * t % stride) + 1 - - ! Add score to tally - call add_to_score(t % scores(score_index, j), score) - end do - - ! reset outgoing energy bin and score index - bins(FILTER_ENERGYOUT) = bin_energyout - score_index = score_index0 + call score_fission_eout(p, t, bins, j) cycle - end if - else - ! For tallies with no outgoing energy filter, the score is - ! calculated normally depending on the quantity specified + else - select case (score_bin) - case (SCORE_FLUX) - score = last_wgt / material_xs % total - case (SCORE_TOTAL) - score = last_wgt - case (SCORE_SCATTER) - score = last_wgt * (material_xs % total - material_xs % absorption) & - / material_xs % total - case (SCORE_ABSORPTION) - score = last_wgt * material_xs % absorption / material_xs % total - case (SCORE_FISSION) - score = last_wgt * material_xs % fission / material_xs % total - case (SCORE_NU_FISSION) - score = last_wgt * material_xs % nu_fission / material_xs % total - end select - end if + score = keff * p % n_bank + + end if + + case default + message = "Invalid score type on tally " // to_str(t % id) // "." + call fatal_error() + end select ! Add score to tally call add_to_score(t % scores(score_index, j), score) @@ -507,7 +448,270 @@ contains end do - end subroutine score_tally + end subroutine score_analog_tally + +!=============================================================================== +! SCORE_FISSION_EOUT +!=============================================================================== + + subroutine score_fission_eout(p, t, bins, j) + + type(Particle), pointer :: p + type(TallyObject), pointer :: t + integer, intent(inout) :: bins(N_FILTER_TYPES) + integer, intent(in) :: j + + integer :: k + integer :: bin_energyout + integer :: score_index + real(8) :: score + real(8) :: E_out + + ! save original outgoing energy bin and score index + bin_energyout = bins(FILTER_ENERGYOUT) + + ! Since the creation of fission sites is weighted such that it is + ! expected to create n_particles sites, we need to multiply the + ! score by keff to get the true nu-fission rate. Otherwise, the sum + ! of all nu-fission rates would be ~1.0. + + score = keff + + ! loop over number of particles banked + do k = 1, p % n_bank + ! determine outgoing energy from fission bank + E_out = fission_bank(n_bank - p % n_bank + k) % E + + ! change outgoing energy bin + bins(FILTER_ENERGYOUT) = binary_search(t % energy_out, & + size(t % energy_out), E_out) + + ! determine scoring index + score_index = sum((bins - 1) * t % stride) + 1 + + ! Add score to tally + call add_to_score(t % scores(score_index, j), score) + end do + + ! reset outgoing energy bin and score index + bins(FILTER_ENERGYOUT) = bin_energyout + + end subroutine score_fission_eout + +!=============================================================================== +! SCORE_TRACKLENGTH_TALLY +!=============================================================================== + + subroutine score_tracklength_tally(p, distance) + + type(Particle), pointer :: p + real(8), intent(in) :: distance + + integer :: i + integer :: j + integer :: bins(N_FILTER_TYPES) + integer :: score_index + integer :: score_bin + real(8) :: flux + real(8) :: score + logical :: found_bin + type(TallyObject), pointer :: t => null() + + ! Determine track-length estimate of flux + flux = p % wgt * distance + + ! A loop over all tallies is necessary because we need to simultaneously + ! determine different filter bins for the same tally in order to score to it + + do i = 1, n_tallies + t => tallies(i) + + ! Surface current tallies are treated separately + if (t % type == TALLY_SURFACE_CURRENT) cycle + + ! Skip tallies that are handled by track-length estimators + if (t % estimator == ESTIMATOR_ANALOG) cycle + + ! ======================================================================= + ! DETERMINE SCORING BIN COMBINATION + + call get_scoring_bins(p, i, bins, found_bin) + if (.not. found_bin) cycle + + ! ======================================================================= + ! CALCULATE SCORES AND ACCUMULATE TALLY + + ! If we have made it here, we have a scoring combination of bins for this + ! tally -- now we need to determine where in the scores array we should + ! be accumulating the tally values + + ! Determine scoring index for this filter combination + score_index = sum((bins - 1) * t % stride) + 1 + + ! Determine score for each bin + do j = 1, t % n_score_bins + ! determine what type of score bin + score_bin = t % score_bins(j) % scalar + + ! Determine cross section + select case(score_bin) + case (SCORE_FLUX) + score = flux + case (SCORE_TOTAL) + score = material_xs % total * flux + case (SCORE_SCATTER) + score = (material_xs % total - material_xs % absorption) * flux + case (SCORE_ABSORPTION) + score = material_xs % absorption * flux + case (SCORE_FISSION) + score = material_xs % fission * flux + case (SCORE_NU_FISSION) + score = material_xs % nu_fission * flux + case default + message = "Invalid score type on tally " // to_str(t % id) // "." + call fatal_error() + end select + + ! Add score to tally + call add_to_score(t % scores(score_index, j), score) + + end do + + ! Reset tally map positioning + position = 0 + + end do + + end subroutine score_tracklength_tally + +!=============================================================================== +! GET_SCORING_BINS +!=============================================================================== + + subroutine get_scoring_bins(p, index_tally, bins, found_bin) + + type(Particle), pointer :: p + integer, intent(in) :: index_tally + integer, intent(out) :: bins(N_FILTER_TYPES) + logical, intent(out) :: found_bin + + integer :: n + integer :: mesh_bin + type(TallyObject), pointer :: t => null() + type(StructuredMesh), pointer :: m => null() + + found_bin = .true. + t => tallies(index_tally) + + ! determine mesh bin + if (t % n_filter_bins(FILTER_MESH) > 0) then + m => meshes(t % mesh) + + ! Determine if we're in the mesh first + call get_mesh_bin(m, p % coord0 % xyz, mesh_bin) + if (mesh_bin == NO_BIN_FOUND) then + found_bin = .false. + return + end if + bins(FILTER_MESH) = mesh_bin + else + bins(FILTER_MESH) = 1 + end if + + ! determine next universe bin + ! TODO: Account for multiple universes when performing this filter + if (t % n_filter_bins(FILTER_UNIVERSE) > 0) then + bins(FILTER_UNIVERSE) = get_next_bin(FILTER_UNIVERSE, & + p % coord % universe, index_tally) + if (bins(FILTER_UNIVERSE) == NO_BIN_FOUND) then + found_bin = .false. + return + end if + else + bins(FILTER_UNIVERSE) = 1 + end if + + ! determine next material bin + if (t % n_filter_bins(FILTER_MATERIAL) > 0) then + bins(FILTER_MATERIAL) = get_next_bin(FILTER_MATERIAL, & + p % material, index_tally) + if (bins(FILTER_MATERIAL) == NO_BIN_FOUND) then + found_bin = .false. + return + end if + else + bins(FILTER_MATERIAL) = 1 + end if + + ! determine next cell bin + ! TODO: Account for cells in multiple levels when performing this filter + if (t % n_filter_bins(FILTER_CELL) > 0) then + bins(FILTER_CELL) = get_next_bin(FILTER_CELL, & + p % coord % cell, index_tally) + if (bins(FILTER_CELL) == NO_BIN_FOUND) then + found_bin = .false. + return + end if + else + bins(FILTER_CELL) = 1 + end if + + ! determine next cellborn bin + if (t % n_filter_bins(FILTER_CELLBORN) > 0) then + bins(FILTER_CELLBORN) = get_next_bin(FILTER_CELLBORN, & + p % cell_born, index_tally) + if (bins(FILTER_CELLBORN) == NO_BIN_FOUND) then + found_bin = .false. + return + end if + else + bins(FILTER_CELLBORN) = 1 + end if + + ! determine next surface bin + if (t % n_filter_bins(FILTER_SURFACE) > 0) then + bins(FILTER_SURFACE) = get_next_bin(FILTER_SURFACE, & + p % surface, index_tally) + if (bins(FILTER_SURFACE) == NO_BIN_FOUND) then + found_bin = .false. + return + end if + else + bins(FILTER_SURFACE) = 1 + end if + + ! determine incoming energy bin + n = t % n_filter_bins(FILTER_ENERGYIN) + if (n > 0) then + ! check if energy of the particle is within energy bins + if (p % last_E < t % energy_in(1) .or. & + p % last_E > t % energy_in(n + 1)) then + found_bin = .false. + return + end if + + ! search to find incoming energy bin + bins(FILTER_ENERGYIN) = binary_search(t % energy_in, n + 1, p % last_E) + else + bins(FILTER_ENERGYIN) = 1 + end if + + ! determine outgoing energy bin + n = t % n_filter_bins(FILTER_ENERGYOUT) + if (n > 0) then + ! check if energy of the particle is within energy bins + if (p % E < t % energy_out(1) .or. p % E > t % energy_out(n + 1)) then + found_bin = .false. + return + end if + + ! search to find incoming energy bin + bins(FILTER_ENERGYOUT) = binary_search(t % energy_out, n + 1, p % E) + else + bins(FILTER_ENERGYOUT) = 1 + end if + + end subroutine get_scoring_bins !=============================================================================== ! SCORE_SURFACE_CURRENT tallies surface crossings in a mesh tally by manually