crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 13:45:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Restructed tally score loop; replaced matching_bins with and filter_weights with new TallyFilterMatch array

Browse source
This commit is contained in:
amandalund 2017-04-05 20:39:12 -05:00
parent fe1f6062a4
commit 1f750bface
9 changed files with 642 additions and 342 deletions
Download patch file Download diff file Expand all files Collapse all files

179
src/cmfd_data.F90
View file

@ -57,7 +57,7 @@ contains
ZERO, ONE, TINY_BIT
use error, only: fatal_error
use global, only: cmfd, n_cmfd_tallies, cmfd_tallies, meshes, &
filters, matching_bins
filters, filter_matches
use mesh, only: mesh_indices_to_bin
use mesh_header, only: RegularMesh
use string, only: to_str
@ -72,6 +72,7 @@ contains
integer :: k ! iteration counter for z
integer :: g ! iteration counter for g
integer :: h ! iteration counter for outgoing groups
integer :: l ! iteration counter for tally filters
integer :: ital ! tally object index
integer :: ijk(3) ! indices for mesh cell
integer :: score_index ! index to pull from tally object
@ -81,6 +82,7 @@ contains
integer :: i_filter_ein ! index for incoming energy filter
integer :: i_filter_eout ! index for outgoing energy filter
integer :: i_filter_surf ! index for surface filter
logical :: energy_filters! energy filters present
real(8) :: flux ! temp variable for flux
type(TallyObject), pointer :: t ! pointer for tally object
type(RegularMesh), pointer :: m ! pointer for mesh object
@ -123,10 +125,14 @@ contains
end select
m => meshes(i_mesh)
i_filter_mesh = t % find_filter(FILTER_MESH)
i_filter_ein = t % find_filter(FILTER_ENERGYIN)
i_filter_eout = t % find_filter(FILTER_ENERGYOUT)
i_filter_surf = t % find_filter(FILTER_SURFACE)
! Check for energy filters
energy_filters = (t % find_filter(FILTER_ENERGYIN) > 0)
i_filter_mesh = t % filter(t % find_filter(FILTER_MESH))
if (energy_filters) then
i_filter_ein = t % filter(t % find_filter(FILTER_ENERGYIN))
i_filter_eout = t % filter(t % find_filter(FILTER_ENERGYOUT))
end if
! Begin loop around space
ZLOOP: do k = 1,nz
@ -149,22 +155,29 @@ contains
TALLY: if (ital == 1) then
! Reset all bins to 1
matching_bins(1:size(t % filter)) = 1
do l = 1, size(t % filter)
call filter_matches(t % filter(l)) % bins % clear()
call filter_matches(t % filter(l)) % bins % push_back(1)
end do
! Set ijk as mesh indices
ijk = (/ i, j, k /)
! Get bin number for mesh indices
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m,ijk)
filter_matches(i_filter_mesh) % bins % data(1) = &
mesh_indices_to_bin(m,ijk)
! Apply energy in filter
if (i_filter_ein > 0) then
matching_bins(i_filter_ein) = ng - h + 1
if (energy_filters) then
filter_matches(i_filter_ein) % bins % data(1) = ng - h + 1
end if
! Calculate score index from bins
score_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t%stride) + 1
score_index = 1
do l = 1, size(t % filter)
score_index = score_index + (filter_matches(t % filter(l)) &
% bins % data(1) - 1) * t % stride(l)
end do
! Get flux
flux = t % results(RESULT_SUM,1,score_index)
@ -193,25 +206,32 @@ contains
INGROUP: do g = 1, ng
! Reset all bins to 1
matching_bins(1:size(t % filter)) = 1
do l = 1, size(t % filter)
call filter_matches(t % filter(l)) % bins % clear()
call filter_matches(t % filter(l)) % bins % push_back(1)
end do
! Set ijk as mesh indices
ijk = (/ i, j, k /)
! Get bin number for mesh indices
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m,ijk)
filter_matches(i_filter_mesh) % bins % data(1) = &
mesh_indices_to_bin(m,ijk)
if (i_filter_ein > 0) then
if (energy_filters) then
! Apply energy in filter
matching_bins(i_filter_ein) = ng - h + 1
filter_matches(i_filter_ein) % bins % data(1) = ng - h + 1
! Set energy out bin
matching_bins(i_filter_eout) = ng - g + 1
filter_matches(i_filter_eout) % bins % data(1) = ng - g + 1
end if
! Calculate score index from bins
score_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t%stride) + 1
score_index = 1
do l = 1, size(t % filter)
score_index = score_index + (filter_matches(t % filter(l)) &
% bins % data(1) - 1) * t % stride(l)
end do
! Get scattering
cmfd % scattxs(h,g,i,j,k) = t % results(RESULT_SUM,1,score_index) /&
@ -231,80 +251,121 @@ contains
else if (ital == 3) then
i_filter_surf = t % filter(t % find_filter(FILTER_SURFACE))
! Initialize and filter for energy
matching_bins(1:size(t % filter)) = 1
if (i_filter_ein > 0) then
matching_bins(i_filter_ein) = ng - h + 1
do l = 1, size(t % filter)
call filter_matches(t % filter(l)) % bins % clear()
call filter_matches(t % filter(l)) % bins % push_back(1)
end do
if (energy_filters) then
filter_matches(i_filter_ein) % bins % data(1) = ng - h + 1
end if
! Get the bin for this mesh cell
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i, j, k /))
filter_matches(i_filter_mesh) % bins % data(1) = &
mesh_indices_to_bin(m, (/ i, j, k /))
! Left surface
matching_bins(i_filter_surf) = OUT_LEFT
score_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = OUT_LEFT
score_index = 1
do l = 1, size(t % filter)
score_index = score_index + (filter_matches(t % filter(l)) &
% bins % data(1) - 1) * t % stride(l)
end do
cmfd % current(1,h,i,j,k) = t % results(RESULT_SUM,1,score_index)
matching_bins(i_filter_surf) = IN_LEFT
score_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = IN_LEFT
score_index = 1
do l = 1, size(t % filter)
score_index = score_index + (filter_matches(t % filter(l)) &
% bins % data(1) - 1) * t % stride(l)
end do
cmfd % current(2,h,i,j,k) = t % results(RESULT_SUM,1,score_index)
! Right surface
matching_bins(i_filter_surf) = IN_RIGHT
score_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = IN_RIGHT
score_index = 1
do l = 1, size(t % filter)
score_index = score_index + (filter_matches(t % filter(l)) &
% bins % data(1) - 1) * t % stride(l)
end do
cmfd % current(3,h,i,j,k) = t % results(RESULT_SUM,1,score_index)
matching_bins(i_filter_surf) = OUT_RIGHT
score_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = OUT_RIGHT
score_index = 1
do l = 1, size(t % filter)
score_index = score_index + (filter_matches(t % filter(l)) &
% bins % data(1) - 1) * t % stride(l)
end do
cmfd % current(4,h,i,j,k) = t % results(RESULT_SUM,1,score_index)
! Back surface
matching_bins(i_filter_surf) = OUT_BACK
score_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = OUT_BACK
score_index = 1
do l = 1, size(t % filter)
score_index = score_index + (filter_matches(t % filter(l)) &
% bins % data(1) - 1) * t % stride(l)
end do
cmfd % current(5,h,i,j,k) = t % results(RESULT_SUM,1,score_index)
matching_bins(i_filter_surf) = IN_BACK
score_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = IN_BACK
score_index = 1
do l = 1, size(t % filter)
score_index = score_index + (filter_matches(t % filter(l)) &
% bins % data(1) - 1) * t % stride(l)
end do
cmfd % current(6,h,i,j,k) = t % results(RESULT_SUM,1,score_index)
! Front surface
matching_bins(i_filter_surf) = IN_FRONT
score_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = IN_FRONT
score_index = 1
do l = 1, size(t % filter)
score_index = score_index + (filter_matches(t % filter(l)) &
% bins % data(1) - 1) * t % stride(l)
end do
cmfd % current(7,h,i,j,k) = t % results(RESULT_SUM,1,score_index)
matching_bins(i_filter_surf) = OUT_FRONT
score_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = OUT_FRONT
score_index = 1
do l = 1, size(t % filter)
score_index = score_index + (filter_matches(t % filter(l)) &
% bins % data(1) - 1) * t % stride(l)
end do
cmfd % current(8,h,i,j,k) = t % results(RESULT_SUM,1,score_index)
! Bottom surface
matching_bins(i_filter_surf) = OUT_BOTTOM
score_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = OUT_BOTTOM
score_index = 1
do l = 1, size(t % filter)
score_index = score_index + (filter_matches(t % filter(l)) &
% bins % data(1) - 1) * t % stride(l)
end do
cmfd % current(9,h,i,j,k) = t % results(RESULT_SUM,1,score_index)
matching_bins(i_filter_surf) = IN_BOTTOM
score_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = IN_BOTTOM
score_index = 1
do l = 1, size(t % filter)
score_index = score_index + (filter_matches(t % filter(l)) &
% bins % data(1) - 1) * t % stride(l)
end do
cmfd % current(10,h,i,j,k) = t % results(RESULT_SUM,1,score_index)
! Top surface
matching_bins(i_filter_surf) = IN_TOP
score_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = IN_TOP
score_index = 1
do l = 1, size(t % filter)
score_index = score_index + (filter_matches(t % filter(l)) &
% bins % data(1) - 1) * t % stride(l)
end do
cmfd % current(11,h,i,j,k) = t % results(RESULT_SUM,1,score_index)
matching_bins(i_filter_surf) = OUT_TOP
score_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = OUT_TOP
score_index = 1
do l = 1, size(t % filter)
score_index = score_index + (filter_matches(t % filter(l)) &
% bins % data(1) - 1) * t % stride(l)
end do
cmfd % current(12,h,i,j,k) = t % results(RESULT_SUM,1,score_index)
end if TALLY

