crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Removed matching_bins from TallyObject type and made it threadprivate. Still have race condition with Lists.

Browse source
This commit is contained in:
Paul Romano 2013-08-16 19:34:51 -05:00
parent 31db282802
commit be05b1ce55
6 changed files with 158 additions and 153 deletions
Download patch file Download diff file Expand all files Collapse all files

85
src/cmfd_data.F90
View file

@ -57,7 +57,8 @@ contains
FILTER_SURFACE, IN_RIGHT, OUT_RIGHT, IN_FRONT, &
OUT_FRONT, IN_TOP, OUT_TOP, CMFD_NOACCEL, ZERO, ONE
use error, only: fatal_error
use global, only: cmfd, message, n_cmfd_tallies, cmfd_tallies, meshes
use global, only: cmfd, message, n_cmfd_tallies, cmfd_tallies, meshes,&
matching_bins
use mesh, only: mesh_indices_to_bin
use mesh_header, only: StructuredMesh
use tally_header, only: TallyObject
@ -137,21 +138,21 @@ contains
TALLY: if (ital == 1) then
! reset all bins to 1
t % matching_bins = 1
matching_bins(1:t%n_filters) = 1
! set ijk as mesh indices
ijk = (/ i, j, k /)
! get bin number for mesh indices
t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m,ijk)
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m,ijk)
! apply energy in filter
if (i_filter_ein > 0) then
t % matching_bins(i_filter_ein) = ng - h + 1
matching_bins(i_filter_ein) = ng - h + 1
end if
! calculate score index from bins
score_index = sum((t % matching_bins - 1) * t%stride) + 1
score_index = sum((matching_bins(1:t%n_filters) - 1) * t%stride) + 1
! get flux
flux = t % results(1,score_index) % sum
@ -180,24 +181,24 @@ contains
INGROUP: do g = 1, ng
! reset all bins to 1
t % matching_bins = 1
matching_bins(1:t%n_filters) = 1
! set ijk as mesh indices
ijk = (/ i, j, k /)
! get bin number for mesh indices
t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m,ijk)
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m,ijk)
if (i_filter_ein > 0) then
! apply energy in filter
t % matching_bins(i_filter_ein) = ng - h + 1
matching_bins(i_filter_ein) = ng - h + 1
! set energy out bin
t % matching_bins(i_filter_eout) = ng - g + 1
matching_bins(i_filter_eout) = ng - g + 1
end if
! calculate score index from bins
score_index = sum((t % matching_bins - 1) * t%stride) + 1
score_index = sum((matching_bins(1:t%n_filters) - 1) * t%stride) + 1
! get scattering
cmfd % scattxs(h,g,i,j,k) = t % results(1,score_index) % sum /&
@ -216,69 +217,69 @@ contains
else if (ital == 3) then
! initialize and filter for energy
t % matching_bins = 1
matching_bins(1:t%n_filters) = 1
if (i_filter_ein > 0) then
t % matching_bins(i_filter_ein) = ng - h + 1
matching_bins(i_filter_ein) = ng - h + 1
end if
! left surface
t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i-1, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_RIGHT
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing
matching_bins(i_filter_surf) = IN_RIGHT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
cmfd % current(1,h,i,j,k) = t % results(1,score_index) % sum
t % matching_bins(i_filter_surf) = OUT_RIGHT
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming
matching_bins(i_filter_surf) = OUT_RIGHT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
cmfd % current(2,h,i,j,k) = t % results(1,score_index) % sum
! right surface
t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_RIGHT
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming
matching_bins(i_filter_surf) = IN_RIGHT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
cmfd % current(3,h,i,j,k) = t % results(1,score_index) % sum
t % matching_bins(i_filter_surf) = OUT_RIGHT
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing
matching_bins(i_filter_surf) = OUT_RIGHT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
cmfd % current(4,h,i,j,k) = t % results(1,score_index) % sum
! back surface
t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i, j-1, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_FRONT
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing
matching_bins(i_filter_surf) = IN_FRONT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
cmfd % current(5,h,i,j,k) = t % results(1,score_index) % sum
t % matching_bins(i_filter_surf) = OUT_FRONT
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming
matching_bins(i_filter_surf) = OUT_FRONT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
cmfd % current(6,h,i,j,k) = t % results(1,score_index) % sum
! front surface
t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_FRONT
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming
matching_bins(i_filter_surf) = IN_FRONT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
cmfd % current(7,h,i,j,k) = t % results(1,score_index) % sum
t % matching_bins(i_filter_surf) = OUT_FRONT
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing
matching_bins(i_filter_surf) = OUT_FRONT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
cmfd % current(8,h,i,j,k) = t % results(1,score_index) % sum
! bottom surface
t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i, j, k-1 /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_TOP
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing
matching_bins(i_filter_surf) = IN_TOP
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
cmfd % current(9,h,i,j,k) = t % results(1,score_index) % sum
t % matching_bins(i_filter_surf) = OUT_TOP
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming
matching_bins(i_filter_surf) = OUT_TOP
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
cmfd % current(10,h,i,j,k) = t % results(1,score_index) % sum
! top surface
t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_TOP
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming
matching_bins(i_filter_surf) = IN_TOP
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
cmfd % current(11,h,i,j,k) = t % results(1,score_index) % sum
t % matching_bins(i_filter_surf) = OUT_TOP
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing
matching_bins(i_filter_surf) = OUT_TOP
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
cmfd % current(12,h,i,j,k) = t % results(1,score_index) % sum
end if TALLY

