diff --git a/src/constants.f90 b/src/constants.f90 index 9d7f46c410..09caced703 100644 --- a/src/constants.f90 +++ b/src/constants.f90 @@ -3,8 +3,8 @@ module constants implicit none ! Versioning numbers - integer, parameter :: VERSION_MAJOR = 0 - integer, parameter :: VERSION_MINOR = 3 + integer, parameter :: VERSION_MAJOR = 0 + integer, parameter :: VERSION_MINOR = 3 integer, parameter :: VERSION_RELEASE = 2 ! Physical constants @@ -163,18 +163,20 @@ module constants ! Tally map bin finding integer, parameter :: NO_BIN_FOUND = -1 - ! Tally map types - integer, parameter :: TALLY_MAP_TYPES = 6 - integer, parameter :: & - MAP_CELL = 1, & - MAP_SURFACE = 2, & - MAP_UNIVERSE = 3, & - MAP_MATERIAL = 4, & - MAP_MESH = 5, & - MAP_BORNIN = 6 + ! Tally filter and map types + integer, parameter :: TALLY_TYPES = 8 + integer, parameter :: & + T_UNIVERSE = 1, & + T_MATERIAL = 2, & + T_CELL = 3, & + T_CELLBORN = 4, & + T_SURFACE = 5, & + T_MESH = 6, & + T_ENERGYIN = 7, & + T_ENERGYOUT = 8 ! Fission neutron emission (nu) type - integer, parameter :: & + integer, parameter :: & NU_NONE = 0, & ! No nu values (non-fissionable) NU_POLYNOMIAL = 1, & ! Nu values given by polynomial NU_TABULAR = 2 ! Nu values given by tabular distribution @@ -186,7 +188,8 @@ module constants integer, parameter :: MAX_WORD_LEN = 150 ! Unit numbers - integer, parameter :: UNIT_LOG = 9 ! unit # for writing log file + integer, parameter :: UNIT_LOG = 11 ! unit # for writing log file + integer, parameter :: UNIT_TALLY = 12 ! unit # for writing tally file end module constants diff --git a/src/initialize.f90 b/src/initialize.f90 index bf32c1307b..ddbcf00bf9 100644 --- a/src/initialize.f90 +++ b/src/initialize.f90 @@ -397,13 +397,13 @@ contains end if ! ======================================================================= - ! ADJUST BORNIN INDICES FOR EACH TALLY + ! ADJUST CELLBORN INDICES FOR EACH TALLY - if (associated(t % bornin_bins)) then - do j = 1, size(t % bornin_bins) - uid = t % bornin_bins(j) % scalar + if (associated(t % cellborn_bins)) then + do j = 1, size(t % cellborn_bins) + uid = t % cellborn_bins(j) % scalar if (dict_has_key(cell_dict, uid)) then - t % bornin_bins(j) % scalar = dict_get_key(cell_dict, uid) + t % cellborn_bins(j) % scalar = dict_get_key(cell_dict, uid) else msg = "Could not find material " // trim(int_to_str(uid)) // & & " specified on tally " // trim(int_to_str(t % uid)) diff --git a/src/input_xml.f90 b/src/input_xml.f90 index b0938fc118..7fb12d6f54 100644 --- a/src/input_xml.f90 +++ b/src/input_xml.f90 @@ -534,6 +534,14 @@ contains do i = 1, n_tallies t => tallies(i) + ! Allocate arrays for number of bins and stride in scores array + allocate(t % n_bins(TALLY_TYPES)) + allocate(t % stride(TALLY_TYPES)) + + ! Initialize number of bins and stride + t % n_bins = 0 + t % stride = 0 + ! Copy material uid t % uid = tally_(i) % id @@ -554,7 +562,7 @@ contains do j = 1, n_words t % cell_bins(j) % scalar = str_to_int(words(j)) end do - t % n_cell_bins = n_words + t % n_bins(T_CELL) = n_words end if ! Read surface filter bins @@ -564,7 +572,7 @@ contains do j = 1, n_words t % surface_bins(j) % scalar = str_to_int(words(j)) end do - t % n_surface_bins = n_words + t % n_bins(T_SURFACE) = n_words end if ! Read universe filter bins @@ -574,7 +582,7 @@ contains do j = 1, n_words t % universe_bins(j) % scalar = str_to_int(words(j)) end do - t % n_universe_bins = n_words + t % n_bins(T_UNIVERSE) = n_words end if ! Read material filter bins @@ -584,7 +592,7 @@ contains do j = 1, n_words t % material_bins(j) % scalar = str_to_int(words(j)) end do - t % n_material_bins = n_words + t % n_bins(T_MATERIAL) = n_words end if ! Read mesh filter bins @@ -599,11 +607,11 @@ contains ! Read birth region filter bins if (len_trim(tally_(i) % filters % bornin) > 0) then call split_string(tally_(i) % filters % bornin, words, n_words) - allocate(t % bornin_bins(n_words)) + allocate(t % cellborn_bins(n_words)) do j = 1, n_words - t % bornin_bins(j) % scalar = str_to_int(words(j)) + t % cellborn_bins(j) % scalar = str_to_int(words(j)) end do - t % n_bornin_bins = n_words + t % n_bins(T_CELLBORN) = n_words end if ! Read incoming energy filter bins @@ -613,7 +621,7 @@ contains do j = 1, n_words t % energy_in(j) = str_to_real(words(j)) end do - t % n_energy_in = n_words - 1 + t % n_bins(T_ENERGYIN) = n_words - 1 end if ! Read outgoing energy filter bins @@ -623,7 +631,7 @@ contains do j = 1, n_words t % energy_out(j) = str_to_real(words(j)) end do - t % n_energy_out = n_words - 1 + t % n_bins(T_ENERGYOUT) = n_words - 1 end if ! Read macro reactions diff --git a/src/output.f90 b/src/output.f90 index 6bea82579c..2723ec8f64 100644 --- a/src/output.f90 +++ b/src/output.f90 @@ -516,9 +516,9 @@ contains write(ou,*) 'Tally ' // int_to_str(t % uid) - if (t % n_cell_bins > 0) then + if (t % n_bins(T_CELL) > 0) then string = "" - do i = 1, t % n_cell_bins + do i = 1, t % n_bins(T_CELL) uid = t % cell_bins(i) % scalar c => cells(uid) string = trim(string) // ' ' // trim(int_to_str(c % uid)) @@ -526,9 +526,9 @@ contains write(ou, *) ' Cell Bins:' // trim(string) end if - if (t % n_surface_bins > 0) then + if (t % n_bins(T_SURFACE) > 0) then string = "" - do i = 1, t % n_surface_bins + do i = 1, t % n_bins(T_SURFACE) uid = t % surface_bins(i) % scalar s => surfaces(uid) string = trim(string) // ' ' // trim(int_to_str(s % uid)) @@ -536,9 +536,9 @@ contains write(ou, *) ' Surface Bins:' // trim(string) end if - if (t % n_universe_bins) then + if (t % n_bins(T_UNIVERSE)) then string = "" - do i = 1, t % n_universe_bins + do i = 1, t % n_bins(T_UNIVERSE) uid = t % universe_bins(i) % scalar u => universes(uid) string = trim(string) // ' ' // trim(int_to_str(u % uid)) @@ -546,9 +546,9 @@ contains write(ou, *) ' Material Bins:' // trim(string) end if - if (t % n_material_bins) then + if (t % n_bins(T_MATERIAL)) then string = "" - do i = 1, t % n_material_bins + do i = 1, t % n_bins(T_MATERIAL) uid = t % material_bins(i) % scalar m => materials(uid) string = trim(string) // ' ' // trim(int_to_str(m % uid)) @@ -565,28 +565,28 @@ contains write(ou, *) ' Mesh Bins:' // trim(string) end if - if (t % n_bornin_bins > 0) then + if (t % n_bins(T_CELLBORN) > 0) then string = "" - do i = 1, t % n_bornin_bins - uid = t % bornin_bins(i) % scalar + do i = 1, t % n_bins(T_CELLBORN) + uid = t % cellborn_bins(i) % scalar c => cells(uid) string = trim(string) // ' ' // trim(int_to_str(c % uid)) end do write(ou, *) ' Birth Region Bins:' // trim(string) end if - if (t % n_energy_in > 0) then + if (t % n_bins(T_ENERGYIN) > 0) then string = "" - do i = 1, t % n_energy_in + 1 + do i = 1, t % n_bins(T_ENERGYIN) + 1 string = trim(string) // ' ' // trim(real_to_str(& t % energy_in(i))) end do write(ou,*) ' Incoming Energy Bins:' // trim(string) end if - if (t % n_energy_out > 0) then + if (t % n_bins(T_ENERGYOUT) > 0) then string = "" - do i = 1, t % n_energy_out + 1 + do i = 1, t % n_bins(T_ENERGYOUT) + 1 string = trim(string) // ' ' // trim(real_to_str(& t % energy_out(i))) end do diff --git a/src/tally.f90 b/src/tally.f90 index d30f26f165..451975287f 100644 --- a/src/tally.f90 +++ b/src/tally.f90 @@ -100,19 +100,20 @@ contains character(MAX_LINE_LEN) :: msg ! output/error message type(TallyObject), pointer :: t => null() - ! allocate tally map array - allocate(tally_maps(TALLY_MAP_TYPES)) + ! allocate tally map array -- note that we don't need a tally map for the + ! energy_in and energy_out filters + allocate(tally_maps(TALLY_TYPES - 2)) ! allocate list of items for each different filter type - allocate(tally_maps(MAP_CELL) % items(n_cells)) - allocate(tally_maps(MAP_SURFACE) % items(n_surfaces)) - allocate(tally_maps(MAP_UNIVERSE) % items(n_universes)) - allocate(tally_maps(MAP_MATERIAL) % items(n_materials)) - allocate(tally_maps(MAP_MESH) % items(100)) ! TODO: Change this - allocate(tally_maps(MAP_BORNIN) % items(n_cells)) + allocate(tally_maps(T_CELL) % items(n_cells)) + allocate(tally_maps(T_SURFACE) % items(n_surfaces)) + allocate(tally_maps(T_UNIVERSE) % items(n_universes)) + allocate(tally_maps(T_MATERIAL) % items(n_materials)) + allocate(tally_maps(T_MESH) % items(100)) ! TODO: Change this + allocate(tally_maps(T_CELLBORN) % items(n_cells)) ! Allocate and initialize tally map positioning for finding bins - allocate(position(TALLY_MAP_TYPES)) + allocate(position(TALLY_TYPES)) position = 0 do i = 1, n_tallies @@ -121,67 +122,74 @@ contains ! initialize number of scoring bins filter_bins = 1 - ! Add map elements for cell bins - n = t % n_cell_bins + ! determine if there are subdivisions for incoming or outgoing energy to + ! adjust the number of filter bins appropriately + n = t % n_bins(T_ENERGYOUT) + t % stride(T_ENERGYOUT) = filter_bins if (n > 0) then - do j = 1, n - index = t % cell_bins(j) % scalar - call add_map_element(tally_maps(MAP_CELL) % items(index), i, j) - end do filter_bins = filter_bins * n end if - ! Add map elements for surface bins - n = t % n_surface_bins + n = t % n_bins(T_ENERGYIN) + t % stride(T_ENERGYIN) = filter_bins if (n > 0) then - do j = 1, n - index = t % surface_bins(j) % scalar - call add_map_element(tally_maps(MAP_SURFACE) % items(index), i, j) - end do - filter_bins = filter_bins * n - end if - - ! Add map elements for universe bins - n = t % n_universe_bins - if (n > 0) then - do j = 1, n - index = t % universe_bins(j) % scalar - call add_map_element(tally_maps(MAP_UNIVERSE) % items(index), i, j) - end do - filter_bins = filter_bins * n - end if - - ! Add map elements for material bins - n = t % n_material_bins - if (n > 0) then - do j = 1, n - index = t % material_bins(j) % scalar - call add_map_element(tally_maps(MAP_MATERIAL) % items(index), i, j) - end do - filter_bins = filter_bins * n - end if - - ! Add map elements for bornin bins - n = t % n_bornin_bins - if (n > 0) then - do j = 1, n - index = t % bornin_bins(j) % scalar - call add_map_element(tally_maps(MAP_BORNIN) % items(index), i, j) - end do filter_bins = filter_bins * n end if ! TODO: Determine size of mesh to increase number of scoring bins - ! determine if there are subdivisions for incoming or outgoing energy to - ! adjust the number of scoring bins appropriately - n = t % n_energy_in + ! Add map elements for surface bins + n = t % n_bins(T_SURFACE) + t % stride(T_SURFACE) = filter_bins if (n > 0) then + do j = 1, n + index = t % surface_bins(j) % scalar + call add_map_element(tally_maps(T_SURFACE) % items(index), i, j) + end do filter_bins = filter_bins * n end if - n = t % n_energy_out + ! Add map elements for cellborn bins + n = t % n_bins(T_CELLBORN) + t % stride(T_CELLBORN) = filter_bins if (n > 0) then + do j = 1, n + index = t % cellborn_bins(j) % scalar + call add_map_element(tally_maps(T_CELLBORN) % items(index), i, j) + end do + filter_bins = filter_bins * n + end if + + ! Add map elements for cell bins + n = t % n_bins(T_CELL) + t % stride(T_CELL) = filter_bins + if (n > 0) then + do j = 1, n + index = t % cell_bins(j) % scalar + call add_map_element(tally_maps(T_CELL) % items(index), i, j) + end do + filter_bins = filter_bins * n + end if + + ! Add map elements for material bins + n = t % n_bins(T_MATERIAL) + t % stride(T_MATERIAL) = filter_bins + if (n > 0) then + do j = 1, n + index = t % material_bins(j) % scalar + call add_map_element(tally_maps(T_MATERIAL) % items(index), i, j) + end do + filter_bins = filter_bins * n + end if + + ! Add map elements for universe bins + n = t % n_bins(T_UNIVERSE) + t % stride(T_UNIVERSE) = filter_bins + if (n > 0) then + do j = 1, n + index = t % universe_bins(j) % scalar + call add_map_element(tally_maps(T_UNIVERSE) % items(index), i, j) + end do filter_bins = filter_bins * n end if @@ -249,13 +257,7 @@ contains integer :: i integer :: j integer :: n - integer :: cell_bin - integer :: surface_bin - integer :: universe_bin - integer :: material_bin - integer :: bornin_bin - integer :: energyin_bin - integer :: energyout_bin + integer :: bins(TALLY_TYPES) integer :: score_index real(8) :: score @@ -271,67 +273,67 @@ contains ! DETERMINE SCORING BIN COMBINATION ! determine next cell bin - if (t % n_cell_bins > 0) then - cell_bin = get_next_bin(MAP_CELL, p % cell, i) - if (cell_bin == NO_BIN_FOUND) cycle + if (t % n_bins(T_CELL) > 0) then + bins(T_CELL) = get_next_bin(T_CELL, p % cell, i) + if (bins(T_CELL) == NO_BIN_FOUND) cycle else - cell_bin = 1 + bins(T_CELL) = 1 end if ! determine next surface bin - if (t % n_surface_bins > 0) then - surface_bin = get_next_bin(MAP_SURFACE, p % surface, i) - if (surface_bin == NO_BIN_FOUND) cycle + if (t % n_bins(T_SURFACE) > 0) then + bins(T_SURFACE) = get_next_bin(T_SURFACE, p % surface, i) + if (bins(T_SURFACE) == NO_BIN_FOUND) cycle else - surface_bin = 1 + bins(T_SURFACE) = 1 end if ! determine next universe bin - if (t % n_universe_bins > 0) then - universe_bin = get_next_bin(MAP_UNIVERSE, p % universe, i) - if (universe_bin == NO_BIN_FOUND) cycle + if (t % n_bins(T_UNIVERSE) > 0) then + bins(T_UNIVERSE) = get_next_bin(T_UNIVERSE, p % universe, i) + if (bins(T_UNIVERSE) == NO_BIN_FOUND) cycle else - universe_bin = 1 + bins(T_UNIVERSE) = 1 end if ! determine next material bin - if (t % n_material_bins > 0) then - material_bin = get_next_bin(MAP_MATERIAL, p % material, i) - if (material_bin == NO_BIN_FOUND) cycle + if (t % n_bins(T_MATERIAL) > 0) then + bins(T_MATERIAL) = get_next_bin(T_MATERIAL, p % material, i) + if (bins(T_MATERIAL) == NO_BIN_FOUND) cycle else - material_bin = 1 + bins(T_MATERIAL) = 1 end if - ! determine next bornin bin - if (t % n_bornin_bins > 0) then - bornin_bin = get_next_bin(MAP_BORNIN, p % cell_born, i) - if (bornin_bin == NO_BIN_FOUND) cycle + ! determine next cellborn bin + if (t % n_bins(T_CELLBORN) > 0) then + bins(T_CELLBORN) = get_next_bin(T_CELLBORN, p % cell_born, i) + if (bins(T_CELLBORN) == NO_BIN_FOUND) cycle else - bornin_bin = 1 + bins(T_CELLBORN) = 1 end if ! determine incoming energy bin - n = t % n_energy_in + n = t % n_bins(T_ENERGYIN) if (n > 0) then ! check if energy of the particle is within energy bins if (p % E < t % energy_in(1) .or. p % E > t % energy_in(n)) cycle ! search to find incoming energy bin - energyin_bin = binary_search(t % energy_in, n, p % E) + bins(T_ENERGYIN) = binary_search(t % energy_in, n, p % E) else - energyin_bin = 1 + bins(T_ENERGYIN) = 1 end if ! determine outgoing energy bin - n = t % n_energy_out + n = t % n_bins(T_ENERGYOUT) if (n > 0) then ! check if energy of the particle is within energy bins if (p % E < t % energy_out(1) .or. p % E > t % energy_out(n)) cycle ! search to find incoming energy bin - energyout_bin = binary_search(t % energy_out, n, p % E) + bins(T_ENERGYOUT) = binary_search(t % energy_out, n, p % E) else - energyout_bin = 1 + bins(T_ENERGYOUT) = 1 end if ! ======================================================================= @@ -341,8 +343,7 @@ contains ! tally -- now we need to determine where in the scores array we should ! be accumulating the tally values - score_index = bins_to_index(t, cell_bin, surface_bin, universe_bin, & - material_bin, bornin_bin, energyin_bin, energyout_bin) + score_index = sum((bins - 1) * t % stride) + 1 ! ======================================================================= ! CALCULATE SCORES AND ACCUMULATE TALLY @@ -427,65 +428,6 @@ contains end function get_next_bin -!=============================================================================== -! BINS_TO_INDEX takes a combination of cell/surface/etc bins and returns the -! index in the scores array for the given tally -!=============================================================================== - - function bins_to_index(t, c_bin, s_bin, u_bin, m_bin, b_bin, & - ei_bin, eo_bin) result (index) - - type(TallyObject), pointer :: t - integer, intent(in) :: c_bin - integer, intent(in) :: s_bin - integer, intent(in) :: u_bin - integer, intent(in) :: m_bin - integer, intent(in) :: b_bin - integer, intent(in) :: ei_bin - integer, intent(in) :: eo_bin - integer :: index - - integer :: product - - index = 0 - product = 1 - - index = index + (c_bin - 1) - if (t % n_cell_bins > 0) then - product = product * t % n_cell_bins - end if - - index = index + (s_bin - 1) * product - if (t % n_surface_bins > 0) then - product = product * t % n_surface_bins - end if - - index = index + (u_bin - 1) * product - if (t % n_universe_bins > 0) then - product = product * t % n_universe_bins - end if - - index = index + (m_bin - 1) * product - if (t % n_material_bins > 0) then - product = product * t % n_material_bins - end if - - index = index + (b_bin - 1) * product - if (t % n_bornin_bins > 0) then - product = product * t % n_bornin_bins - end if - - index = index + (ei_bin - 1) * product - if (t % n_energy_in > 9) then - product = product * t % n_energy_in - end if - - ! We need to shift the index by one since the array in fortran starts at - ! unity, not zero - index = index + (eo_bin - 1) * product + 1 - - end function bins_to_index - !=============================================================================== ! ADD_TO_SCORE accumulates a scoring contribution to a specific tally bin and ! specific response function. Note that we don't need to add the square of the diff --git a/src/tally_header.f90 b/src/tally_header.f90 index 92b39c5177..c55c7e52fb 100644 --- a/src/tally_header.f90 +++ b/src/tally_header.f90 @@ -67,27 +67,29 @@ module tally_header ! Tally bin specifications - type(TallyFilter), pointer :: cell_bins(:) => null() - type(TallyFilter), pointer :: surface_bins(:) => null() type(TallyFilter), pointer :: universe_bins(:) => null() type(TallyFilter), pointer :: material_bins(:) => null() - type(TallyFilter), pointer :: mesh_bins(:) => null() - type(TallyFilter), pointer :: bornin_bins(:) => null() + type(TallyFilter), pointer :: cell_bins(:) => null() + type(TallyFilter), pointer :: cellborn_bins(:) => null() + type(TallyFilter), pointer :: surface_bins(:) => null() + type(TallyFilter), pointer :: mesh_bins(:) => null() real(8), allocatable :: energy_in(:) real(8), allocatable :: energy_out(:) ! Number of bins for each filter - - integer :: n_cell_bins = 0 - integer :: n_surface_bins = 0 - integer :: n_universe_bins = 0 - integer :: n_material_bins = 0 - integer :: n_mesh_bins = 0 - integer :: n_bornin_bins = 0 - integer :: n_energy_in = 0 - integer :: n_energy_out = 0 integer :: n_total_bins = 0 + ! The following attributes do not necessarily need to be stored but they + ! greatly simplify logic in many places. n_bins gives the number of bins + ! for each filter type, e.g. n_bins(T_CELL) would be the size of + ! cell_bins. The stride attribute is used for determining the index in the + ! scores array for a bin combination. Since multiple dimensions are mapped + ! onto one dimension in the scores array, the stride attribute gives the + ! stride for a given filter type within the scores array + + integer, allocatable :: n_bins(:) + integer, allocatable :: stride(:) + ! Macroscopic properties to score type(TallyFilter), pointer :: macro_bins(:) => null()