11
src/global.F90
View file

@ -22,7 +22,7 @@ module global
use surface_header, only: SurfaceContainer
use source_header, only: SourceDistribution
use tally_header, only: TallyObject, TallyDerivative
use tally_filter_header, only: TallyFilterContainer
use tally_filter_header, only: TallyFilterContainer, TallyFilterMatch
use trigger_header, only: KTrigger
use timer_header, only: Timer
use volume_header, only: VolumeCalculation
@ -143,8 +143,7 @@ module global
type(RegularMesh), allocatable, target :: meshes(:)
type(TallyObject), allocatable, target :: tallies(:)
type(TallyFilterContainer), allocatable, target :: filters(:)
integer, allocatable :: matching_bins(:)
real(8), allocatable :: filter_weights(:)
type(TallyFilterMatch), allocatable :: filter_matches(:)
! Pointers for different tallies
type(TallyObject), pointer :: user_tallies(:) => null()
@ -438,8 +437,7 @@ module global
type(Nuclide0K), allocatable, target :: nuclides_0K(:) ! 0K nuclides info
!$omp threadprivate(micro_xs, material_xs, fission_bank, n_bank, &
!$omp& trace, thread_id, current_work, matching_bins, &
!$omp& filter_weights)
!$omp& trace, thread_id, current_work, filter_matches)
contains
@ -500,8 +498,7 @@ contains
if (allocated(meshes)) deallocate(meshes)
if (allocated(filters)) deallocate(filters)
if (allocated(tallies)) deallocate(tallies)
if (allocated(matching_bins)) deallocate(matching_bins)
if (allocated(filter_weights)) deallocate(filter_weights)
if (allocated(filter_matches)) deallocate(filter_matches)
! Deallocate fission and source bank and entropy
!$omp parallel

2
src/initialize.F90
View file

@ -107,7 +107,7 @@ contains
call logarithmic_grid()
end if
! Allocate and setup tally stride, matching_bins, and tally maps
! Allocate and setup tally stride, filter_matches, and tally maps
call configure_tallies()
! Set up tally procedure pointers

165
src/output.F90
View file

