diff --git a/src/constants.f90 b/src/constants.f90 index dd2cdfbd00..9d7f46c410 100644 --- a/src/constants.f90 +++ b/src/constants.f90 @@ -160,6 +160,9 @@ module constants MACRO_FISSION = -5, & ! total fission rate MACRO_NU_FISSION = -6 ! total neutron production rate + ! Tally map bin finding + integer, parameter :: NO_BIN_FOUND = -1 + ! Tally map types integer, parameter :: TALLY_MAP_TYPES = 6 integer, parameter :: & diff --git a/src/physics.f90 b/src/physics.f90 index 0ffe000281..9f50ab1d79 100644 --- a/src/physics.f90 +++ b/src/physics.f90 @@ -278,7 +278,6 @@ contains real(8) :: total ! total macroscopic xs for material real(8) :: prob ! cumulative probability real(8) :: cutoff ! random number - real(8) :: flux ! collision estimator of flux real(8) :: atom_density ! atom density of nuclide in atom/b-cm character(MAX_LINE_LEN) :: msg ! output/error message type(Material), pointer :: mat @@ -287,6 +286,12 @@ contains mat => materials(p % material) + ! Score collision estimator tallies for any macro tallies -- we can do this + ! before sampling the nuclide since + if (tallies_on) then + ! call score_tally(p) + end if + ! sample nuclide i = 0 total = material_xs % total @@ -353,12 +358,6 @@ contains call warning(msg) end select - ! Add collision estimator tallies - if (tallies_on) then - flux = p%wgt / SIGMA - call score_tally(p, flux) - end if - ! check for very low energy if (p % E < 1.0e-100_8) then p % alive = .false. diff --git a/src/tally.f90 b/src/tally.f90 index a4cdf7a920..e2bfc8f555 100644 --- a/src/tally.f90 +++ b/src/tally.f90 @@ -92,8 +92,9 @@ contains integer :: i integer :: j integer :: index + integer :: n + integer :: score_bins type(TallyObject), pointer :: t => null() - type(Cell), pointer :: c => null() ! allocate tally map array allocate(tally_maps(TALLY_MAP_TYPES)) @@ -109,46 +110,79 @@ contains do i = 1, n_tallies t => tallies(i) + ! initialize number of scoring bins + score_bins = 1 + ! Add map elements for cell bins if (associated(t % cell_bins)) then - do j = 1, size(t % cell_bins) + n = size(t % cell_bins) + do j = 1, n index = t % cell_bins(j) % scalar call add_map_element(tally_maps(MAP_CELL) % items(index), i, j) end do + score_bins = score_bins * n end if ! Add map elements for surface bins if (associated(t % surface_bins)) then - do j = 1, size(t % surface_bins) + n = size(t % surface_bins) + do j = 1, n index = t % surface_bins(j) % scalar call add_map_element(tally_maps(MAP_SURFACE) % items(index), i, j) end do + score_bins = score_bins * n end if ! Add map elements for universe bins if (associated(t % universe_bins)) then - do j = 1, size(t % universe_bins) + n = size(t % universe_bins) + do j = 1, n index = t % universe_bins(j) % scalar call add_map_element(tally_maps(MAP_UNIVERSE) % items(index), i, j) end do + score_bins = score_bins * n end if ! Add map elements for material bins if (associated(t % material_bins)) then - do j = 1, size(t % material_bins) + n = size(t % material_bins) + do j = 1, n index = t % material_bins(j) % scalar call add_map_element(tally_maps(MAP_MATERIAL) % items(index), i, j) end do + score_bins = score_bins * n end if ! Add map elements for bornin bins if (associated(t % bornin_bins)) then + n = size(t % bornin_bins) do j = 1, size(t % bornin_bins) index = t % bornin_bins(j) % scalar call add_map_element(tally_maps(MAP_BORNIN) % items(index), i, j) end do + score_bins = score_bins * n end if + ! TODO: Determine size of mesh to increase number of scoring bins + + ! determine if there are subdivisions for incoming or outgoing energy to + ! adjust the number of scoring bins appropriately + if (allocated(t % energy_in)) then + score_bins = score_bins * (size(t % energy_in) - 1) + end if + + if (allocated(t % energy_out)) then + score_bins = score_bins * (size(t % energy_out) - 1) + end if + + ! Finally add scoring bins for the macro tallies + if (associated(t % macro_bins)) then + score_bins = score_bins * size(t % macro_bins) + end if + + ! Allocate scores for tally + allocate(t % score(score_bins)) + end do end subroutine create_tally_map @@ -193,189 +227,154 @@ contains ! SCORE_TALLY !=============================================================================== - subroutine score_tally(p, flux) + subroutine score_tally(p) type(Particle), pointer :: p ! particle - real(8), intent(in) :: flux ! estimator of flux -!!$ integer :: i ! index for tallies in cell -!!$ integer :: j ! index for reactions in tally -!!$ integer :: e_bin ! energy bin -!!$ integer :: c_bin ! cell bin -!!$ integer :: r_bin ! reaction bin -!!$ integer :: n_energy ! number of energy bins -!!$ integer :: n_reaction ! number of reactions -!!$ integer :: type ! cell or reaction type -!!$ integer :: MT ! reaction MT value -!!$ real(8) :: E ! energy of particle -!!$ real(8) :: val ! value to score -!!$ real(8) :: Sigma ! macroscopic cross section of reaction -!!$ character(MAX_LINE_LEN) :: msg ! output/error message -!!$ type(Cell), pointer :: c => null() -!!$ type(Tally), pointer :: t => null() -!!$ type(Material), pointer :: mat => null() -!!$ -!!$ ! ========================================================================== -!!$ ! HANDLE LOCAL TALLIES -!!$ -!!$ c => cells(p % cell) -!!$ if (.not. allocated(c % tallies)) return -!!$ -!!$ ! if so, loop over each tally -!!$ do i = 1, size(c % tallies) -!!$ t => tallies(c % tallies(i)) -!!$ -!!$ ! ======================================================================= -!!$ ! DETERMINE ENERGY BIN -!!$ if (allocated(t % energies)) then -!!$ E = p % E -!!$ n_energy = size(t % energies) -!!$ if (E < t % energies(1) .or. E > t % energies(n_energy)) then -!!$ ! Energy outside of specified range -!!$ cycle -!!$ else -!!$ e_bin = binary_search(t % energies, n_energy, E) -!!$ end if -!!$ end if -!!$ -!!$ ! ======================================================================= -!!$ ! DETERMINE CELL BIN -!!$ type = t % cell_type -!!$ if (type == TALLY_SUM) then -!!$ ! Sum tallies from separate cells into one bin -!!$ c_bin = 1 -!!$ elseif (type == TALLY_BINS) then -!!$ ! Need to determine cell bin -!!$ do j = 1, size(t % cells) -!!$ if (p % cell == t % cells(j)) then -!!$ c_bin = j -!!$ exit -!!$ end if -!!$ end do -!!$ else -!!$ msg = "Invalid type for cell bins in tally " // int_to_str(t % uid) -!!$ call fatal_error(msg) -!!$ end if -!!$ -!!$ ! ======================================================================= -!!$ ! DETERMINE REACTION BIN AND ADD VALUE TO SCORE -!!$ type = t % reaction_type -!!$ if (type == TALLY_FLUX) then -!!$ ! Tally flux only -!!$ val = flux -!!$ call add_to_score(t % score(r_bin, c_bin, e_bin), val) -!!$ -!!$ elseif (type == TALLY_ALL) then -!!$ ! Tally total reaction rate -!!$ val = ONE -!!$ r_bin = 1 -!!$ call add_to_score(t % score(r_bin, c_bin, e_bin), val) -!!$ -!!$ elseif (type == TALLY_BINS) then -!!$ ! Individually bin reactions -!!$ n_reaction = t % reaction_type -!!$ -!!$ r_bin = 0 -!!$ do j = 1, n_reaction -!!$ MT = t % reactions(j) -!!$ mat => materials(p % material) -!!$ Sigma = get_macro_xs(p, mat, MT) -!!$ val = Sigma * flux -!!$ r_bin = r_bin + 1 -!!$ call add_to_score(t % score(r_bin, c_bin, e_bin), & -!!$ & val) -!!$ end do -!!$ elseif (type == TALLY_SUM) then -!!$ ! Tally reactions in one bin -!!$ n_reaction = t % reaction_type -!!$ -!!$ r_bin = 1 -!!$ do j = 1, n_reaction -!!$ MT = t % reactions(j) -!!$ mat => materials(p % material) -!!$ Sigma = get_macro_xs(p, mat, MT) -!!$ val = Sigma * flux -!!$ call add_to_score(t % score(r_bin, c_bin, e_bin), & -!!$ & val) -!!$ end do -!!$ end if -!!$ end do -!!$ -!!$ ! ========================================================================== -!!$ ! HANDLE GLOBAL TALLIES -!!$ -!!$ do i = 1, n_tallies_global -!!$ t => tallies_global(i) -!!$ -!!$ ! ======================================================================= -!!$ ! DETERMINE ENERGY BIN -!!$ if (allocated(t % energies)) then -!!$ E = p % E -!!$ n_energy = size(t % energies) -!!$ if (E < t % energies(1) .or. E > t % energies(n_energy)) then -!!$ ! Energy outside of specified range -!!$ cycle -!!$ else -!!$ e_bin = binary_search(t % energies, n_energy, E) -!!$ end if -!!$ end if -!!$ -!!$ ! Since it's a global tally, the cell bin is unity -!!$ c_bin = 1 -!!$ -!!$ ! ======================================================================= -!!$ ! DETERMINE REACTION BIN AND ADD VALUE TO SCORE -!!$ type = t % reaction_type -!!$ if (type == TALLY_FLUX) then -!!$ ! Tally flux only -!!$ val = flux -!!$ call add_to_score(t % score(r_bin, c_bin, e_bin), val) -!!$ -!!$ elseif (type == TALLY_ALL) then -!!$ ! Tally total reaction rate -!!$ val = ONE -!!$ r_bin = 1 -!!$ call add_to_score(t % score(r_bin, c_bin, e_bin), val) -!!$ -!!$ elseif (type == TALLY_BINS) then -!!$ ! Individually bin reactions -!!$ n_reaction = t % reaction_type -!!$ -!!$ r_bin = 0 -!!$ do j = 1, n_reaction -!!$ MT = t % reactions(j) -!!$ mat => materials(p % material) -!!$ Sigma = get_macro_xs(p, mat, MT) -!!$ val = Sigma * flux -!!$ r_bin = r_bin + 1 -!!$ call add_to_score(t % score(r_bin, c_bin, e_bin), & -!!$ & val) -!!$ end do -!!$ elseif (type == TALLY_SUM) then -!!$ ! Tally reactions in one bin -!!$ n_reaction = t % reaction_type -!!$ -!!$ r_bin = 1 -!!$ do j = 1, n_reaction -!!$ MT = t % reactions(j) -!!$ mat => materials(p % material) -!!$ Sigma = get_macro_xs(p, mat, MT) -!!$ val = Sigma * flux -!!$ call add_to_score(t % score(r_bin, c_bin, e_bin), & -!!$ & val) -!!$ end do -!!$ end if -!!$ -!!$ end do -!!$ -!!$ ! ========================================================================== -!!$ ! TODO: Add lattice tallies -!!$ -!!$ ! ========================================================================== -!!$ ! TODO: Add mesh tallies + integer :: i + integer :: n + integer :: cell_bin + integer :: surface_bin + integer :: universe_bin + integer :: material_bin + integer :: bornin_bin + integer :: energyin_bin + integer :: energyout_bin + type(TallyObject), pointer :: t + + ! 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) + + ! determine next cell bin + if (associated(t % cell_bins)) then + cell_bin = get_next_bin(MAP_CELL, p % cell, i) + if (cell_bin == NO_BIN_FOUND) cycle + else + cell_bin = 1 + end if + + ! determine next surface bin + if (associated(t % surface_bins)) then + surface_bin = get_next_bin(MAP_SURFACE, p % surface, i) + if (surface_bin == NO_BIN_FOUND) cycle + else + surface_bin = 1 + end if + + ! determine next universe bin + if (associated(t % universe_bins)) then + universe_bin = get_next_bin(MAP_UNIVERSE, p % universe, i) + if (universe_bin == NO_BIN_FOUND) cycle + else + universe_bin = 1 + end if + + ! determine next material bin + if (associated(t % material_bins)) then + material_bin = get_next_bin(MAP_MATERIAL, p % material, i) + if (material_bin == NO_BIN_FOUND) cycle + else + material_bin = 1 + end if + + ! determine next bornin bin + if (associated(t % bornin_bins)) then + ! bornin_bin = get_next_bin(MAP_BORNIN, p % bornin, i) + if (bornin_bin == NO_BIN_FOUND) cycle + else + bornin_bin = 1 + end if + + ! determine incoming energy bin + if (allocated(t % energy_in)) then + ! check if energy of the particle is within energy bins + n = size(t % energy_in) + if (p % E < t % energy_in(1) .or. p % E > t % energy_in(n)) cycle + + ! search to find incoming energy bin + energyin_bin = binary_search(t % energy_in, n, p % E) + else + energyin_bin = 1 + end if + + ! determine outgoing energy bin + if (allocated(t % energy_out)) then + ! check if energy of the particle is within energy bins + n = size(t % energy_out) + if (p % E < t % energy_out(1) .or. p % E > t % energy_out(n)) cycle + + ! search to find incoming energy bin + energyout_bin = binary_search(t % energy_out, n, p % E) + else + energyout_bin = 1 + end if + + end do end subroutine score_tally +!=============================================================================== +! GET_NEXT_BIN +!=============================================================================== + + function get_next_bin(i_map, i_cell, i_tally) result(bin) + + integer, intent(in) :: i_map + integer, intent(in) :: i_cell + integer, intent(in) :: i_tally + integer :: bin + + integer :: index_tally + integer :: index_bin + integer, allocatable, save :: position(:) + logical, save :: first_entry = .true. + integer :: n + + ! initialize position to zero + if (first_entry) then + allocate(position(TALLY_MAP_TYPES)) + position = 0 + first_entry = .false. + end if + + n = size(tally_maps(i_map) % items(i_cell) % elements) + + do + ! Increment position in elements + position(i_map) = position(i_map) + 1 + + ! If we've reached the end of the array, there is no more bin to score to + if (position(i_map) > n) then + position(i_map) = 0 + bin = NO_BIN_FOUND + return + end if + + index_tally = tally_maps(i_map) % items(i_cell) % & + elements(position(i_map)) % index_tally + index_bin = tally_maps(i_map) % items(i_cell) % & + elements(position(i_map)) % index_bin + + if (index_tally > i_tally) then + ! Since the index being checked against is greater than the index we + ! need (and the tally indices were added to elements sequentially), we + ! know that no more bins will be scoring bins for this tally + position(i_map) = 0 + bin = NO_BIN_FOUND + return + elseif (index_tally == i_tally) then + ! Found a match + bin = index_bin + return + end if + + end do + + end function get_next_bin + !=============================================================================== ! ADD_TO_SCORE !===============================================================================