diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index a2a642fc7a..2e304a9096 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -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 diff --git a/src/global.F90 b/src/global.F90 index dd013b3a89..ee308bbef3 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -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 diff --git a/src/initialize.F90 b/src/initialize.F90 index 566b5aad7d..e9718fee10 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -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 diff --git a/src/output.F90 b/src/output.F90 index 344b6c4607..7dc4ef4751 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -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)), & diff --git a/src/tally.F90 b/src/tally.F90 index 02ab07339e..fab30485dc 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -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 diff --git a/src/tally_filter_header.F90 b/src/tally_filter_header.F90 index d934f63ffc..350070e574 100644 --- a/src/tally_filter_header.F90 +++ b/src/tally_filter_header.F90 @@ -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 diff --git a/src/tally_initialize.F90 b/src/tally_initialize.F90 index 8baabcd828..46c0283642 100644 --- a/src/tally_initialize.F90 +++ b/src/tally_initialize.F90 @@ -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 diff --git a/src/trigger.F90 b/src/trigger.F90 index 32f2549c8f..9eb39e9cc2 100644 --- a/src/trigger.F90 +++ b/src/trigger.F90 @@ -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 diff --git a/tests/test_diff_tally/test_diff_tally.py b/tests/test_diff_tally/test_diff_tally.py index 51e4cd5ebd..84aa181568 100644 --- a/tests/test_diff_tally/test_diff_tally.py +++ b/tests/test_diff_tally/test_diff_tally.py @@ -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')