@ -690,6 +690,7 @@ contains
integer :: k ! loop index for scoring bins
integer :: n ! loop index for nuclides
integer :: l ! loop index for user scores
integer :: h ! loop index for tally filters
integer :: indent ! number of spaces to preceed output
integer :: filter_index ! index in results array for filters
integer :: score_index ! scoring bin index
@ -802,7 +803,10 @@ contains
! to be used for a given tally.
! Initialize bins, filter level, and indentation
matching_bins(1:size(t % filter)) = 0
do h = 1, size(t % filter)
call filter_matches(t % filter(h)) % bins % clear()
call filter_matches(t % filter(h)) % bins % push_back(0)
end do
j = 1
indent = 0
@ -812,16 +816,18 @@ contains
if (size(t % filter) == 0) exit find_bin
! Increment bin combination
matching_bins(j) = matching_bins(j) + 1
filter_matches(j) % bins % data(1) = &
filter_matches(j) % bins % data(1) + 1
! =================================================================
! REACHED END OF BINS FOR THIS FILTER, MOVE TO NEXT FILTER
if (matching_bins(j) > filters(t % filter(j)) % obj % n_bins) then
if (filter_matches(j) % bins % data(1) > &
filters(t % filter(j)) % obj % n_bins) then
! If this is the first filter, then exit
if (j == 1) exit print_bin
matching_bins(j) = 0
filter_matches(j) % bins % data(1) = 0
j = j - 1
indent = indent - 2
@ -835,7 +841,7 @@ contains
! Print current filter information
write(UNIT=unit_tally, FMT='(1X,2A)') repeat(" ", indent), &
trim(filters(t % filter(j)) % obj % &
text_label(matching_bins(j)))
text_label(filter_matches(j) % bins % data(1)))
indent = indent + 2
j = j + 1
end if
@ -846,19 +852,18 @@ contains
if (size(t % filter) > 0) then
write(UNIT=unit_tally, FMT='(1X,2A)') repeat(" ", indent), &
trim(filters(t % filter(j)) % obj % &
text_label(matching_bins(j)))
text_label(filter_matches(j) % bins % data(1)))
end if
! Determine scoring index for this bin combination -- note that unlike
! in the score_tally subroutine, we have to use max(bins,1) since all
! bins below the lowest filter level will be zeros
if (size(t % filter) > 0) then
filter_index = sum((max(matching_bins(1:size(t % filter)),1) - 1) &
* t % stride) + 1
else
filter_index = 1
end if
filter_index = 1
do h = 1, size(t % filter)
filter_index = filter_index + (max(filter_matches(t % filter(h)) &
% bins % data(1),1) - 1) * t % stride(h)
end do
! Write results for this filter bin combination
score_index = 0
@ -958,6 +963,7 @@ contains
integer, intent(in) :: unit_tally
integer :: i ! mesh index
integer :: j ! loop index over tally filters
integer :: ijk(3) ! indices of mesh cells
integer :: n_dim ! number of mesh dimensions
integer :: n_cells ! number of mesh cells
@ -968,25 +974,30 @@ contains
integer :: n ! number of incoming energy bins
integer :: filter_index ! index in results array for filters
logical :: print_ebin ! should incoming energy bin be displayed?
logical :: energy_filters ! energy filters present
character(MAX_LINE_LEN) :: string
type(RegularMesh), pointer :: m
! Get pointer to mesh
i_filter_mesh = t % find_filter(FILTER_MESH)
i_filter_surf = t % find_filter(FILTER_SURFACE)
select type(filt => filters(t % filter(i_filter_mesh)) % obj)
i_filter_mesh = t % filter(t % find_filter(FILTER_MESH))
i_filter_surf = t % filter(t % find_filter(FILTER_SURFACE))
select type(filt => filters(i_filter_mesh) % obj)
type is (MeshFilter)
m => meshes(filt % mesh)
end select
! initialize bins array
matching_bins(1:size(t % filter)) = 1
do j = 1, size(t % filter)
call filter_matches(t % filter(j)) % bins % clear()
call filter_matches(t % filter(j)) % bins % push_back(1)
end do
! determine how many energy in bins there are
i_filter_ein = t % find_filter(FILTER_ENERGYIN)
if (i_filter_ein > 0) then
energy_filters = (t % find_filter(FILTER_ENERGYIN) > 0)
if (energy_filters) then
print_ebin = .true.
n = filters(t % filter(i_filter_ein)) % obj % n_bins
i_filter_ein = t % filter(t % find_filter(FILTER_ENERGYIN))
n = filters(i_filter_ein) % obj % n_bins
else
print_ebin = .false.
n = 1
@ -1001,7 +1012,7 @@ contains
! Get the indices for this cell
call bin_to_mesh_indices(m, i, ijk)
matching_bins(i_filter_mesh) = i
filter_matches(i_filter_mesh) % bins % data(1) = i
! Write the header for this cell
if (n_dim == 1) then
@ -1019,43 +1030,55 @@ contains
do l = 1, n
if (print_ebin) then
! Set incoming energy bin
matching_bins(i_filter_ein) = l
filter_matches(i_filter_ein) % bins % data(1) = l
! Write incoming energy bin
write(UNIT=unit_tally, FMT='(3X,A)') &
trim(filters(t % filter(i_filter_ein)) % obj % text_label( &
matching_bins(i_filter_ein)))
trim(filters(i_filter_ein) % obj % text_label( &
filter_matches(i_filter_ein) % bins % data(1)))
end if
! Left Surface
matching_bins(i_filter_surf) = OUT_LEFT
filter_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = OUT_LEFT
filter_index = 1
do j = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % filter(j)) &
% bins % data(1) - 1) * t % stride(j)
end do
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current on Left", &
to_str(t % results(RESULT_SUM,1,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index)))
matching_bins(i_filter_surf) = IN_LEFT
filter_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = IN_LEFT
filter_index = 1
do j = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % filter(j)) &
% bins % data(1) - 1) * t % stride(j)
end do
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current on Left", &
to_str(t % results(RESULT_SUM,1,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index)))
! Right Surface
matching_bins(i_filter_surf) = OUT_RIGHT
filter_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = OUT_RIGHT
filter_index = 1
do j = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % filter(j)) &
% bins % data(1) - 1) * t % stride(j)
end do
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current on Right", &
to_str(t % results(RESULT_SUM,1,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index)))
matching_bins(i_filter_surf) = IN_RIGHT
filter_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = IN_RIGHT
filter_index = 1
do j = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % filter(j)) &
% bins % data(1) - 1) * t % stride(j)
end do
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current on Right", &
to_str(t % results(RESULT_SUM,1,filter_index)), &
@ -1064,34 +1087,46 @@ contains
if (n_dim >= 2) then
! Back Surface
matching_bins(i_filter_surf) = OUT_BACK
filter_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = OUT_BACK
filter_index = 1
do j = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % filter(j)) &
% bins % data(1) - 1) * t % stride(j)
end do
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current on Back", &
to_str(t % results(RESULT_SUM,1,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index)))
matching_bins(i_filter_surf) = IN_BACK
filter_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = IN_BACK
filter_index = 1
do j = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % filter(j)) &
% bins % data(1) - 1) * t % stride(j)
end do
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current on Back", &
to_str(t % results(RESULT_SUM,1,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index)))
! Front Surface
matching_bins(i_filter_surf) = OUT_FRONT
filter_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = OUT_FRONT
filter_index = 1
do j = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % filter(j)) &
% bins % data(1) - 1) * t % stride(j)
end do
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Net Current on Front", &
to_str(t % results(RESULT_SUM,1,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index)))
matching_bins(i_filter_surf) = IN_FRONT
filter_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = IN_FRONT
filter_index = 1
do j = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % filter(j)) &
% bins % data(1) - 1) * t % stride(j)
end do
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Net Current on Front", &
to_str(t % results(RESULT_SUM,1,filter_index)), &
@ -1100,34 +1135,46 @@ contains
if (n_dim == 3) then
! Bottom Surface
matching_bins(i_filter_surf) = OUT_BOTTOM
filter_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = OUT_BOTTOM
filter_index = 1
do j = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % filter(j)) &
% bins % data(1) - 1) * t % stride(j)
end do
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current on Bottom", &
to_str(t % results(RESULT_SUM,1,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index)))
matching_bins(i_filter_surf) = IN_BOTTOM
filter_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = IN_BOTTOM
filter_index = 1
do j = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % filter(j)) &
% bins % data(1) - 1) * t % stride(j)
end do
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current on Bottom", &
to_str(t % results(RESULT_SUM,1,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index)))
! Top Surface
matching_bins(i_filter_surf) = OUT_TOP
filter_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = OUT_TOP
filter_index = 1
do j = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % filter(j)) &
% bins % data(1) - 1) * t % stride(j)
end do
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current on Top", &
to_str(t % results(RESULT_SUM,1,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index)))
matching_bins(i_filter_surf) = IN_TOP
filter_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = IN_TOP
filter_index = 1
do j = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % filter(j)) &
% bins % data(1) - 1) * t % stride(j)
end do
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current on Top", &
to_str(t % results(RESULT_SUM,1,filter_index)), &

525
src/tally.F90
View file