3
src/global.F90
View file

@ -97,6 +97,7 @@ module global
type(StructuredMesh), allocatable, target :: meshes(:)
type(TallyObject), allocatable, target :: tallies(:)
integer, allocatable :: matching_bins(:)
! Pointers for different tallies
type(TallyObject), pointer :: user_tallies(:) => null()
@ -384,7 +385,7 @@ module global
logical :: output_tallies = .true.
!$omp threadprivate(micro_xs, material_xs, fission_bank, n_bank, message, &
!$omp& trace, thread_id, current_work)
!$omp& trace, thread_id, current_work, matching_bins)
contains

76
src/output.F90
View file

@ -1644,7 +1644,7 @@ contains
! to be used for a given tally.
! Initialize bins, filter level, and indentation
t % matching_bins = 0
matching_bins(1:t%n_filters) = 0
j = 1
indent = 0
@ -1654,16 +1654,16 @@ contains
if (t % n_filters == 0) exit find_bin
! Increment bin combination
t % matching_bins(j) = t % matching_bins(j) + 1
matching_bins(j) = matching_bins(j) + 1
! =================================================================
! REACHED END OF BINS FOR THIS FILTER, MOVE TO NEXT FILTER
if (t % matching_bins(j) > t % filters(j) % n_bins) then
if (matching_bins(j) > t % filters(j) % n_bins) then
! If this is the first filter, then exit
if (j == 1) exit print_bin
t % matching_bins(j) = 0
matching_bins(j) = 0
j = j - 1
indent = indent - 2
@ -1696,7 +1696,7 @@ contains
! bins below the lowest filter level will be zeros
if (t % n_filters > 0) then
filter_index = sum((max(t % matching_bins,1) - 1) * t % stride) + 1
filter_index = sum((max(matching_bins(1:t%n_filters),1) - 1) * t % stride) + 1
else
filter_index = 1
end if
@ -1805,7 +1805,7 @@ contains
m => meshes(t % filters(i_filter_mesh) % int_bins(1))
! initialize bins array
t % matching_bins = 1
matching_bins(1:t%n_filters) = 1
! determine how many energy in bins there are
i_filter_ein = t % find_filter(FILTER_ENERGYIN)
@ -1831,7 +1831,7 @@ contains
do l = 1, n
if (print_ebin) then
! Set incoming energy bin
t % matching_bins(i_filter_ein) = l
matching_bins(i_filter_ein) = l
! Write incoming energy bin
write(UNIT=UNIT_TALLY, FMT='(3X,A,1X,A)') &
@ -1839,102 +1839,102 @@ contains
end if
! Left Surface
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i-1, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_RIGHT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
matching_bins(i_filter_surf) = IN_RIGHT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Left", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
t % matching_bins(i_filter_surf) = OUT_RIGHT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
matching_bins(i_filter_surf) = OUT_RIGHT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Left", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
! Right Surface
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_RIGHT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
matching_bins(i_filter_surf) = IN_RIGHT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Right", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
t % matching_bins(i_filter_surf) = OUT_RIGHT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
matching_bins(i_filter_surf) = OUT_RIGHT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Right", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
! Back Surface
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j-1, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_FRONT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
matching_bins(i_filter_surf) = IN_FRONT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Back", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
t % matching_bins(i_filter_surf) = OUT_FRONT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
matching_bins(i_filter_surf) = OUT_FRONT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Back", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
! Front Surface
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_FRONT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
matching_bins(i_filter_surf) = IN_FRONT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Front", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
t % matching_bins(i_filter_surf) = OUT_FRONT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
matching_bins(i_filter_surf) = OUT_FRONT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Front", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
! Bottom Surface
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k-1 /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_TOP
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
matching_bins(i_filter_surf) = IN_TOP
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Bottom", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
t % matching_bins(i_filter_surf) = OUT_TOP
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
matching_bins(i_filter_surf) = OUT_TOP
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Bottom", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
! Top Surface
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_TOP
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
matching_bins(i_filter_surf) = IN_TOP
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Top", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
t % matching_bins(i_filter_surf) = OUT_TOP
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
matching_bins(i_filter_surf) = OUT_TOP
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Top", &
to_str(t % results(1,filter_index) % sum), &
@ -1965,7 +1965,7 @@ contains
real(8) :: E1 ! upper bound for energy bin
type(StructuredMesh), pointer :: m => null()
bin = t % matching_bins(i_filter)
bin = matching_bins(i_filter)
select case(t % filters(i_filter) % type)
case (FILTER_UNIVERSE)

