Greatly modified score_tally subroutine and added get_next_bin subroutine.

This commit is contained in:
Paul Romano 2011-09-20 14:04:48 -04:00
parent eb97f30e4b
commit 51666b4d52
3 changed files with 190 additions and 189 deletions

View file

@ -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 :: &

View file

@ -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.

View file

@ -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
!===============================================================================