@ -2204,14 +2204,16 @@ contains
type(Particle), intent(in) :: p
integer :: i
integer :: i, j
integer :: i_tally
integer :: i_filt
integer :: i_bin
integer :: k ! loop index for nuclide bins
! position during the loop
integer :: filter_index ! single index for single bin
integer :: i_nuclide ! index in nuclides array
integer :: matching_bin ! next valid filter bin
real(8) :: filter_weight ! combined weight of all filters
logical :: finished ! found all valid bin combinations
type(TallyObject), pointer :: t
! A loop over all tallies is necessary because we need to simultaneously
@ -2222,15 +2224,29 @@ contains
i_tally = active_analog_tallies % data(i)
t => tallies(i_tally)
! Find the first bin in each filter. There may be more than one matching
! bin per filter, but we'll deal with those later.
do i_filt = 1, size(t % filter)
call filters(t % filter(i_filt)) % obj % get_next_bin(p, &
t % estimator, NO_BIN_FOUND, matching_bins(i_filt), &
filter_weights(i_filt))
! Find all valid bins in each filter if they have not already been found
! for a previous tally.
do j = 1, size(t % filter)
i_filt = t % filter(j)
if (.not. filter_matches(i_filt) % bins_present) then
call filter_matches(i_filt) % bins % clear()
call filter_matches(i_filt) % weights % clear()
matching_bin = NO_BIN_FOUND
do
call filters(i_filt) % obj % get_next_bin(p, t % estimator, &
matching_bin, matching_bin, filter_weight)
if (matching_bin == NO_BIN_FOUND) exit
call filter_matches(i_filt) % bins % push_back(matching_bin)
call filter_matches(i_filt) % weights % push_back(filter_weight)
end do
filter_matches(i_filt) % bins_present = .true.
end if
! If there are no valid bins for this filter, then there is nothing to
! score and we can move on to the next tally.
if (matching_bins(i_filt) == NO_BIN_FOUND) cycle TALLY_LOOP
if (filter_matches(i_filt) % bins % size() == 0) cycle TALLY_LOOP
! Set the index of the bin used in the first filter combination
filter_matches(i_filt) % i_bin = 1
end do
! ========================================================================
@ -2238,10 +2254,19 @@ contains
FILTER_LOOP: do
! Reset scoring index and weight
filter_index = 1
filter_weight = 1
! Determine scoring index and weight for this filter combination
filter_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_weight = product(filter_weights(:size(t % filter)))
do j = 1, size(t % filter)
i_filt = t % filter(j)
i_bin = filter_matches(i_filt) % i_bin
filter_index = filter_index + (filter_matches(i_filt) % bins % &
data(i_bin) - 1) * t % stride(j)
filter_weight = filter_weight * filter_matches(i_filt) % weights % &
data(i_bin)
end do
! ======================================================================
! Nuclide logic
@ -2294,31 +2319,25 @@ contains
! If there are no filters, then we are done.
if (size(t % filter) == 0) exit FILTER_LOOP
! Increment the filter bins, starting with the last filter. If we get a
! NO_BIN_FOUND for the last filter, it means we finished all valid bins
! for that filter, but next-to-last filter might have more than one
! valid bin so we need to increment that one as well, and so on.
do i_filt = size(t % filter), 1, -1
call filters(t % filter(i_filt)) % obj % get_next_bin(p, &
t % estimator, matching_bins(i_filt), matching_bins(i_filt), &
filter_weights(i_filt))
if (matching_bins(i_filt) /= NO_BIN_FOUND) exit
end do
! If we got all NO_BIN_FOUNDs, then we have finished all valid bins for
! each of the filters. Exit the loop.
if (all(matching_bins(:size(t % filter)) == NO_BIN_FOUND)) &
exit FILTER_LOOP
! Reset all the filters with NO_BIN_FOUND. This will set them back to
! their first valid bin.
do i_filt = 1, size(t % filter)
if (matching_bins(i_filt) == NO_BIN_FOUND) then
call filters(t % filter(i_filt)) % obj % get_next_bin(p, &
t % estimator, matching_bins(i_filt), matching_bins(i_filt), &
filter_weights(i_filt))
! Increment the filter bins, starting with the last filter to find the
! next valid bin combination
finished = .true.
do j = size(t % filter), 1, -1
i_filt = t % filter(j)
if (filter_matches(i_filt) % i_bin < filter_matches(i_filt) % &
bins % size()) then
filter_matches(i_filt) % i_bin = filter_matches(i_filt) % i_bin + 1
finished = .false.
exit
else
filter_matches(i_filt) % i_bin = 1
end if
end do
! Once we have finished all valid bins for each of the filters, exit
! the loop.
if (finished) exit FILTER_LOOP
end do FILTER_LOOP
! If the user has specified that we can assume all tallies are spatially
@ -2330,6 +2349,11 @@ contains
end do TALLY_LOOP
! Reset filter matches flag
do i = 1, n_filters
filter_matches(i) % bins_present = .false.
end do
! Reset tally map positioning
position = 0
@ -2339,15 +2363,17 @@ contains
type(Particle), intent(in) :: p
integer :: i, m
integer :: i, j, m
integer :: i_tally
integer :: i_filt
integer :: i_bin
integer :: k ! loop index for nuclide bins
! position during the loop
integer :: filter_index ! single index for single bin
integer :: i_nuclide ! index in nuclides array
integer :: matching_bin ! next valid filter bin
real(8) :: filter_weight ! combined weight of all filters
real(8) :: atom_density
logical :: finished ! found all valid bin combinations
type(TallyObject), pointer :: t
type(Material), pointer :: mat
@ -2363,15 +2389,29 @@ contains
! nuclides are in the material
mat => materials(p % material)
! Find the first bin in each filter. There may be more than one matching
! bin per filter, but we'll deal with those later.
do i_filt = 1, size(t % filter)
call filters(t % filter(i_filt)) % obj % get_next_bin(p, &
t % estimator, NO_BIN_FOUND, matching_bins(i_filt), &
filter_weights(i_filt))
! Find all valid bins in each filter if they have not already been found
! for a previous tally.
do j = 1, size(t % filter)
i_filt = t % filter(j)
if (.not. filter_matches(i_filt) % bins_present) then
call filter_matches(i_filt) % bins % clear()
call filter_matches(i_filt) % weights % clear()
matching_bin = NO_BIN_FOUND
do
call filters(i_filt) % obj % get_next_bin(p, t % estimator, &
matching_bin, matching_bin, filter_weight)
if (matching_bin == NO_BIN_FOUND) exit
call filter_matches(i_filt) % bins % push_back(matching_bin)
call filter_matches(i_filt) % weights % push_back(filter_weight)
end do
filter_matches(i_filt) % bins_present = .true.
end if
! If there are no valid bins for this filter, then there is nothing to
! score and we can move on to the next tally.
if (matching_bins(i_filt) == NO_BIN_FOUND) cycle TALLY_LOOP
if (filter_matches(i_filt) % bins % size() == 0) cycle TALLY_LOOP
! Set the index of the bin used in the first filter combination
filter_matches(i_filt) % i_bin = 1
end do
! ========================================================================
@ -2379,10 +2419,19 @@ contains
FILTER_LOOP: do
! Reset scoring index and weight
filter_index = 1
filter_weight = 1
! Determine scoring index and weight for this filter combination
filter_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_weight = product(filter_weights(:size(t % filter)))
do j = 1, size(t % filter)
i_filt = t % filter(j)
i_bin = filter_matches(i_filt) % i_bin
filter_index = filter_index + (filter_matches(i_filt) % bins % &
data(i_bin) - 1) * t % stride(j)
filter_weight = filter_weight * filter_matches(i_filt) % weights % &
data(i_bin)
end do
! ======================================================================
! Nuclide logic
@ -2423,31 +2472,25 @@ contains
! If there are no filters, then we are done.
if (size(t % filter) == 0) exit FILTER_LOOP
! Increment the filter bins, starting with the last filter. If we get a
! NO_BIN_FOUND for the last filter, it means we finished all valid bins
! for that filter, but next-to-last filter might have more than one
! valid bin so we need to increment that one as well, and so on.
do i_filt = size(t % filter), 1, -1
call filters(t % filter(i_filt)) % obj % get_next_bin(p, &
t % estimator, matching_bins(i_filt), matching_bins(i_filt), &
filter_weights(i_filt))
if (matching_bins(i_filt) /= NO_BIN_FOUND) exit
end do
! If we got all NO_BIN_FOUNDs, then we have finished all valid bins for
! each of the filters. Exit the loop.
if (all(matching_bins(:size(t % filter)) == NO_BIN_FOUND)) &
exit FILTER_LOOP
! Reset all the filters with NO_BIN_FOUND. This will set them back to
! their first valid bin.
do i_filt = 1, size(t % filter)
if (matching_bins(i_filt) == NO_BIN_FOUND) then
call filters(t % filter(i_filt)) % obj % get_next_bin(p, &
t % estimator, matching_bins(i_filt), matching_bins(i_filt), &
filter_weights(i_filt))
! Increment the filter bins, starting with the last filter to find the
! next valid bin combination
finished = .true.
do j = size(t % filter), 1, -1
i_filt = t % filter(j)
if (filter_matches(i_filt) % i_bin < filter_matches(i_filt) % &
bins % size()) then
filter_matches(i_filt) % i_bin = filter_matches(i_filt) % i_bin + 1
finished = .false.
exit
else
filter_matches(i_filt) % i_bin = 1
end if
end do
! Once we have finished all valid bins for each of the filters, exit
! the loop.
if (finished) exit FILTER_LOOP
end do FILTER_LOOP
! If the user has specified that we can assume all tallies are spatially
@ -2459,6 +2502,11 @@ contains
end do TALLY_LOOP
! Reset filter matches flag
do i = 1, n_filters
filter_matches(i) % bins_present = .false.
end do
! Reset tally map positioning
position = 0
@ -2485,6 +2533,10 @@ contains
integer :: d_bin ! delayed group bin index
integer :: n ! number of energies on filter
integer :: k ! loop index for bank sites
integer :: l ! loop index for tally filters
integer :: f ! index in filters array
integer :: b ! index of filter bin
integer :: i_bin ! index of matching filter bin
integer :: bin_energyout ! original outgoing energy bin
integer :: i_filter ! index for matching filter bin combination
real(8) :: filter_weight ! combined weight of all filters
@ -2493,11 +2545,12 @@ contains
integer :: g_out ! energy group of fission bank site
! save original outgoing energy bin and score index
i = t % find_filter(FILTER_ENERGYOUT)
bin_energyout = matching_bins(i)
i = t % filter(t % find_filter(FILTER_ENERGYOUT))
i_bin = filter_matches(i) % i_bin
bin_energyout = filter_matches(i) % bins % data(i_bin)
! declare the energyout filter type
select type(eo_filt => filters(t % filter(i)) % obj)
select type(eo_filt => filters(i) % obj)
type is (EnergyoutFilter)
! Get number of energies on filter
@ -2534,7 +2587,7 @@ contains
g_out = size(eo_filt % bins) - g_out
! change outgoing energy bin
matching_bins(i) = g_out
filter_matches(i) % bins % data(i_bin) = g_out
else
@ -2550,7 +2603,8 @@ contains
if (E_out < eo_filt % bins(1) .or. E_out > eo_filt % bins(n)) cycle
! change outgoing energy bin
matching_bins(i) = binary_search(eo_filt % bins, n, E_out)
filter_matches(i) % bins % data(i_bin) = &
binary_search(eo_filt % bins, n, E_out)
end if
@ -2559,8 +2613,11 @@ contains
(score_bin == SCORE_PROMPT_NU_FISSION .and. g == 0)) then
! determine scoring index and weight for this filter combination
i_filter = sum((matching_bins(1:size(t % filter)) - 1) * &
t % stride) + 1
i_filter = 1
do l = 1, size(t % filter)
i_filter = i_filter + (filter_matches(t % filter(l)) % bins % &
data(filter_matches(t % filter(l)) % i_bin) - 1) * t % stride(l)
end do
! Add score to tally
!$omp atomic
@ -2590,11 +2647,20 @@ contains
! the delayed group of this bin
if (d == g) then
! Reset scoring index and filter weight
i_filter = 1
filter_weight = 1
! determine scoring index and weight for this filter
! combination
i_filter = sum((matching_bins(1:size(t % filter)) - 1) * &
t % stride) + 1
filter_weight = product(filter_weights(:size(t % filter)))
do l = 1, size(t % filter)
f = t % filter(l)
b = filter_matches(f) % i_bin
i_filter = i_filter + (filter_matches(f) % bins % &
data(b) - 1) * t % stride(l)
filter_weight = filter_weight * filter_matches(f) % &
weights % data(b)
end do
call score_fission_delayed_dg(t, d_bin, &
score * filter_weight, i_score)
@ -2605,10 +2671,19 @@ contains
! if the delayed group filter is not present, add score to tally
else
! Reset scoring index and filter weight
i_filter = 1
filter_weight = 1
! determine scoring index and weight for this filter combination
i_filter = sum((matching_bins(1:size(t % filter)) - 1) * &
t % stride) + 1
filter_weight = product(filter_weights(:size(t % filter)))
do l = 1, size(t % filter)
f = t % filter(l)
b = filter_matches(f) % i_bin
i_filter = i_filter + (filter_matches(f) % bins % data(b) - 1) &
* t % stride(l)
filter_weight = filter_weight * filter_matches(f) % weights % &
data(b)
end do
! Add score to tally
!$omp atomic
@ -2620,7 +2695,7 @@ contains
end select
! reset outgoing energy bin and score index
matching_bins(i) = bin_energyout
filter_matches(i) % bins % data(i_bin) = bin_energyout
end subroutine score_fission_eout
@ -2636,23 +2711,31 @@ contains
real(8), intent(in) :: score ! actual score
integer, intent(in) :: score_index ! index for score
integer :: i ! loop over tally filters
integer :: i_filt ! index in filters array
integer :: i_bin ! index of matching filter bin
integer :: bin_original ! original bin index
integer :: filter_index ! index for matching filter bin combination
! save original delayed group bin
bin_original = matching_bins(t % find_filter(FILTER_DELAYEDGROUP))
matching_bins(t % find_filter(FILTER_DELAYEDGROUP)) = d_bin
i_filt = t % filter(t % find_filter(FILTER_DELAYEDGROUP))
i_bin = filter_matches(i_filt) % i_bin
bin_original = filter_matches(i_filt) % bins % data(i_bin)
filter_matches(i_filt) % bins % data(i_bin) = d_bin
! determine scoring index and weight on the modified matching_bins
filter_index = sum((matching_bins(1:size(t % filter)) - 1) * &
t % stride) + 1
! determine scoring index and weight on the modified matching bins
filter_index = 1
do i = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % filter(i)) % bins % &
data(filter_matches(t % filter(i)) % i_bin) - 1) * t % stride(i)
end do
!$omp atomic
t % results(RESULT_VALUE, score_index, filter_index) = &
t % results(RESULT_VALUE, score_index, filter_index) + score
! reset original delayed group bin
matching_bins(t % find_filter(FILTER_DELAYEDGROUP)) = bin_original
filter_matches(i_filt) % bins % data(i_bin) = bin_original
end subroutine score_fission_delayed_dg
@ -2671,13 +2754,16 @@ contains
integer :: i
integer :: i_tally
integer :: i_filt
integer :: i_bin
integer :: j ! loop index for scoring bins
integer :: k ! loop index for nuclide bins
integer :: filter_index ! single index for single bin
integer :: i_nuclide ! index in nuclides array (from bins)
integer :: matching_bin ! next valid filter bin
real(8) :: flux ! tracklength estimate of flux
real(8) :: atom_density ! atom density of single nuclide in atom/b-cm
real(8) :: filter_weight ! combined weight of all filters
logical :: finished ! found all valid bin combinations
type(TallyObject), pointer :: t
type(Material), pointer :: mat
@ -2692,15 +2778,29 @@ contains
i_tally = active_tracklength_tallies % data(i)
t => tallies(i_tally)
! Find the first bin in each filter. There may be more than one matching
! bin per filter, but we'll deal with those later.
do i_filt = 1, size(t % filter)
call filters(t % filter(i_filt)) % obj % get_next_bin(p, &
t % estimator, NO_BIN_FOUND, matching_bins(i_filt), &
filter_weights(i_filt))
! Find all valid bins in each filter if they have not already been found
! for a previous tally.
do j = 1, size(t % filter)
i_filt = t % filter(j)
if (.not. filter_matches(i_filt) % bins_present) then
call filter_matches(i_filt) % bins % clear()
call filter_matches(i_filt) % weights % clear()
matching_bin = NO_BIN_FOUND
do
call filters(i_filt) % obj % get_next_bin(p, t % estimator, &
matching_bin, matching_bin, filter_weight)
if (matching_bin == NO_BIN_FOUND) exit
call filter_matches(i_filt) % bins % push_back(matching_bin)
call filter_matches(i_filt) % weights % push_back(filter_weight)
end do
filter_matches(i_filt) % bins_present = .true.
end if
! If there are no valid bins for this filter, then there is nothing to
! score and we can move on to the next tally.
if (matching_bins(i_filt) == NO_BIN_FOUND) cycle TALLY_LOOP
if (filter_matches(i_filt) % bins % size() == 0) cycle TALLY_LOOP
! Set the index of the bin used in the first filter combination
filter_matches(i_filt) % i_bin = 1
end do
! ========================================================================
@ -2708,10 +2808,19 @@ contains
FILTER_LOOP: do
! Reset scoring index and weight
filter_index = 1
filter_weight = 1
! Determine scoring index and weight for this filter combination
filter_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_weight = product(filter_weights(:size(t % filter)))
do j = 1, size(t % filter)
i_filt = t % filter(j)
i_bin = filter_matches(i_filt) % i_bin
filter_index = filter_index + (filter_matches(i_filt) % bins % &
data(i_bin) - 1) * t % stride(j)
filter_weight = filter_weight * filter_matches(i_filt) % weights % &
data(i_bin)
end do
! ======================================================================
! Nuclide logic
@ -2765,31 +2874,25 @@ contains
! If there are no filters, then we are done.
if (size(t % filter) == 0) exit FILTER_LOOP
! Increment the filter bins, starting with the last filter. If we get a
! NO_BIN_FOUND for the last filter, it means we finished all valid bins
! for that filter, but next-to-last filter might have more than one
! valid bin so we need to increment that one as well, and so on.
do i_filt = size(t % filter), 1, -1
call filters(t % filter(i_filt)) % obj % get_next_bin(p, &
t % estimator, matching_bins(i_filt), matching_bins(i_filt), &
filter_weights(i_filt))
if (matching_bins(i_filt) /= NO_BIN_FOUND) exit
end do
! If we got all NO_BIN_FOUNDs, then we have finished all valid bins for
! each of the filters. Exit the loop.
if (all(matching_bins(:size(t % filter)) == NO_BIN_FOUND)) &
exit FILTER_LOOP
! Reset all the filters with NO_BIN_FOUND. This will set them back to
! their first valid bin.
do i_filt = 1, size(t % filter)
if (matching_bins(i_filt) == NO_BIN_FOUND) then
call filters(t % filter(i_filt)) % obj % get_next_bin(p, &
t % estimator, matching_bins(i_filt), matching_bins(i_filt), &
filter_weights(i_filt))
! Increment the filter bins, starting with the last filter to find the
! next valid bin combination
finished = .true.
do j = size(t % filter), 1, -1
i_filt = t % filter(j)
if (filter_matches(i_filt) % i_bin < filter_matches(i_filt) % &
bins % size()) then
filter_matches(i_filt) % i_bin = filter_matches(i_filt) % i_bin + 1
finished = .false.
exit
else
filter_matches(i_filt) % i_bin = 1
end if
end do
! Once we have finished all valid bins for each of the filters, exit
! the loop.
if (finished) exit FILTER_LOOP
end do FILTER_LOOP
! If the user has specified that we can assume all tallies are spatially
@ -2801,6 +2904,11 @@ contains
end do TALLY_LOOP
! Reset filter matches flag
do i = 1, n_filters
filter_matches(i) % bins_present = .false.
end do
! Reset tally map positioning
position = 0
@ -2821,14 +2929,17 @@ contains
integer :: i
integer :: i_tally
integer :: i_filt
integer :: i_bin
integer :: j ! loop index for scoring bins
integer :: k ! loop index for nuclide bins
integer :: filter_index ! single index for single bin
integer :: i_nuclide ! index in nuclides array (from bins)
integer :: matching_bin ! next valid filter bin
real(8) :: flux ! collision estimate of flux
real(8) :: atom_density ! atom density of single nuclide
! in atom/b-cm
real(8) :: filter_weight ! combined weight of all filters
logical :: finished ! found all valid bin combinations
type(TallyObject), pointer :: t
type(Material), pointer :: mat
@ -2848,15 +2959,29 @@ contains
i_tally = active_collision_tallies % data(i)
t => tallies(i_tally)
! Find the first bin in each filter. There may be more than one matching
! bin per filter, but we'll deal with those later.
do i_filt = 1, size(t % filter)
call filters(t % filter(i_filt)) % obj % get_next_bin(p, &
t % estimator, NO_BIN_FOUND, matching_bins(i_filt), &
filter_weights(i_filt))
! Find all valid bins in each filter if they have not already been found
! for a previous tally.
do j = 1, size(t % filter)
i_filt = t % filter(j)
if (.not. filter_matches(i_filt) % bins_present) then
call filter_matches(i_filt) % bins % clear()
call filter_matches(i_filt) % weights % clear()
matching_bin = NO_BIN_FOUND
do
call filters(i_filt) % obj % get_next_bin(p, t % estimator, &
matching_bin, matching_bin, filter_weight)
if (matching_bin == NO_BIN_FOUND) exit
call filter_matches(i_filt) % bins % push_back(matching_bin)
call filter_matches(i_filt) % weights % push_back(filter_weight)
end do
filter_matches(i_filt) % bins_present = .true.
end if
! If there are no valid bins for this filter, then there is nothing to
! score and we can move on to the next tally.
if (matching_bins(i_filt) == NO_BIN_FOUND) cycle TALLY_LOOP
if (filter_matches(i_filt) % bins % size() == 0) cycle TALLY_LOOP
! Set the index of the bin used in the first filter combination
filter_matches(i_filt) % i_bin = 1
end do
! ========================================================================
@ -2864,10 +2989,19 @@ contains
FILTER_LOOP: do
! Reset scoring index and weight
filter_index = 1
filter_weight = 1
! Determine scoring index and weight for this filter combination
filter_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_weight = product(filter_weights(:size(t % filter)))
do j = 1, size(t % filter)
i_filt = t % filter(j)
i_bin = filter_matches(i_filt) % i_bin
filter_index = filter_index + (filter_matches(i_filt) % bins % &
data(i_bin) - 1) * t % stride(j)
filter_weight = filter_weight * filter_matches(i_filt) % weights % &
data(i_bin)
end do
! ======================================================================
! Nuclide logic
@ -2921,31 +3055,25 @@ contains
! If there are no filters, then we are done.
if (size(t % filter) == 0) exit FILTER_LOOP
! Increment the filter bins, starting with the last filter. If we get a
! NO_BIN_FOUND for the last filter, it means we finished all valid bins
! for that filter, but next-to-last filter might have more than one
! valid bin so we need to increment that one as well, and so on.
do i_filt = size(t % filter), 1, -1
call filters(t % filter(i_filt)) % obj % get_next_bin(p, &
t % estimator, matching_bins(i_filt), matching_bins(i_filt), &
filter_weights(i_filt))
if (matching_bins(i_filt) /= NO_BIN_FOUND) exit
end do
! If we got all NO_BIN_FOUNDs, then we have finished all valid bins for
! each of the filters. Exit the loop.
if (all(matching_bins(:size(t % filter)) == NO_BIN_FOUND)) &
exit FILTER_LOOP
! Reset all the filters with NO_BIN_FOUND. This will set them back to
! their first valid bin.
do i_filt = 1, size(t % filter)
if (matching_bins(i_filt) == NO_BIN_FOUND) then
call filters(t % filter(i_filt)) % obj % get_next_bin(p, &
t % estimator, matching_bins(i_filt), matching_bins(i_filt), &
filter_weights(i_filt))
! Increment the filter bins, starting with the last filter to find the
! next valid bin combination
finished = .true.
do j = size(t % filter), 1, -1
i_filt = t % filter(j)
if (filter_matches(i_filt) % i_bin < filter_matches(i_filt) % &
bins % size()) then
filter_matches(i_filt) % i_bin = filter_matches(i_filt) % i_bin + 1
finished = .false.
exit
else
filter_matches(i_filt) % i_bin = 1
end if
end do
! Once we have finished all valid bins for each of the filters, exit
! the loop.
if (finished) exit FILTER_LOOP
end do FILTER_LOOP
! If the user has specified that we can assume all tallies are spatially
@ -2957,6 +3085,11 @@ contains
end do TALLY_LOOP
! Reset filter matches flag
do i = 1, n_filters
filter_matches(i) % bins_present = .false.
end do
! Reset tally map positioning
position = 0
@ -2973,7 +3106,7 @@ contains
integer :: i
integer :: i_tally
integer :: j ! loop indices
integer :: j, k ! loop indices
integer :: n_dim ! num dimensions of the mesh
integer :: d1 ! dimension index
integer :: d2 ! dimension index
@ -2992,9 +3125,11 @@ contains
real(8) :: d(3) ! distance to each bounding surface
real(8) :: distance ! actual distance traveled
real(8) :: filt_score ! score applied by filters
integer :: matching_bin ! next valid filter bin
logical :: start_in_mesh ! particle's starting xyz in mesh?
logical :: end_in_mesh ! particle's ending xyz in mesh?
logical :: cross_surface ! whether the particle crosses a surface
logical :: energy_filter ! energy filter present
type(TallyObject), pointer :: t
type(RegularMesh), pointer :: m
@ -3007,13 +3142,25 @@ contains
i_tally = active_current_tallies % data(i)
t => tallies(i_tally)
! Check for energy filter
energy_filter = (t % find_filter(FILTER_ENERGYIN) > 0)
! Get index for mesh, surface, and energy filters
i_filter_mesh = t % find_filter(FILTER_MESH)
i_filter_surf = t % find_filter(FILTER_SURFACE)
i_filter_energy = t % find_filter(FILTER_ENERGYIN)
i_filter_mesh = t % filter(t % find_filter(FILTER_MESH))
i_filter_surf = t % filter(t % find_filter(FILTER_SURFACE))
if (energy_filter) then
i_filter_energy = t % filter(t % find_filter(FILTER_ENERGYIN))
end if
! Reset the matching bins arrays
call filter_matches(i_filter_mesh) % bins % resize(1)
call filter_matches(i_filter_surf) % bins % resize(1)
if (energy_filter) then
call filter_matches(i_filter_energy) % bins % resize(1)
end if
! Get pointer to mesh
select type(filt => filters(t % filter(i_filter_mesh)) % obj)
select type(filt => filters(i_filter_mesh) % obj)
type is (MeshFilter)
m => meshes(filt % mesh)
end select
@ -3047,11 +3194,11 @@ contains
! Determine incoming energy bin. We need to tell the energy filter this
! is a tracklength tally so it uses the pre-collision energy.
if (i_filter_energy > 0) then
call filters(t % filter(i_filter_energy)) % obj % get_next_bin(p, &
ESTIMATOR_TRACKLENGTH, NO_BIN_FOUND, &
matching_bins(i_filter_energy), filt_score)
if (matching_bins(i_filter_energy) == NO_BIN_FOUND) cycle
if (energy_filter) then
call filters(i_filter_energy) % obj % get_next_bin(p, &
ESTIMATOR_TRACKLENGTH, NO_BIN_FOUND, matching_bin, filt_score)
if (matching_bin == NO_BIN_FOUND) cycle
filter_matches(i_filter_energy) % bins % data(1) = matching_bin
end if
! Bounding coordinates
@ -3065,7 +3212,7 @@ contains
do j = 1, n_cross
! Reset scoring bin index
matching_bins(i_filter_surf) = 0
filter_matches(i_filter_surf) % bins % data(1) = 0
! Set the distances to infinity
d = INFINITY
@ -3107,11 +3254,14 @@ contains
! Outward current on d1 max surface
if (all(ijk0(:n_dim) >= 1) .and. &
all(ijk0(:n_dim) <= m % dimension)) then
matching_bins(i_filter_surf) = d1 * 4 - 1
matching_bins(i_filter_mesh) = &
filter_matches(i_filter_surf) % bins % data(1) = d1 * 4 - 1
filter_matches(i_filter_mesh) % bins % data(1) = &
mesh_indices_to_bin(m, ijk0)
filter_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_index = 1
do k = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % &
filter(k)) % bins % data(1) - 1) * t % stride(k)
end do
!$omp atomic
t % results(RESULT_VALUE, 1, filter_index) = &
t % results(RESULT_VALUE, 1, filter_index) + p % wgt
@ -3143,11 +3293,14 @@ contains
! If the particle crossed the surface, tally the current
if (cross_surface) then
ijk0(d1) = ijk0(d1) + 1
matching_bins(i_filter_surf) = d1 * 4 - 2
matching_bins(i_filter_mesh) = &
filter_matches(i_filter_surf) % bins % data(1) = d1 * 4 - 2
filter_matches(i_filter_mesh) % bins % data(1) = &
mesh_indices_to_bin(m, ijk0)
filter_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_index = 1
do k = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % &
filter(k)) % bins % data(1) - 1) * t % stride(k)
end do
!$omp atomic
t % results(RESULT_VALUE, 1, filter_index) = &
t % results(RESULT_VALUE, 1, filter_index) + p % wgt
@ -3163,11 +3316,14 @@ contains
! Outward current on d1 min surface
if (all(ijk0(:n_dim) >= 1) .and. &
all(ijk0(:n_dim) <= m % dimension)) then
matching_bins(i_filter_surf) = d1 * 4 - 3
matching_bins(i_filter_mesh) = &
filter_matches(i_filter_surf) % bins % data(1) = d1 * 4 - 3
filter_matches(i_filter_mesh) % bins % data(1) = &
mesh_indices_to_bin(m, ijk0)
filter_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_index = 1
do k = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % &
filter(k)) % bins % data(1) - 1) * t % stride(k)
end do
!$omp atomic
t % results(RESULT_VALUE, 1, filter_index) = &
t % results(RESULT_VALUE, 1, filter_index) + p % wgt
@ -3199,11 +3355,14 @@ contains
! If the particle crossed the surface, tally the current
if (cross_surface) then
ijk0(d1) = ijk0(d1) - 1
matching_bins(i_filter_surf) = d1 * 4
matching_bins(i_filter_mesh) = &
filter_matches(i_filter_surf) % bins % data(1) = d1 * 4
filter_matches(i_filter_mesh) % bins % data(1) = &
mesh_indices_to_bin(m, ijk0)
filter_index = sum((matching_bins(1:size(t % filter)) - 1) &
* t % stride) + 1
filter_index = 1
do k = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % &
filter(k)) % bins % data(1) - 1) * t % stride(k)
end do
!$omp atomic
t % results(RESULT_VALUE, 1, filter_index) = &
t % results(RESULT_VALUE, 1, filter_index) + p % wgt

