diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 12589fa11c..14709b6eab 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -71,6 +71,7 @@ initialize.o: mpi_routines.o initialize.o: output.o initialize.o: source.o initialize.o: string.o +initialize.o: tally.o initialize.o: timing.o input_xml.o: constants.o diff --git a/src/constants.f90 b/src/constants.f90 index fa7681a0af..dd2cdfbd00 100644 --- a/src/constants.f90 +++ b/src/constants.f90 @@ -160,6 +160,16 @@ module constants MACRO_FISSION = -5, & ! total fission rate MACRO_NU_FISSION = -6 ! total neutron production rate + ! 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 + ! Fission neutron emission (nu) type integer, parameter :: & NU_NONE = 0, & ! No nu values (non-fissionable) diff --git a/src/global.f90 b/src/global.f90 index bc9288aace..2b93b412a3 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -9,7 +9,7 @@ module global use material_header, only: Material use particle_header, only: Particle use source_header, only: ExtSource - use tally_header, only: TallyObject + use tally_header, only: TallyObject, TallyMap use timing, only: Timer #ifdef MPI @@ -56,6 +56,9 @@ module global type(NuclideMicroXS), allocatable :: micro_xs(:) type(MaterialMacroXS) :: material_xs + ! Tally map structure + type(TallyMap), allocatable :: tally_maps(:) + ! unionized energy grid integer :: n_grid ! number of points on unionized grid real(8), allocatable :: e_grid(:) ! energies on unionized grid diff --git a/src/initialize.f90 b/src/initialize.f90 index fc4f55573d..bf32c1307b 100644 --- a/src/initialize.f90 +++ b/src/initialize.f90 @@ -18,6 +18,7 @@ module initialize print_particle, header use source, only: initialize_source use string, only: int_to_str, starts_with, ends_with + use tally, only: create_tally_map, TallyObject use timing, only: timer_start, timer_stop contains @@ -92,6 +93,9 @@ contains call unionized_grid() call original_indices() + ! Create tally map + call create_tally_map() + ! create source particles call initialize_source() @@ -240,11 +244,11 @@ contains integer :: j ! index over surface list integer :: k integer :: index ! index in surfaces/materials array - integer :: universe_num ! user-specified universe number - integer :: surf_num ! user-specified surface number + integer :: uid ! user-specified uid character(MAX_LINE_LEN) :: msg ! output/error message - type(Cell), pointer :: c => null() - type(Lattice), pointer :: l => null() + type(Cell), pointer :: c => null() + type(Lattice), pointer :: l => null() + type(TallyObject), pointer :: t => null() do i = 1, n_cells ! ======================================================================= @@ -252,26 +256,27 @@ contains c => cells(i) do j = 1, c % n_surfaces - surf_num = c % surfaces(j) - if (surf_num < OP_DIFFERENCE) then - index = dict_get_key(surface_dict, abs(surf_num)) - if (index == DICT_NULL) then - surf_num = abs(surf_num) - msg = "Could not find surface " // trim(int_to_str(surf_num)) // & + uid = c % surfaces(j) + if (uid < OP_DIFFERENCE) then + if (dict_has_key(surface_dict, abs(uid))) then + index = dict_get_key(surface_dict, abs(uid)) + c % surfaces(j) = sign(index, uid) + else + msg = "Could not find surface " // trim(int_to_str(abs(uid))) // & & " specified on cell " // trim(int_to_str(c % uid)) call fatal_error(msg) end if - c % surfaces(j) = sign(index, surf_num) end if end do ! ======================================================================= ! ADJUST UNIVERSE INDEX FOR EACH CELL - if (dict_has_key(universe_dict, c % universe)) then - c % universe = dict_get_key(universe_dict, c % universe) + uid = c % universe + if (dict_has_key(universe_dict, uid)) then + c % universe = dict_get_key(universe_dict, uid) else - msg = "Could not find universe " // trim(int_to_str(c % universe)) // & + msg = "Could not find universe " // trim(int_to_str(uid)) // & " specified on cell " // trim(int_to_str(c % uid)) call fatal_error(msg) end if @@ -279,25 +284,26 @@ contains ! ======================================================================= ! ADJUST MATERIAL/FILL POINTERS FOR EACH CELL - if (c % material /= 0) then - if (dict_has_key(material_dict, c % material)) then + uid = c % material + if (uid /= 0) then + if (dict_has_key(material_dict, uid)) then c % type = CELL_NORMAL - c % material = dict_get_key(material_dict, c % material) + c % material = dict_get_key(material_dict, uid) else - msg = "Could not find material " // trim(int_to_str(c % material)) // & + msg = "Could not find material " // trim(int_to_str(uid)) // & " specified on cell " // trim(int_to_str(c % uid)) call fatal_error(msg) end if else - universe_num = c % fill - if (dict_has_key(universe_dict, universe_num)) then + uid = c % fill + if (dict_has_key(universe_dict, uid)) then c % type = CELL_FILL - c % fill = dict_get_key(universe_dict, universe_num) - elseif (dict_has_key(lattice_dict, universe_num)) then + c % fill = dict_get_key(universe_dict, uid) + elseif (dict_has_key(lattice_dict, uid)) then c % type = CELL_LATTICE - c % fill = dict_get_key(lattice_dict, universe_num) + c % fill = dict_get_key(lattice_dict, uid) else - msg = "Specified fill " // trim(int_to_str(c % fill)) // " on cell " // & + msg = "Specified fill " // trim(int_to_str(uid)) // " on cell " // & trim(int_to_str(c % uid)) // " is neither a universe nor a lattice." call fatal_error(msg) end if @@ -311,11 +317,11 @@ contains l => lattices(i) do j = 1, l % n_x do k = 1, l % n_y - universe_num = l % element(j,k) - if (dict_has_key(universe_dict, universe_num)) then - l % element(j,k) = dict_get_key(universe_dict, universe_num) + uid = l % element(j,k) + if (dict_has_key(universe_dict, uid)) then + l % element(j,k) = dict_get_key(universe_dict, uid) else - msg = "Invalid universe number " // trim(int_to_str(universe_num)) & + msg = "Invalid universe number " // trim(int_to_str(uid)) & // " specified on lattice " // trim(int_to_str(l % uid)) call fatal_error(msg) end if @@ -323,6 +329,90 @@ contains end do end do + do i = 1, n_tallies + t => tallies(i) + + ! ======================================================================= + ! ADJUST CELL INDICES FOR EACH TALLY + + if (associated(t % cell_bins)) then + do j = 1, size(t % cell_bins) + uid = t % cell_bins(j) % scalar + if (dict_has_key(cell_dict, uid)) then + t % cell_bins(j) % scalar = dict_get_key(cell_dict, uid) + else + msg = "Could not find cell " // trim(int_to_str(uid)) // & + & " specified on tally " // trim(int_to_str(t % uid)) + call fatal_error(msg) + end if + end do + end if + + ! ======================================================================= + ! ADJUST SURFACE INDICES FOR EACH TALLY + + if (associated(t % surface_bins)) then + do j = 1, size(t % surface_bins) + uid = t % surface_bins(j) % scalar + if (dict_has_key(surface_dict, uid)) then + t % surface_bins(j) % scalar = dict_get_key(surface_dict, uid) + else + msg = "Could not find surface " // trim(int_to_str(uid)) // & + & " specified on tally " // trim(int_to_str(t % uid)) + call fatal_error(msg) + end if + end do + end if + + ! ======================================================================= + ! ADJUST UNIVERSE INDICES FOR EACH TALLY + + if (associated(t % universe_bins)) then + do j = 1, size(t % universe_bins) + uid = t % universe_bins(j) % scalar + if (dict_has_key(universe_dict, uid)) then + t % universe_bins(j) % scalar = dict_get_key(universe_dict, uid) + else + msg = "Could not find universe " // trim(int_to_str(uid)) // & + & " specified on tally " // trim(int_to_str(t % uid)) + call fatal_error(msg) + end if + end do + end if + + ! ======================================================================= + ! ADJUST MATERIAL INDICES FOR EACH TALLY + + if (associated(t % material_bins)) then + do j = 1, size(t % material_bins) + uid = t % material_bins(j) % scalar + if (dict_has_key(material_dict, uid)) then + t % material_bins(j) % scalar = dict_get_key(material_dict, uid) + else + msg = "Could not find material " // trim(int_to_str(uid)) // & + & " specified on tally " // trim(int_to_str(t % uid)) + call fatal_error(msg) + end if + end do + end if + + ! ======================================================================= + ! ADJUST BORNIN INDICES FOR EACH TALLY + + if (associated(t % bornin_bins)) then + do j = 1, size(t % bornin_bins) + uid = t % bornin_bins(j) % scalar + if (dict_has_key(cell_dict, uid)) then + t % bornin_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)) + call fatal_error(msg) + end if + end do + end if + end do + end subroutine adjust_indices !=============================================================================== diff --git a/src/output.f90 b/src/output.f90 index f4bdea7ffc..ca8a00800f 100644 --- a/src/output.f90 +++ b/src/output.f90 @@ -507,15 +507,21 @@ contains type(TallyObject), pointer :: tal integer :: i + integer :: uid character(MAX_LINE_LEN) :: string + type(Cell), pointer :: c => null() + type(Surface), pointer :: s => null() + type(Universe), pointer :: u => null() + type(Material), pointer :: m => null() write(ou,*) 'Tally ' // int_to_str(tal % uid) if (associated(tal % cell_bins)) then string = "" do i = 1, size(tal % cell_bins) - string = trim(string) // ' ' // trim(int_to_str(& - tal % cell_bins(i) % scalar)) + uid = tal % cell_bins(i) % scalar + c => cells(uid) + string = trim(string) // ' ' // trim(int_to_str(c % uid)) end do write(ou, *) ' Cell Bins:' // trim(string) end if @@ -523,17 +529,29 @@ contains if (associated(tal % surface_bins)) then string = "" do i = 1, size(tal % surface_bins) - string = trim(string) // ' ' // trim(int_to_str(& - tal % surface_bins(i) % scalar)) + uid = tal % surface_bins(i) % scalar + s => surfaces(uid) + string = trim(string) // ' ' // trim(int_to_str(s % uid)) end do write(ou, *) ' Surface Bins:' // trim(string) end if + if (associated(tal % universe_bins)) then + string = "" + do i = 1, size(tal % universe_bins) + uid = tal % universe_bins(i) % scalar + u => universes(uid) + string = trim(string) // ' ' // trim(int_to_str(u % uid)) + end do + write(ou, *) ' Material Bins:' // trim(string) + end if + if (associated(tal % material_bins)) then string = "" do i = 1, size(tal % material_bins) - string = trim(string) // ' ' // trim(int_to_str(& - tal % material_bins(i) % scalar)) + uid = tal % material_bins(i) % scalar + m => materials(uid) + string = trim(string) // ' ' // trim(int_to_str(m % uid)) end do write(ou, *) ' Material Bins:' // trim(string) end if @@ -550,8 +568,9 @@ contains if (associated(tal % bornin_bins)) then string = "" do i = 1, size(tal % bornin_bins) - string = trim(string) // ' ' // trim(int_to_str(& - tal % bornin_bins(i) % scalar)) + uid = tal % bornin_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 diff --git a/src/tally.f90 b/src/tally.f90 index 9fd969b610..a499f03cf3 100644 --- a/src/tally.f90 +++ b/src/tally.f90 @@ -82,6 +82,33 @@ contains end subroutine calculate_keff +!=============================================================================== +! CREATE_TALLY_MAP +!=============================================================================== + + subroutine create_tally_map() + + integer :: i + type(TallyObject), pointer :: t => null() + + ! allocate tally map array + allocate(tally_maps(TALLY_MAP_TYPES)) + + ! 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)) + + do i = 1, n_tallies + t => tallies(i) + + end do + + end subroutine create_tally_map + !=============================================================================== ! SCORE_TALLY !=============================================================================== diff --git a/src/tally_header.f90 b/src/tally_header.f90 index 9997aaee1a..0171cc2429 100644 --- a/src/tally_header.f90 +++ b/src/tally_header.f90 @@ -2,6 +2,33 @@ module tally_header implicit none +!=============================================================================== +! TALLYMAPELEMENT gives an index to a tally which is to be scored and the +! corresponding bin for the filter variable +!=============================================================================== + + type TallyMapElement + integer :: index_tally + integer :: index_bin + end type TallyMapElement + +!=============================================================================== +! TALLYMAPITEM +!=============================================================================== + + type TallyMapItem + type(TallyMapElement), pointer :: elements => null() + end type TallyMapItem + +!=============================================================================== +! TALLYMAP contains a list of pairs of indices to tallies and the corresponding +! bin for a given filter. +!=============================================================================== + + type TallyMap + type(TallyMapItem), allocatable :: items(:) + end type TallyMap + !=============================================================================== ! TALLYSCORE provides accumulation of scores in a particular tally bin !===============================================================================