130
src/tally.F90
View file

@ -85,7 +85,7 @@ contains
! be accumulating the tally values
! Determine scoring index for this filter combination
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
! Check for nuclide bins
k = 0
@ -388,7 +388,7 @@ contains
! save original outgoing energy bin and score index
i = t % find_filter(FILTER_ENERGYOUT)
bin_energyout = t % matching_bins(i)
bin_energyout = matching_bins(i)
! Get number of energies on filter
n = size(t % filters(i) % real_bins)
@ -411,10 +411,10 @@ contains
E_out > t % filters(i) % real_bins(n)) cycle
! change outgoing energy bin
t % matching_bins(i) = binary_search(t % filters(i) % real_bins, n, E_out)
matching_bins(i) = binary_search(t % filters(i) % real_bins, n, E_out)
! determine scoring index
i_filter = sum((t % matching_bins - 1) * t % stride) + 1
i_filter = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
! Add score to tally
!$omp critical
@ -424,7 +424,7 @@ contains
end do
! reset outgoing energy bin and score index
t % matching_bins(i) = bin_energyout
matching_bins(i) = bin_energyout
end subroutine score_fission_eout
@ -495,7 +495,7 @@ contains
! be accumulating the tally values
! Determine scoring index for this filter combination
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
if (t % all_nuclides) then
call score_all_nuclides(p, i_tally, flux, filter_index)
@ -1005,7 +1005,7 @@ contains
!$omp threadprivate(t, m, mat, coord)
t => tallies(i_tally)
t % matching_bins = 1
matching_bins(1:t%n_filters) = 1
! ==========================================================================
! CHECK IF THIS TRACK INTERSECTS THE MESH
@ -1045,11 +1045,11 @@ contains
case (FILTER_UNIVERSE)
! determine next universe bin
! TODO: Account for multiple universes when performing this filter
t % matching_bins(i) = get_next_bin(FILTER_UNIVERSE, &
matching_bins(i) = get_next_bin(FILTER_UNIVERSE, &
p % coord % universe, i_tally)
case (FILTER_MATERIAL)
t % matching_bins(i) = get_next_bin(FILTER_MATERIAL, &
matching_bins(i) = get_next_bin(FILTER_MATERIAL, &
p % material, i_tally)
case (FILTER_CELL)
@ -1057,21 +1057,21 @@ contains
coord => p % coord0
do while(associated(coord))
position(FILTER_CELL) = 0
t % matching_bins(i) = get_next_bin(FILTER_CELL, &
matching_bins(i) = get_next_bin(FILTER_CELL, &
coord % cell, i_tally)
if (t % matching_bins(i) /= NO_BIN_FOUND) exit
if (matching_bins(i) /= NO_BIN_FOUND) exit
coord => coord % next
end do
nullify(coord)
case (FILTER_CELLBORN)
! determine next cellborn bin
t % matching_bins(i) = get_next_bin(FILTER_CELLBORN, &
matching_bins(i) = get_next_bin(FILTER_CELLBORN, &
p % cell_born, i_tally)
case (FILTER_SURFACE)
! determine next surface bin
t % matching_bins(i) = get_next_bin(FILTER_SURFACE, &
matching_bins(i) = get_next_bin(FILTER_SURFACE, &
p % surface, i_tally)
case (FILTER_ENERGYIN)
@ -1081,17 +1081,17 @@ contains
! check if energy of the particle is within energy bins
if (p % E < t % filters(i) % real_bins(1) .or. &
p % E > t % filters(i) % real_bins(k + 1)) then
t % matching_bins(i) = NO_BIN_FOUND
matching_bins(i) = NO_BIN_FOUND
else
! search to find incoming energy bin
t % matching_bins(i) = binary_search(t % filters(i) % real_bins, &
matching_bins(i) = binary_search(t % filters(i) % real_bins, &
k + 1, p % E)
end if
end select
! Check if no matching bin was found
if (t % matching_bins(i) == NO_BIN_FOUND) return
if (matching_bins(i) == NO_BIN_FOUND) return
end do FILTER_LOOP
@ -1161,10 +1161,10 @@ contains
flux = p % wgt * distance
! Determine mesh bin
t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, ijk_cross)
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, ijk_cross)
! Determining scoring index
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
if (t % all_nuclides) then
! Score reaction rates for each nuclide in material
@ -1301,7 +1301,7 @@ contains
found_bin = .true.
t => tallies(i_tally)
t % matching_bins = 1
matching_bins(1:t%n_filters) = 1
FILTER_LOOP: do i = 1, t % n_filters
@ -1311,16 +1311,16 @@ contains
m => meshes(t % filters(i) % int_bins(1))
! Determine if we're in the mesh first
call get_mesh_bin(m, p % coord0 % xyz, t % matching_bins(i))
call get_mesh_bin(m, p % coord0 % xyz, matching_bins(i))
case (FILTER_UNIVERSE)
! determine next universe bin
! TODO: Account for multiple universes when performing this filter
t % matching_bins(i) = get_next_bin(FILTER_UNIVERSE, &
matching_bins(i) = get_next_bin(FILTER_UNIVERSE, &
p % coord % universe, i_tally)
case (FILTER_MATERIAL)
t % matching_bins(i) = get_next_bin(FILTER_MATERIAL, &
matching_bins(i) = get_next_bin(FILTER_MATERIAL, &
p % material, i_tally)
case (FILTER_CELL)
@ -1328,21 +1328,21 @@ contains
coord => p % coord0
do while(associated(coord))
position(FILTER_CELL) = 0
t % matching_bins(i) = get_next_bin(FILTER_CELL, &
matching_bins(i) = get_next_bin(FILTER_CELL, &
coord % cell, i_tally)
if (t % matching_bins(i) /= NO_BIN_FOUND) exit
if (matching_bins(i) /= NO_BIN_FOUND) exit
coord => coord % next
end do
nullify(coord)
case (FILTER_CELLBORN)
! determine next cellborn bin
t % matching_bins(i) = get_next_bin(FILTER_CELLBORN, &
matching_bins(i) = get_next_bin(FILTER_CELLBORN, &
p % cell_born, i_tally)
case (FILTER_SURFACE)
! determine next surface bin
t % matching_bins(i) = get_next_bin(FILTER_SURFACE, &
matching_bins(i) = get_next_bin(FILTER_SURFACE, &
p % surface, i_tally)
case (FILTER_ENERGYIN)
@ -1359,10 +1359,10 @@ contains
! check if energy of the particle is within energy bins
if (E < t % filters(i) % real_bins(1) .or. &
E > t % filters(i) % real_bins(n + 1)) then
t % matching_bins(i) = NO_BIN_FOUND
matching_bins(i) = NO_BIN_FOUND
else
! search to find incoming energy bin
t % matching_bins(i) = binary_search(t % filters(i) % real_bins, &
matching_bins(i) = binary_search(t % filters(i) % real_bins, &
n + 1, E)
end if
@ -1373,17 +1373,17 @@ contains
! check if energy of the particle is within energy bins
if (p % E < t % filters(i) % real_bins(1) .or. &
p % E > t % filters(i) % real_bins(n + 1)) then
t % matching_bins(i) = NO_BIN_FOUND
matching_bins(i) = NO_BIN_FOUND
else
! search to find incoming energy bin
t % matching_bins(i) = binary_search(t % filters(i) % real_bins, &
matching_bins(i) = binary_search(t % filters(i) % real_bins, &
n + 1, p % E)
end if
end select
! If the current filter didn't match, exit this subroutine
if (t % matching_bins(i) == NO_BIN_FOUND) then
if (matching_bins(i) == NO_BIN_FOUND) then
found_bin = .false.
return
end if
@ -1475,7 +1475,7 @@ contains
end if
! search to find incoming energy bin
t % matching_bins(j) = binary_search(t % filters(j) % real_bins, &
matching_bins(j) = binary_search(t % filters(j) % real_bins, &
n + 1, p % E)
end if
@ -1492,10 +1492,10 @@ contains
do j = ijk0(3), ijk1(3) - 1
ijk0(3) = j
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = OUT_TOP
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_surf) = OUT_TOP
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
!$omp critical
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
@ -1506,10 +1506,10 @@ contains
do j = ijk0(3) - 1, ijk1(3), -1
ijk0(3) = j
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = IN_TOP
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_surf) = IN_TOP
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
!$omp critical
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
@ -1524,10 +1524,10 @@ contains
do j = ijk0(2), ijk1(2) - 1
ijk0(2) = j
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = OUT_FRONT
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_surf) = OUT_FRONT
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
!$omp critical
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
@ -1538,10 +1538,10 @@ contains
do j = ijk0(2) - 1, ijk1(2), -1
ijk0(2) = j
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = IN_FRONT
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_surf) = IN_FRONT
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
!$omp critical
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
@ -1556,10 +1556,10 @@ contains
do j = ijk0(1), ijk1(1) - 1
ijk0(1) = j
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = OUT_RIGHT
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_surf) = OUT_RIGHT
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
!$omp critical
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
@ -1570,10 +1570,10 @@ contains
do j = ijk0(1) - 1, ijk1(1), -1
ijk0(1) = j
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = IN_RIGHT
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_surf) = IN_RIGHT
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
!$omp critical
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
@ -1598,7 +1598,7 @@ contains
do k = 1, n_cross
! Reset scoring bin index
t % matching_bins(i_filter_surf) = 0
matching_bins(i_filter_surf) = 0
! Calculate distance to each bounding surface. We need to treat
! special case where the cosine of the angle is zero since this would
@ -1625,8 +1625,8 @@ contains
! Crossing into right mesh cell -- this is treated as outgoing
! current from (i,j,k)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = OUT_RIGHT
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_surf) = OUT_RIGHT
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
end if
ijk0(1) = ijk0(1) + 1
@ -1637,8 +1637,8 @@ contains
ijk0(1) = ijk0(1) - 1
xyz_cross(1) = xyz_cross(1) - m % width(1)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = IN_RIGHT
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_surf) = IN_RIGHT
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
end if
end if
@ -1647,8 +1647,8 @@ contains
! Crossing into front mesh cell -- this is treated as outgoing
! current in (i,j,k)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = OUT_FRONT
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_surf) = OUT_FRONT
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
end if
ijk0(2) = ijk0(2) + 1
@ -1659,8 +1659,8 @@ contains
ijk0(2) = ijk0(2) - 1
xyz_cross(2) = xyz_cross(2) - m % width(2)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = IN_FRONT
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_surf) = IN_FRONT
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
end if
end if
@ -1669,8 +1669,8 @@ contains
! Crossing into top mesh cell -- this is treated as outgoing
! current in (i,j,k)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = OUT_TOP
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_surf) = OUT_TOP
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
end if
ijk0(3) = ijk0(3) + 1
@ -1681,16 +1681,16 @@ contains
ijk0(3) = ijk0(3) - 1
xyz_cross(3) = xyz_cross(3) - m % width(3)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = IN_TOP
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_surf) = IN_TOP
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
end if
end if
end if
! Determine scoring index
if (t % matching_bins(i_filter_surf) > 0) then
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
if (matching_bins(i_filter_surf) > 0) then
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
! Check for errors
if (filter_index <= 0 .or. filter_index > &

3
src/tally_header.F90
View file

@ -86,7 +86,6 @@ module tally_header
! mapped onto one dimension in the results array, the stride attribute gives
! the stride for a given filter type within the results array
integer, allocatable :: matching_bins(:)
integer, allocatable :: stride(:)
! This array provides a way to lookup what index in the filters array a
@ -169,8 +168,6 @@ module tally_header
deallocate(this % filters)
end if
if (allocated(this % matching_bins)) &
deallocate(this % matching_bins)
if (allocated(this % stride)) &
deallocate(this % stride)

14
src/tally_initialize.F90
View file

@ -31,9 +31,10 @@ contains
subroutine setup_tally_arrays()
integer :: i ! loop index for tallies
integer :: j ! loop index for filters
integer :: n ! temporary stride
integer :: i ! loop index for tallies
integer :: j ! loop index for filters
integer :: n ! temporary stride
integer :: max_n_filters = 0 ! maximum number of filters
type(TallyObject), pointer :: t => null()
TALLY_LOOP: do i = 1, n_tallies
@ -42,7 +43,7 @@ contains
! Allocate stride and matching_bins arrays
allocate(t % stride(t % n_filters))
allocate(t % matching_bins(t % n_filters))
max_n_filters = max(max_n_filters, t % n_filters)
! The filters are traversed in opposite order so that the last filter has
! the shortest stride in memory and the first filter has the largest
@ -63,6 +64,11 @@ contains
end do TALLY_LOOP
! Allocate array for matching filter bins
!$omp parallel
allocate(matching_bins(max_n_filters))
!$omp end parallel
end subroutine setup_tally_arrays
!===============================================================================