14
src/tally_filter_header.F90
View file

@ -2,11 +2,25 @@ module tally_filter_header
use constants, only: MAX_LINE_LEN
use particle_header, only: Particle
use stl_vector, only: VectorInt, VectorReal
use hdf5
implicit none
!===============================================================================
! TALLYFILTERMATCH stores every valid bin and weight for a filter
!===============================================================================
type TallyFilterMatch
integer :: i_bin
type(VectorInt) :: bins
type(VectorReal) :: weights
! Indicates whether all valid bins for this filter have been found
logical :: bins_present = .false.
end type TallyFilterMatch
!===============================================================================
! TALLYFILTER describes a filter that limits what events score to a tally. For
! example, a cell filter indicates that only particles in a specified cell

10
src/tally_initialize.F90
View file

@ -29,7 +29,7 @@ contains
!===============================================================================
! SETUP_TALLY_ARRAYS allocates and populates several member arrays of the
! TallyObject derived type, including stride, matching_bins, and results.
! TallyObject derived type, including stride, filter_matches, and results.
!===============================================================================
subroutine setup_tally_arrays()
@ -38,21 +38,18 @@ contains
integer :: j ! loop index for filters
integer :: n ! temporary stride
integer :: i_filt ! filter index
integer :: max_n_filters = 0 ! maximum number of filters
type(TallyObject), pointer :: t
TALLY_LOOP: do i = 1, n_tallies
! Get pointer to tally
t => tallies(i)
! Allocate stride and matching_bins arrays
! Allocate stride
allocate(t % stride(size(t % filter)))
max_n_filters = max(max_n_filters, size(t % filter))
! 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
! stride
n = 1
STRIDE: do j = size(t % filter), 1, -1
i_filt = t % filter(j)
@ -72,8 +69,7 @@ contains
! Allocate array for matching filter bins
!$omp parallel
allocate(matching_bins(max_n_filters))
allocate(filter_weights(max_n_filters))
allocate(filter_matches(n_filters))
!$omp end parallel
end subroutine setup_tally_arrays

77
src/trigger.F90
View file

@ -94,6 +94,7 @@ contains
character(len=52), intent(inout) :: name ! "eigenvalue" or tally score
integer :: i ! index in tallies array
integer :: j ! index in tally filters
integer :: n ! loop index for nuclides
integer :: s ! loop index for triggers
integer :: filter_index ! index in results array for filters
@ -167,7 +168,10 @@ contains
else
! Initialize bins, filter level
matching_bins(1:size(t % filter)) = 0
do j = 1, size(t % filter)
call filter_matches(t % filter(j)) % bins % clear()
call filter_matches(t % filter(j)) % bins % push_back(0)
end do
FILTER_LOOP: do filter_index = 1, t % total_filter_bins
@ -284,6 +288,7 @@ contains
subroutine compute_tally_current(t, trigger)
integer :: i ! mesh index
integer :: j ! loop index for tally filters
integer :: ijk(3) ! indices of mesh cells
integer :: n_dim ! number of mesh dimensions
integer :: n_cells ! number of mesh cells
@ -301,21 +306,25 @@ contains
type(RegularMesh), pointer :: m ! surface current mesh
! Get pointer to mesh
i_filter_mesh = t % find_filter(FILTER_MESH)
i_filter_surf = t % find_filter(FILTER_SURFACE)
select type(filt => filters(t % filter(i_filter_mesh)) % obj)
i_filter_mesh = t % filter(t % find_filter(FILTER_MESH))
i_filter_surf = t % filter(t % find_filter(FILTER_SURFACE))
select type(filt => filters(i_filter_mesh) % obj)
type is (MeshFilter)
m => meshes(filt % mesh)
end select
! initialize bins array
matching_bins(1:size(t % filter)) = 1
do j = 1, size(t % filter)
call filter_matches(t % filter(j)) % bins % clear()
call filter_matches(t % filter(j)) % bins % push_back(1)
end do
! determine how many energyin bins there are
i_filter_ein = t % find_filter(FILTER_ENERGYIN)
if (i_filter_ein > 0) then
print_ebin = .true.
n = filters(t % filter(i_filter_ein)) % obj % n_bins
i_filter_ein = t % filter(i_filter_ein)
else
print_ebin = .false.
n = 1
@ -330,18 +339,21 @@ contains
! Get the indices for this cell
call bin_to_mesh_indices(m, i, ijk)
matching_bins(i_filter_mesh) = i
filter_matches(i_filter_mesh) % bins % data(1) = i
do l = 1, n
if (print_ebin) then
matching_bins(i_filter_ein) = l
filter_matches(i_filter_ein) % bins % data(1) = l
end if
! Left Surface
matching_bins(i_filter_surf) = OUT_LEFT
filter_index = &
sum((matching_bins(1:size(t % filter)) - 1) * t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = OUT_LEFT
filter_index = 1
do j = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % filter(j)) % &
bins % data(1) - 1) * t % stride(j)
end do
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -352,9 +364,12 @@ contains
trigger % variance = std_dev**2
! Right Surface
matching_bins(i_filter_surf) = OUT_RIGHT
filter_index = &
sum((matching_bins(1:size(t % filter)) - 1) * t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = OUT_RIGHT
filter_index = 1
do j = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % filter(j)) % &
bins % data(1) - 1) * t % stride(j)
end do
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -365,9 +380,12 @@ contains
trigger % variance = trigger % std_dev**2
! Back Surface
matching_bins(i_filter_surf) = OUT_BACK
filter_index = &
sum((matching_bins(1:size(t % filter)) - 1) * t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = OUT_BACK
filter_index = 1
do j = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % filter(j)) % &
bins % data(1) - 1) * t % stride(j)
end do
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -378,9 +396,12 @@ contains
trigger % variance = trigger % std_dev**2
! Front Surface
matching_bins(i_filter_surf) = OUT_FRONT
filter_index = &
sum((matching_bins(1:size(t % filter)) - 1) * t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = OUT_FRONT
filter_index = 1
do j = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % filter(j)) % &
bins % data(1) - 1) * t % stride(j)
end do
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -391,9 +412,12 @@ contains
trigger % variance = trigger % std_dev**2
! Bottom Surface
matching_bins(i_filter_surf) = OUT_BOTTOM
filter_index = &
sum((matching_bins(1:size(t % filter)) - 1) * t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = OUT_BOTTOM
filter_index = 1
do j = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % filter(j)) % &
bins % data(1) - 1) * t % stride(j)
end do
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -404,9 +428,12 @@ contains
trigger % variance = trigger % std_dev**2
! Top Surface
matching_bins(i_filter_surf) = OUT_TOP
filter_index = &
sum((matching_bins(1:size(t % filter)) - 1) * t % stride) + 1
filter_matches(i_filter_surf) % bins % data(1) = OUT_TOP
filter_index = 1
do j = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % filter(j)) % &
bins % data(1) - 1) * t % stride(j)
end do
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev

1
tests/test_diff_tally/test_diff_tally.py
View file

@ -127,7 +127,6 @@ class DiffTallyTestHarness(PyAPITestHarness):
# Extract the relevant data as a CSV string.
cols = ('d_material', 'd_nuclide', 'd_variable', 'score', 'mean',
'std. dev.')
print(df.to_csv(None, columns=cols, index=False, float_format='%.7e'))
return df.to_csv(None, columns=cols, index=False, float_format='%.7e')