diff --git a/src/constants.F90 b/src/constants.F90 index 58685a5677..648a1d3703 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -224,47 +224,65 @@ module constants ! ============================================================================ ! TALLY-RELATED CONSTANTS - ! Tally macro reactions - integer, parameter :: N_MACRO_TYPES = 15 + ! Tally type integer, parameter :: & - MACRO_FLUX = -1, & ! flux - MACRO_TOTAL = -2, & ! total reaction rate - MACRO_SCATTER = -3, & ! scattering rate - MACRO_NU_SCATTER = -4, & ! scattering production rate - MACRO_SCATTER_1 = -5, & ! first scattering moment - MACRO_SCATTER_2 = -6, & ! second scattering moment - MACRO_SCATTER_3 = -7, & ! third scattering moment - MACRO_N_1N = -8, & ! (n,1n) rate - MACRO_N_2N = -9, & ! (n,2n) rate - MACRO_N_3N = -10, & ! (n,3n) rate - MACRO_N_4N = -11, & ! (n,4n) rate - MACRO_ABSORPTION = -12, & ! absorption rate - MACRO_FISSION = -13, & ! fission rate - MACRO_NU_FISSION = -14, & ! neutron production rate - MACRO_CURRENT = -15 ! partial current + TALLY_VOLUME = 1, & + TALLY_SURFACE_CURRENT = 2 + + ! Tally estimator types + integer, parameter :: & + ESTIMATOR_ANALOG = 1, & + ESTIMATOR_TRACKLENGTH = 2 + + ! Event types for tallies + integer, parameter :: & + EVENT_SURFACE = -2, & + EVENT_LATTICE = -1, & + EVENT_SCATTER = 1, & + EVENT_ABSORB = 2, & + EVENT_FISSION = 3 + + ! Tally score type + integer, parameter :: N_SCORE_TYPES = 15 + integer, parameter :: & + SCORE_FLUX = -1, & ! flux + SCORE_TOTAL = -2, & ! total reaction rate + SCORE_SCATTER = -3, & ! scattering rate + SCORE_NU_SCATTER = -4, & ! scattering production rate + SCORE_SCATTER_1 = -5, & ! first scattering moment + SCORE_SCATTER_2 = -6, & ! second scattering moment + SCORE_SCATTER_3 = -7, & ! third scattering moment + SCORE_N_1N = -8, & ! (n,1n) rate + SCORE_N_2N = -9, & ! (n,2n) rate + SCORE_N_3N = -10, & ! (n,3n) rate + SCORE_N_4N = -11, & ! (n,4n) rate + SCORE_ABSORPTION = -12, & ! absorption rate + SCORE_FISSION = -13, & ! fission rate + SCORE_NU_FISSION = -14, & ! neutron production rate + SCORE_CURRENT = -15 ! partial current ! Tally map bin finding integer, parameter :: NO_BIN_FOUND = -1 ! Tally filter and map types - integer, parameter :: TALLY_TYPES = 8 + integer, parameter :: N_FILTER_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 + FILTER_UNIVERSE = 1, & + FILTER_MATERIAL = 2, & + FILTER_CELL = 3, & + FILTER_CELLBORN = 4, & + FILTER_SURFACE = 5, & + FILTER_MESH = 6, & + FILTER_ENERGYIN = 7, & + FILTER_ENERGYOUT = 8 ! Filter types for surface current tallies integer, parameter :: & - TS_MESH_X = 1, & - TS_MESH_Y = 2, & - TS_MESH_Z = 3, & - TS_ENERGYIN = 4, & - TS_SURFACE = 5 + SURF_FILTER_MESH_X = 1, & + SURF_FILTER_MESH_Y = 2, & + SURF_FILTER_MESH_Z = 3, & + SURF_FILTER_ENERGYIN = 4, & + SURF_FILTER_SURFACE = 5 ! Tally surface current directions integer, parameter :: & diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index e4f077d9b1..481765757c 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -110,6 +110,9 @@ contains call hdf5_write_geometry() call hdf5_write_materials() + if (n_tallies > 0) then + call hdf5_write_tallies() + end if end subroutine hdf5_write_summary @@ -422,6 +425,197 @@ contains end subroutine hdf5_write_materials +!=============================================================================== +! HDF5_WRITE_TALLIES +!=============================================================================== + + subroutine hdf5_write_tallies() + + integer :: i, j + integer(SIZE_T) :: num_elements + integer(HSIZE_T) :: dims(1) + integer(HSIZE_T) :: coord(1,1) + integer(HSIZE_T) :: size_string = 12 + integer(HID_T) :: string_type + integer(HID_T) :: tallies_group + integer(HID_T) :: temp_group + integer(HID_T) :: dspace + integer(HID_T) :: subspace + integer(HID_T) :: dset + character(12) :: string + type(TallyObject), pointer :: t => null() + type(StructuredMesh), pointer :: m => null() + + ! Create group for tallies + call h5gcreate_f(hdf5_output_file, "/tallies", tallies_group, hdf5_err) + + ! Use H5LT interface to write number of tallies + call hdf5_make_integer(tallies_group, "n_tallies", n_tallies) + + ! Write information on each material + do i = 1, n_tallies + t => tallies(i) + + ! Create group for i-th universe + call h5gcreate_f(tallies_group, "tally " // trim(to_str(t % id)), & + temp_group, hdf5_err) + + ! ======================================================================= + ! WRITE INFORMATION ON TALLY FILTERS + + ! Write filters + call hdf5_make_integer(temp_group, "n_filters", t % n_filters) + if (t % n_filters > 0) then + dims(1) = t % n_filters + call h5ltmake_dataset_int_f(temp_group, "filters", 1, & + dims, t % filters, hdf5_err) + end if + + ! Write universe_bins if present + if (t % n_filter_bins(FILTER_UNIVERSE) > 0) then + dims(1) = t % n_filter_bins(FILTER_UNIVERSE) + call h5ltmake_dataset_int_f(temp_group, "universe_bins", 1, & + dims, t % universe_bins(:) % scalar, hdf5_err) + end if + + ! Write material_bins if present + if (t % n_filter_bins(FILTER_MATERIAL) > 0) then + dims(1) = t % n_filter_bins(FILTER_MATERIAL) + call h5ltmake_dataset_int_f(temp_group, "material_bins", 1, & + dims, t % material_bins(:) % scalar, hdf5_err) + end if + + ! Write cell_bins if present + if (t % n_filter_bins(FILTER_CELL) > 0) then + dims(1) = t % n_filter_bins(FILTER_CELL) + call h5ltmake_dataset_int_f(temp_group, "cell_bins", 1, & + dims, t % cell_bins(:) % scalar, hdf5_err) + end if + + ! Write cellborn_bins if present + if (t % n_filter_bins(FILTER_CELLBORN) > 0) then + dims(1) = t % n_filter_bins(FILTER_CELLBORN) + call h5ltmake_dataset_int_f(temp_group, "cellborn_bins", 1, & + dims, t % cellborn_bins(:) % scalar, hdf5_err) + end if + + ! Write surface_bins if present + if (t % n_filter_bins(FILTER_SURFACE) > 0) then + dims(1) = t % n_filter_bins(FILTER_SURFACE) + call h5ltmake_dataset_int_f(temp_group, "surface_bins", 1, & + dims, t % surface_bins(:) % scalar, hdf5_err) + end if + + ! Write incoming energy filter + if (t % n_filter_bins(FILTER_ENERGYIN) > 0) then + dims(1) = t % n_filter_bins(FILTER_ENERGYIN) + call h5ltmake_dataset_double_f(temp_group, "energy_in", 1, & + dims, t % energy_in, hdf5_err) + end if + + ! Write outgoing energy filter + if (t % n_filter_bins(FILTER_ENERGYOUT) > 0) then + dims(1) = t % n_filter_bins(FILTER_ENERGYOUT) + call h5ltmake_dataset_double_f(temp_group, "energy_out", 1, & + dims, t % energy_out, hdf5_err) + end if + + ! Write mesh information + if (t % n_filter_bins(FILTER_MESH) > 0) then + m => meshes(t % mesh) + dims(1) = m % n_dimension + + ! Write mesh dimensions + call h5ltmake_dataset_int_f(temp_group, "mesh_dimensions", 1, & + dims, m % dimension, hdf5_err) + + ! Write mesh lower-left corner + call h5ltmake_dataset_double_f(temp_group, "mesh_lower_left", 1, & + dims, m % origin, hdf5_err) + + ! Write mesh element width + call h5ltmake_dataset_double_f(temp_group, "mesh_element_width", 1, & + dims, m % width, hdf5_err) + end if + + ! ======================================================================= + ! WRITE INFORMATION ON TALLY SCORES + + ! Create string type of length 12 + call h5tcopy_f(H5T_C_S1, string_type, hdf5_err) + call h5tset_size_f(string_type, size_string, hdf5_err) + + ! Create dataspace and dataset for writing nuclides + dims(1) = t % n_score_bins + call h5screate_simple_f(1, dims, dspace, hdf5_err) + call h5dcreate_f(temp_group, "scores", string_type, dspace, & + dset, hdf5_err) + + ! Create subspace for writing one element + dims(1) = 1 + call h5screate_simple_f(1, dims, subspace, hdf5_err) + + ! Write list of nuclides + num_elements = 1 + do j = 1, t % n_score_bins + ! Determine string for this scoring bin + select case (t % score_bins(j) % scalar) + case (SCORE_FLUX) + string = 'flux' + case (SCORE_TOTAL) + string = 'total' + case (SCORE_SCATTER) + string = 'scatter' + case (SCORE_NU_SCATTER) + string = 'nu-scatter' + case (SCORE_SCATTER_1) + string = '1st moment' + case (SCORE_SCATTER_2) + string = '2nd moment' + case (SCORE_SCATTER_3) + string = '3rd moment' + case (SCORE_N_1N) + string = '(n,1n)' + case (SCORE_N_2N) + string = '(n,2n)' + case (SCORE_N_3N) + string = '(n,3n)' + case (SCORE_N_4N) + string = '(n,4n)' + case (SCORE_ABSORPTION) + string = 'absorption' + case (SCORE_FISSION) + string = 'fission' + case (SCORE_NU_FISSION) + string = 'nu-fission' + case (SCORE_CURRENT) + string = 'current' + end select + + ! Select block within array + coord(1,1) = j + call h5sselect_elements_f(dspace, H5S_SELECT_SET_F, 1, & + num_elements, coord, hdf5_err) + + ! Write data to partial array + call h5dwrite_f(dset, string_type, string, dims, hdf5_err, & + subspace, dspace) + end do + + ! Close dataspace and dataset for nuclides + call h5dclose_f(dset, hdf5_err) + call h5sclose_f(subspace, hdf5_err) + call h5sclose_f(dspace, hdf5_err) + + ! Close group for i-th material + call h5gclose_f(temp_group, hdf5_err) + end do + + ! Close tallies group + call h5gclose_f(tallies_group, hdf5_err) + + end subroutine hdf5_write_tallies + !=============================================================================== ! HDF5_WRITE_TIMING !=============================================================================== diff --git a/src/initialize.F90 b/src/initialize.F90 index f008f216b7..e8dda1d3d1 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -459,8 +459,8 @@ contains ! ======================================================================= ! ADJUST CELL INDICES FOR EACH TALLY - if (t % n_bins(T_CELL) > 0) then - do j = 1, t % n_bins(T_CELL) + if (t % n_filter_bins(FILTER_CELL) > 0) then + do j = 1, t % n_filter_bins(FILTER_CELL) id = t % cell_bins(j) % scalar if (dict_has_key(cell_dict, id)) then t % cell_bins(j) % scalar = dict_get_key(cell_dict, id) @@ -475,8 +475,8 @@ contains ! ======================================================================= ! ADJUST SURFACE INDICES FOR EACH TALLY - if (t % n_bins(T_SURFACE) > 0) then - do j = 1, t % n_bins(T_SURFACE) + if (t % n_filter_bins(FILTER_SURFACE) > 0) then + do j = 1, t % n_filter_bins(FILTER_SURFACE) id = t % surface_bins(j) % scalar if (dict_has_key(surface_dict, id)) then t % surface_bins(j) % scalar = dict_get_key(surface_dict, id) @@ -491,8 +491,8 @@ contains ! ======================================================================= ! ADJUST UNIVERSE INDICES FOR EACH TALLY - if (t % n_bins(T_UNIVERSE) > 0) then - do j = 1, t % n_bins(T_UNIVERSE) + if (t % n_filter_bins(FILTER_UNIVERSE) > 0) then + do j = 1, t % n_filter_bins(FILTER_UNIVERSE) id = t % universe_bins(j) % scalar if (dict_has_key(universe_dict, id)) then t % universe_bins(j) % scalar = dict_get_key(universe_dict, id) @@ -507,8 +507,8 @@ contains ! ======================================================================= ! ADJUST MATERIAL INDICES FOR EACH TALLY - if (t % n_bins(T_MATERIAL) > 0) then - do j = 1, t % n_bins(T_MATERIAL) + if (t % n_filter_bins(FILTER_MATERIAL) > 0) then + do j = 1, t % n_filter_bins(FILTER_MATERIAL) id = t % material_bins(j) % scalar if (dict_has_key(material_dict, id)) then t % material_bins(j) % scalar = dict_get_key(material_dict, id) @@ -523,8 +523,8 @@ contains ! ======================================================================= ! ADJUST CELLBORN INDICES FOR EACH TALLY - if (t % n_bins(T_CELLBORN) > 0) then - do j = 1, t % n_bins(T_CELLBORN) + if (t % n_filter_bins(FILTER_CELLBORN) > 0) then + do j = 1, t % n_filter_bins(FILTER_CELLBORN) id = t % cellborn_bins(j) % scalar if (dict_has_key(cell_dict, id)) then t % cellborn_bins(j) % scalar = dict_get_key(cell_dict, id) @@ -539,7 +539,7 @@ contains ! ======================================================================= ! ADJUST MESH INDICES FOR EACH TALLY - if (t % n_bins(T_MESH) > 0) then + if (t % n_filter_bins(FILTER_MESH) > 0) then id = t % mesh if (dict_has_key(mesh_dict, id)) then t % mesh = dict_get_key(mesh_dict, id) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index b9babeb557..f0c1452596 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -617,6 +617,8 @@ contains integer :: i_mesh ! index in meshes array integer :: n ! size of arrays in mesh specification integer :: n_words ! number of words read + integer :: n_filters ! number of filters + integer :: filters(N_FILTER_TYPES) ! temporary list of filters logical :: file_exists ! does tallies.xml file exist? character(MAX_LINE_LEN) :: filename character(MAX_WORD_LEN) :: word @@ -730,13 +732,28 @@ contains 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)) + allocate(t % n_filter_bins(N_FILTER_TYPES)) + allocate(t % stride(N_FILTER_TYPES)) ! Initialize number of bins and stride - t % n_bins = 0 + t % n_filter_bins = 0 t % stride = 0 + ! Initialize filters + n_filters = 0 + filters = 0 + + ! Set tally type to volume by default + t % type = TALLY_VOLUME + + ! It's desirable to use a track-length esimator for tallies since + ! generally more events will score to the tally, reducing the + ! variance. However, for tallies that require information on + ! post-collision parameters (e.g. tally with an energyout filter) the + ! analog esimator must be used. + + t % estimator = ESTIMATOR_TRACKLENGTH + ! Copy material id t % id = tally_(i) % id @@ -757,7 +774,10 @@ contains do j = 1, n_words t % cell_bins(j) % scalar = int(str_to_int(words(j)),4) end do - t % n_bins(T_CELL) = n_words + t % n_filter_bins(FILTER_CELL) = n_words + + n_filters = n_filters + 1 + filters(n_filters) = FILTER_CELL end if ! Read surface filter bins @@ -767,7 +787,10 @@ contains do j = 1, n_words t % surface_bins(j) % scalar = int(str_to_int(words(j)),4) end do - t % n_bins(T_SURFACE) = n_words + t % n_filter_bins(FILTER_SURFACE) = n_words + + n_filters = n_filters + 1 + filters(n_filters) = FILTER_SURFACE end if ! Read universe filter bins @@ -777,7 +800,10 @@ contains do j = 1, n_words t % universe_bins(j) % scalar = int(str_to_int(words(j)),4) end do - t % n_bins(T_UNIVERSE) = n_words + t % n_filter_bins(FILTER_UNIVERSE) = n_words + + n_filters = n_filters + 1 + filters(n_filters) = FILTER_UNIVERSE end if ! Read material filter bins @@ -787,7 +813,10 @@ contains do j = 1, n_words t % material_bins(j) % scalar = int(str_to_int(words(j)),4) end do - t % n_bins(T_MATERIAL) = n_words + t % n_filter_bins(FILTER_MATERIAL) = n_words + + n_filters = n_filters + 1 + filters(n_filters) = FILTER_MATERIAL end if ! Read mesh filter bins @@ -804,7 +833,10 @@ contains call fatal_error() end if - t % n_bins(T_MESH) = t % n_bins(T_MESH) + product(m % dimension) + t % n_filter_bins(FILTER_MESH) = t % n_filter_bins(FILTER_MESH) + product(m % dimension) + + n_filters = n_filters + 1 + filters(n_filters) = FILTER_MESH end if ! Read birth region filter bins @@ -814,7 +846,10 @@ contains do j = 1, n_words t % cellborn_bins(j) % scalar = int(str_to_int(words(j)),4) end do - t % n_bins(T_CELLBORN) = n_words + t % n_filter_bins(FILTER_CELLBORN) = n_words + + n_filters = n_filters + 1 + filters(n_filters) = FILTER_CELLBORN end if ! Read incoming energy filter bins @@ -824,7 +859,10 @@ contains do j = 1, n_words t % energy_in(j) = str_to_real(words(j)) end do - t % n_bins(T_ENERGYIN) = n_words - 1 + t % n_filter_bins(FILTER_ENERGYIN) = n_words - 1 + + n_filters = n_filters + 1 + filters(n_filters) = FILTER_ENERGYIN end if ! Read outgoing energy filter bins @@ -834,88 +872,117 @@ contains do j = 1, n_words t % energy_out(j) = str_to_real(words(j)) end do - t % n_bins(T_ENERGYOUT) = n_words - 1 + t % n_filter_bins(FILTER_ENERGYOUT) = n_words - 1 + + ! Set tally estimator to analog + t % estimator = ESTIMATOR_ANALOG + + n_filters = n_filters + 1 + filters(n_filters) = FILTER_ENERGYOUT end if + ! Allocate and set filters + t % n_filters = n_filters + allocate(t % filters(n_filters)) + t % filters = filters(1:n_filters) + ! Read macro reactions if (len_trim(tally_(i) % macros) > 0) then call split_string(tally_(i) % macros, words, n_words) - allocate(t % macro_bins(n_words)) + allocate(t % score_bins(n_words)) do j = 1, n_words word = words(j) call lower_case(word) select case (trim(word)) case ('flux') - t % macro_bins(j) % scalar = MACRO_FLUX - if (t % n_bins(T_ENERGYOUT) > 0) then + t % score_bins(j) % scalar = SCORE_FLUX + if (t % n_filter_bins(FILTER_ENERGYOUT) > 0) then message = "Cannot tally flux with an outgoing energy filter." call fatal_error() end if case ('total') - t % macro_bins(j) % scalar = MACRO_TOTAL - if (t % n_bins(T_ENERGYOUT) > 0) then + t % score_bins(j) % scalar = SCORE_TOTAL + if (t % n_filter_bins(FILTER_ENERGYOUT) > 0) then message = "Cannot tally total reaction rate with an & &outgoing energy filter." call fatal_error() end if case ('scatter') - t % macro_bins(j) % scalar = MACRO_SCATTER + t % score_bins(j) % scalar = SCORE_SCATTER case ('nu-scatter') - t % macro_bins(j) % scalar = MACRO_NU_SCATTER + t % score_bins(j) % scalar = SCORE_NU_SCATTER + + ! Set tally estimator to analog + t % estimator = ESTIMATOR_ANALOG case ('scatter-1') - t % macro_bins(j) % scalar = MACRO_SCATTER_1 + t % score_bins(j) % scalar = SCORE_SCATTER_1 + + ! Set tally estimator to analog + t % estimator = ESTIMATOR_ANALOG case ('scatter-2') - t % macro_bins(j) % scalar = MACRO_SCATTER_2 + t % score_bins(j) % scalar = SCORE_SCATTER_2 + + ! Set tally estimator to analog + t % estimator = ESTIMATOR_ANALOG case ('scatter-3') - t % macro_bins(j) % scalar = MACRO_SCATTER_3 + t % score_bins(j) % scalar = SCORE_SCATTER_3 + + ! Set tally estimator to analog + t % estimator = ESTIMATOR_ANALOG case ('n1n') - t % macro_bins(j) % scalar = MACRO_N_1N + t % score_bins(j) % scalar = SCORE_N_1N + + ! Set tally estimator to analog + t % estimator = ESTIMATOR_ANALOG case ('n2n') - t % macro_bins(j) % scalar = MACRO_N_2N + t % score_bins(j) % scalar = SCORE_N_2N + + ! Set tally estimator to analog + t % estimator = ESTIMATOR_ANALOG case ('n3n') - t % macro_bins(j) % scalar = MACRO_N_3N + t % score_bins(j) % scalar = SCORE_N_3N + + ! Set tally estimator to analog + t % estimator = ESTIMATOR_ANALOG case ('n4n') - t % macro_bins(j) % scalar = MACRO_N_4N + t % score_bins(j) % scalar = SCORE_N_4N + + ! Set tally estimator to analog + t % estimator = ESTIMATOR_ANALOG case ('absorption') - t % macro_bins(j) % scalar = MACRO_ABSORPTION - if (t % n_bins(T_ENERGYOUT) > 0) then + t % score_bins(j) % scalar = SCORE_ABSORPTION + if (t % n_filter_bins(FILTER_ENERGYOUT) > 0) then message = "Cannot tally absorption rate with an outgoing & &energy filter." call fatal_error() end if case ('fission') - t % macro_bins(j) % scalar = MACRO_FISSION - if (t % n_bins(T_ENERGYOUT) > 0) then + t % score_bins(j) % scalar = SCORE_FISSION + if (t % n_filter_bins(FILTER_ENERGYOUT) > 0) then message = "Cannot tally fission rate with an outgoing & &energy filter." call fatal_error() end if case ('nu-fission') - t % macro_bins(j) % scalar = MACRO_NU_FISSION + t % score_bins(j) % scalar = SCORE_NU_FISSION case ('current') - t % macro_bins(j) % scalar = MACRO_CURRENT - t % surface_current = .true. + t % score_bins(j) % scalar = SCORE_CURRENT + t % type = TALLY_SURFACE_CURRENT ! Check to make sure that current is the only desired response ! for this tally if (n_words > 1) then - message = "Cannot tally other macro reactions in the same & - &tally as surface currents. Separate other macro & - &reactions into a distinct tally." + message = "Cannot tally other scoring functions in the same & + &tally as surface currents. Separate other scoring & + &functions into a distinct tally." call fatal_error() end if - ! Check to make sure that only the mesh filter was specified -!!$ if (t % mesh == 0 .or. t % n_bins(T_MESH) /= & -!!$ product(t % n_bins, t % n_bins > 0)) then -!!$ message = "Surface currents must be used with a mesh filter only." -!!$ call fatal_error() -!!$ end if - ! Since the number of bins for the mesh filter was already set ! assuming it was a flux tally, we need to adjust the number of ! bins - t % n_bins(T_MESH) = t % n_bins(T_MESH) - product(m % dimension) + t % n_filter_bins(FILTER_MESH) = t % n_filter_bins(FILTER_MESH) & + - product(m % dimension) ! Get pointer to mesh id = t % mesh @@ -925,19 +992,19 @@ contains ! We need to increase the dimension by one since we also need ! currents coming into and out of the boundary mesh cells. if (size(m % dimension) == 2) then - t % n_bins(T_MESH) = t % n_bins(T_MESH) + & - product(m % dimension + 1) * 4 + t % n_filter_bins(FILTER_MESH) = t % n_filter_bins(FILTER_MESH) & + + product(m % dimension + 1) * 4 elseif (size(m % dimension) == 3) then - t % n_bins(T_MESH) = t % n_bins(T_MESH) + & - product(m % dimension + 1) * 6 + t % n_filter_bins(FILTER_MESH) = t % n_filter_bins(FILTER_MESH) & + + product(m % dimension + 1) * 6 end if case default - message = "Unknown macro reaction: " // trim(words(j)) + message = "Unknown scoring function: " // trim(words(j)) call fatal_error() end select end do - t % n_macro_bins = n_words + t % n_score_bins = n_words end if end do diff --git a/src/intercycle.F90 b/src/intercycle.F90 index 50301927ce..44c2cbf1f9 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -299,7 +299,7 @@ contains t => tallies(i) n = t % n_total_bins - m = t % n_macro_bins + m = t % n_score_bins n_bins = n*m allocate(tally_temp(n,m)) diff --git a/src/output.F90 b/src/output.F90 index e6dcf1f952..fb390f005d 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -575,9 +575,9 @@ contains write(unit_,*) 'Tally ' // to_str(t % id) - if (t % n_bins(T_CELL) > 0) then + if (t % n_filter_bins(FILTER_CELL) > 0) then string = "" - do i = 1, t % n_bins(T_CELL) + do i = 1, t % n_filter_bins(FILTER_CELL) id = t % cell_bins(i) % scalar c => cells(id) string = trim(string) // ' ' // trim(to_str(c % id)) @@ -585,9 +585,9 @@ contains write(unit_, *) ' Cell Bins:' // trim(string) end if - if (t % n_bins(T_SURFACE) > 0) then + if (t % n_filter_bins(FILTER_SURFACE) > 0) then string = "" - do i = 1, t % n_bins(T_SURFACE) + do i = 1, t % n_filter_bins(FILTER_SURFACE) id = t % surface_bins(i) % scalar s => surfaces(id) string = trim(string) // ' ' // trim(to_str(s % id)) @@ -595,9 +595,9 @@ contains write(unit_, *) ' Surface Bins:' // trim(string) end if - if (t % n_bins(T_UNIVERSE) > 0) then + if (t % n_filter_bins(FILTER_UNIVERSE) > 0) then string = "" - do i = 1, t % n_bins(T_UNIVERSE) + do i = 1, t % n_filter_bins(FILTER_UNIVERSE) id = t % universe_bins(i) % scalar u => universes(id) string = trim(string) // ' ' // trim(to_str(u % id)) @@ -605,9 +605,9 @@ contains write(unit_, *) ' Material Bins:' // trim(string) end if - if (t % n_bins(T_MATERIAL) > 0) then + if (t % n_filter_bins(FILTER_MATERIAL) > 0) then string = "" - do i = 1, t % n_bins(T_MATERIAL) + do i = 1, t % n_filter_bins(FILTER_MATERIAL) id = t % material_bins(i) % scalar m => materials(id) string = trim(string) // ' ' // trim(to_str(m % id)) @@ -615,7 +615,7 @@ contains write(unit_, *) ' Material Bins:' // trim(string) end if - if (t % n_bins(T_MESH) > 0) then + if (t % n_filter_bins(FILTER_MESH) > 0) then string = "" id = t % mesh sm => meshes(id) @@ -626,9 +626,9 @@ contains write(unit_, *) ' Mesh Bins:' // trim(string) end if - if (t % n_bins(T_CELLBORN) > 0) then + if (t % n_filter_bins(FILTER_CELLBORN) > 0) then string = "" - do i = 1, t % n_bins(T_CELLBORN) + do i = 1, t % n_filter_bins(FILTER_CELLBORN) id = t % cellborn_bins(i) % scalar c => cells(id) string = trim(string) // ' ' // trim(to_str(c % id)) @@ -636,39 +636,39 @@ contains write(unit_, *) ' Birth Region Bins:' // trim(string) end if - if (t % n_bins(T_ENERGYIN) > 0) then + if (t % n_filter_bins(FILTER_ENERGYIN) > 0) then string = "" - do i = 1, t % n_bins(T_ENERGYIN) + 1 + do i = 1, t % n_filter_bins(FILTER_ENERGYIN) + 1 string = trim(string) // ' ' // trim(to_str(& t % energy_in(i))) end do write(unit_,*) ' Incoming Energy Bins:' // trim(string) end if - if (t % n_bins(T_ENERGYOUT) > 0) then + if (t % n_filter_bins(FILTER_ENERGYOUT) > 0) then string = "" - do i = 1, t % n_bins(T_ENERGYOUT) + 1 + do i = 1, t % n_filter_bins(FILTER_ENERGYOUT) + 1 string = trim(string) // ' ' // trim(to_str(& t % energy_out(i))) end do write(unit_,*) ' Outgoing Energy Bins:' // trim(string) end if - if (t % n_macro_bins > 0) then + if (t % n_score_bins > 0) then string = "" - do i = 1, t % n_macro_bins - select case (t % macro_bins(i) % scalar) - case (MACRO_FLUX) + do i = 1, t % n_score_bins + select case (t % score_bins(i) % scalar) + case (SCORE_FLUX) string = trim(string) // ' flux' - case (MACRO_TOTAL) + case (SCORE_TOTAL) string = trim(string) // ' total' - case (MACRO_SCATTER) + case (SCORE_SCATTER) string = trim(string) // ' scatter' - case (MACRO_ABSORPTION) + case (SCORE_ABSORPTION) string = trim(string) // ' absorption' - case (MACRO_FISSION) + case (SCORE_FISSION) string = trim(string) // ' fission' - case (MACRO_NU_FISSION) + case (SCORE_NU_FISSION) string = trim(string) // ' nu-fission' end select end do diff --git a/src/particle_header.F90 b/src/particle_header.F90 index 8f81fbeb1f..076f6d1bcb 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -1,6 +1,6 @@ module particle_header - use constants, only: NEUTRON, ONE, NONE + use constants, only: NEUTRON, ONE, NONE, ZERO use geometry_header, only: BASE_UNIVERSE implicit none @@ -23,9 +23,8 @@ module particle_header real(8) :: xyz(3) real(8) :: uvw(3) - ! Pointers to next (lower) and previous (higher) universe + ! Pointer to next (more local) set of coordinates type(LocalCoord), pointer :: next => null() - type(LocalCoord), pointer :: prev => null() end type LocalCoord !=============================================================================== @@ -53,6 +52,9 @@ module particle_header real(8) :: last_wgt ! last particle weight real(8) :: last_E ! last energy + ! What event last took place + integer :: event + ! Post-collision physical data integer :: n_bank ! number of fission sites banked diff --git a/src/physics.F90 b/src/physics.F90 index 2d07d2f4e9..9e6bffef31 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -16,7 +16,8 @@ module physics use random_lcg, only: prn use search, only: binary_search use string, only: to_str - use tally, only: score_tally, score_surface_current + use tally, only: score_analog_tally, score_tracklength_tally, & + score_surface_current implicit none @@ -45,7 +46,6 @@ contains ! Particle couldn't be located if (.not. found_cell) then message = "Could not locate particle " // trim(to_str(p % id)) - print *, p % coord0 % xyz call fatal_error() end if @@ -67,7 +67,6 @@ contains ! Initialize number of events to zero n_event = 0 - ! find energy index, interpolation factor do while (p % alive) ! Calculate microscopic and macroscopic cross sections -- note: if the @@ -96,15 +95,22 @@ contains coord => coord % next end do + ! Score track-length tallies + if (tallies_on) call score_tracklength_tally(p, distance) + if (d_collision > d_boundary) then last_cell = p % coord % cell p % coord % cell = NONE if (lattice_crossed /= NONE) then + ! Particle crosses lattice boundary p % surface = NONE call cross_lattice(p, lattice_crossed) + p % event = EVENT_LATTICE else + ! Particle crosses surface p % surface = surface_crossed call cross_surface(p, last_cell) + p % event = EVENT_SURFACE end if else ! collision @@ -149,10 +155,6 @@ contains type(Particle), pointer :: p - integer :: MT ! ENDF reaction number - logical :: scattered ! was this a scattering reaction? - logical :: fissioned ! was this a fission reaction? - ! Store pre-collision particle properties p % last_wgt = p % wgt p % last_E = p % E @@ -167,13 +169,13 @@ contains if (tallies_on) call score_surface_current(p) ! Sample nuclide/reaction for the material the particle is in - call sample_reaction(p, MT) + call sample_reaction(p) ! Display information about collision if (verbosity >= 10 .or. trace) then - message = " " // trim(reaction_name(MT)) // ". Energy = " // & - trim(to_str(p % E * 1e6_8)) // " eV." - call write_message() +!!$ message = " " // trim(reaction_name(MT)) // ". Energy = " // & +!!$ trim(to_str(p % E * 1e6_8)) // " eV." +!!$ call write_message() end if ! check for very low energy @@ -183,21 +185,12 @@ contains call warning() end if - ! Check if particle scattered or fissioned - if (survival_biasing) then - fissioned = .false. - scattered = .true. - else - fissioned = is_fission(MT) - scattered = is_scatter(MT) - end if - ! Score collision estimator tallies for any macro tallies -- this is done ! after a collision has occurred rather than before because we need ! information on the outgoing energy for any tallies with an outgoing energy ! filter - if (tallies_on) call score_tally(p, scattered, fissioned) + if (tallies_on) call score_analog_tally(p) ! Reset number of particles banked during collision p % n_bank = 0 @@ -212,10 +205,9 @@ contains ! disappearance are treated implicitly. !=============================================================================== - subroutine sample_reaction(p, MT) + subroutine sample_reaction(p) type(Particle), pointer :: p - integer, intent(out) :: MT ! ENDF MT number of reaction that occured integer :: i ! index over nuclides in a material integer :: index_nuclide ! index in nuclides array @@ -282,7 +274,7 @@ contains ! See if disappearance reaction happens if (prob > cutoff) then p % alive = .false. - MT = N_DISAPPEAR + p % event = EVENT_ABSORB return end if end if @@ -338,16 +330,24 @@ contains if (micro_xs(index_nuclide) % use_ptable) then + ! In the case that the particle is in the unresolved resonance + ! region and probability tables are being used, we should only + ! check the first fission reaction + prob = prob + micro_xs(index_nuclide) % fission if (prob > cutoff) then rxn => nuc % reactions(nuc % index_fission(1)) call create_fission_sites(p, index_nuclide, rxn, .true.) p % alive = .false. - MT = rxn % MT + p % event = EVENT_FISSION return end if + else + ! With no probability tables, we need to loop through each fission + ! reaction type + do i = 1, nuc % n_fission rxn => nuc % reactions(nuc % index_fission(i)) @@ -366,7 +366,7 @@ contains if (prob > cutoff) then call create_fission_sites(p, index_nuclide, rxn, .true.) p % alive = .false. - MT = rxn % MT + p % event = EVENT_FISSION return end if end do @@ -419,8 +419,6 @@ contains end if - ! Set MT to be returned - MT = 2 else ! ======================================================================= ! INELASTIC SCATTERING @@ -461,10 +459,12 @@ contains ! Perform collision physics for inelastics scattering call inelastic_scatter(p, nuc, rxn) - ! Set MT to be returned - MT = rxn % MT end if + ! If we made it this far, it means that a scattering reaction took place + ! since the absorption and fission blocks had return statements. + p % event = EVENT_SCATTER + end subroutine sample_reaction !=============================================================================== diff --git a/src/tally.F90 b/src/tally.F90 index a1d66674ec..67393d46ac 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -41,17 +41,17 @@ contains ! 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 - 3)) + allocate(tally_maps(N_FILTER_TYPES - 3)) ! allocate list of items for each different filter type - allocate(tally_maps(T_UNIVERSE) % items(n_universes)) - allocate(tally_maps(T_MATERIAL) % items(n_materials)) - allocate(tally_maps(T_CELL) % items(n_cells)) - allocate(tally_maps(T_CELLBORN) % items(n_cells)) - allocate(tally_maps(T_SURFACE) % items(n_surfaces)) + allocate(tally_maps(FILTER_UNIVERSE) % items(n_universes)) + allocate(tally_maps(FILTER_MATERIAL) % items(n_materials)) + allocate(tally_maps(FILTER_CELL) % items(n_cells)) + allocate(tally_maps(FILTER_CELLBORN) % items(n_cells)) + allocate(tally_maps(FILTER_SURFACE) % items(n_surfaces)) ! Allocate and initialize tally map positioning for finding bins - allocate(position(TALLY_TYPES - 3)) + allocate(position(N_FILTER_TYPES - 3)) position = 0 do i = 1, n_tallies @@ -61,10 +61,10 @@ contains filter_bins = 1 ! handle surface current tallies separately - if (t % surface_current) then + if (t % type == TALLY_SURFACE_CURRENT) then m => meshes(t % mesh) - t % stride(TS_SURFACE) = filter_bins + t % stride(SURF_FILTER_SURFACE) = filter_bins ! Set stride for surface/direction if (m % n_dimension == 2) then filter_bins = filter_bins * 4 @@ -73,27 +73,27 @@ contains end if ! Add filter for incoming energy - n = t % n_bins(T_ENERGYIN) - t % stride(TS_ENERGYIN) = filter_bins + n = t % n_filter_bins(FILTER_ENERGYIN) + t % stride(SURF_FILTER_ENERGYIN) = filter_bins if (n > 0) then filter_bins = filter_bins * n end if ! account for z direction - t % stride(TS_MESH_Z) = filter_bins + t % stride(SURF_FILTER_MESH_Z) = filter_bins filter_bins = filter_bins * (m % dimension(3) + 1) ! account for y direction - t % stride(TS_MESH_Y) = filter_bins + t % stride(SURF_FILTER_MESH_Y) = filter_bins filter_bins = filter_bins * (m % dimension(2) + 1) ! account for z direction - t % stride(TS_MESH_X) = filter_bins + t % stride(SURF_FILTER_MESH_X) = filter_bins filter_bins = filter_bins * (m % dimension(1) + 1) - ! Finally add scoring bins for the macro tallies and allocate scores - ! for tally - score_bins = t % n_macro_bins + ! Finally add scoring bins for the tallies and allocate the main + ! scores array for tally + score_bins = t % n_score_bins t % n_total_bins = filter_bins allocate(t % scores(filter_bins, score_bins)) @@ -104,86 +104,86 @@ contains ! 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 + n = t % n_filter_bins(FILTER_ENERGYOUT) + t % stride(FILTER_ENERGYOUT) = filter_bins if (n > 0) then filter_bins = filter_bins * n end if - n = t % n_bins(T_ENERGYIN) - t % stride(T_ENERGYIN) = filter_bins + n = t % n_filter_bins(FILTER_ENERGYIN) + t % stride(FILTER_ENERGYIN) = filter_bins if (n > 0) then filter_bins = filter_bins * n end if ! Add map elements for mesh bins - n = t % n_bins(T_MESH) - t % stride(T_MESH) = filter_bins + n = t % n_filter_bins(FILTER_MESH) + t % stride(FILTER_MESH) = filter_bins if (n > 0) then filter_bins = filter_bins * n end if ! Add map elements for surface bins - n = t % n_bins(T_SURFACE) - t % stride(T_SURFACE) = filter_bins + n = t % n_filter_bins(FILTER_SURFACE) + t % stride(FILTER_SURFACE) = filter_bins if (n > 0) then do j = 1, n i_item = t % surface_bins(j) % scalar - call add_map_element(tally_maps(T_SURFACE) % items(i_item), i, j) + call add_map_element(tally_maps(FILTER_SURFACE) % items(i_item), i, j) end do filter_bins = filter_bins * n end if ! Add map elements for cellborn bins - n = t % n_bins(T_CELLBORN) - t % stride(T_CELLBORN) = filter_bins + n = t % n_filter_bins(FILTER_CELLBORN) + t % stride(FILTER_CELLBORN) = filter_bins if (n > 0) then do j = 1, n i_item = t % cellborn_bins(j) % scalar - call add_map_element(tally_maps(T_CELLBORN) % items(i_item), i, j) + call add_map_element(tally_maps(FILTER_CELLBORN) % items(i_item), 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 + n = t % n_filter_bins(FILTER_CELL) + t % stride(FILTER_CELL) = filter_bins if (n > 0) then do j = 1, n i_item = t % cell_bins(j) % scalar - call add_map_element(tally_maps(T_CELL) % items(i_item), i, j) + call add_map_element(tally_maps(FILTER_CELL) % items(i_item), 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 + n = t % n_filter_bins(FILTER_MATERIAL) + t % stride(FILTER_MATERIAL) = filter_bins if (n > 0) then do j = 1, n i_item = t % material_bins(j) % scalar - call add_map_element(tally_maps(T_MATERIAL) % items(i_item), i, j) + call add_map_element(tally_maps(FILTER_MATERIAL) % items(i_item), 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 + n = t % n_filter_bins(FILTER_UNIVERSE) + t % stride(FILTER_UNIVERSE) = filter_bins if (n > 0) then do j = 1, n i_item = t % universe_bins(j) % scalar - call add_map_element(tally_maps(T_UNIVERSE) % items(i_item), i, j) + call add_map_element(tally_maps(FILTER_UNIVERSE) % items(i_item), i, j) end do filter_bins = filter_bins * n end if - ! Finally add scoring bins for the macro tallies - n = t % n_macro_bins + ! Finally add scoring bins for the tally + n = t % n_score_bins if (n > 0) then score_bins = n else - message = "Must have macro tally bins!" + message = "Must have scoring bins!" call fatal_error() end if @@ -232,34 +232,26 @@ contains end subroutine add_map_element !=============================================================================== -! SCORE_TALLY contains the main logic for scoring user-specified tallies +! SCORE_ANALOG_TALLY keeps track of how many events occur in a specified cell, +! energy range, etc. Note that since these are "analog" tallies, they are only +! triggered at every collision, not every event !=============================================================================== - subroutine score_tally(p, scattered, fissioned) + subroutine score_analog_tally(p) type(Particle), pointer :: p - logical, intent(in) :: scattered - logical, intent(in) :: fissioned integer :: i integer :: j - integer :: k - integer :: n - integer :: bins(TALLY_TYPES) - integer :: bin_energyout + integer :: bins(N_FILTER_TYPES) integer :: score_index - integer :: score_index0 - integer :: macro_bin - integer :: mesh_bin + integer :: score_bin real(8) :: score real(8) :: last_wgt real(8) :: wgt real(8) :: mu - real(8) :: E_out - logical :: has_energyout_bin - logical :: analog - type(TallyObject), pointer :: t - type(StructuredMesh), pointer :: m + logical :: found_bin + type(TallyObject), pointer :: t => null() ! Copy particle's pre- and post-collision weight and angle last_wgt = p % last_wgt @@ -273,89 +265,16 @@ contains t => tallies(i) ! Surface current tallies are treated separately - if (t % surface_current) cycle + if (t % type == TALLY_SURFACE_CURRENT) cycle + + ! Skip tallies that are handled by track-length estimators + if (t % estimator == ESTIMATOR_TRACKLENGTH) cycle ! ======================================================================= ! DETERMINE SCORING BIN COMBINATION - ! determine mesh bin - if (t % n_bins(T_MESH) > 0) then - m => meshes(t % mesh) - - ! Determine if we're in the mesh first - call get_mesh_bin(m, p % coord0 % xyz, mesh_bin) - if (mesh_bin == NO_BIN_FOUND) cycle - bins(T_MESH) = mesh_bin - else - bins(T_MESH) = 1 - end if - - ! determine next universe bin - if (t % n_bins(T_UNIVERSE) > 0) then - bins(T_UNIVERSE) = get_next_bin(T_UNIVERSE, p % coord % universe, i) - if (bins(T_UNIVERSE) == NO_BIN_FOUND) cycle - else - bins(T_UNIVERSE) = 1 - end if - - ! determine next material bin - 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 - bins(T_MATERIAL) = 1 - end if - - ! determine next cell bin - if (t % n_bins(T_CELL) > 0) then - bins(T_CELL) = get_next_bin(T_CELL, p % coord % cell, i) - if (bins(T_CELL) == NO_BIN_FOUND) cycle - else - bins(T_CELL) = 1 - end if - - ! 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 - bins(T_CELLBORN) = 1 - end if - - ! determine next surface bin - 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 - bins(T_SURFACE) = 1 - end if - - ! determine incoming energy bin - n = t % n_bins(T_ENERGYIN) - if (n > 0) then - ! check if energy of the particle is within energy bins - if (p % last_E < t % energy_in(1) .or. & - p % last_E > t % energy_in(n + 1)) cycle - - ! search to find incoming energy bin - bins(T_ENERGYIN) = binary_search(t % energy_in, n + 1, p % last_E) - else - bins(T_ENERGYIN) = 1 - end if - - ! determine outgoing energy bin - 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 + 1)) cycle - - ! search to find incoming energy bin - bins(T_ENERGYOUT) = binary_search(t % energy_out, n + 1, p % E) - has_energyout_bin = .true. - else - bins(T_ENERGYOUT) = 1 - has_energyout_bin = .false. - end if + call get_scoring_bins(p, i, bins, found_bin) + if (.not. found_bin) cycle ! ======================================================================= ! CALCULATE SCORES AND ACCUMULATE TALLY @@ -368,134 +287,161 @@ contains score_index = sum((bins - 1) * t % stride) + 1 ! Determine score for each bin - do j = 1, t % n_macro_bins - ! determine what type of macro bin - macro_bin = t % macro_bins(j) % scalar + do j = 1, t % n_score_bins + ! determine what type of score bin + score_bin = t % score_bins(j) % scalar - ! determine if we need outgoing angle - analog = (macro_bin == MACRO_NU_SCATTER .or. & - macro_bin == MACRO_SCATTER_1 .or. macro_bin == MACRO_SCATTER_2 .or. & - macro_bin == MACRO_SCATTER_3 .or. macro_bin == MACRO_N_1N .or. & - macro_bin == MACRO_N_2N .or. macro_bin == MACRO_N_3N .or. & - macro_bin == MACRO_N_4N) + select case (score_bin) + case (SCORE_FLUX) + ! All events score to a flux bin. We actually use a collision + ! estimator since there is no way to count 'events' exactly for the + ! flux - if (has_energyout_bin .or. analog) then - ! If this tally has an outgoing energy filter, the only supported - ! reaction is scattering or nu-fission. Note that some macro - ! quantities can only be scored in analog such as scattering - ! moments and (n,xn) + score = last_wgt / material_xs % total - if ((.not. scattered) .and. (.not. fissioned)) cycle + case (SCORE_TOTAL) + ! All events will score to the total reaction rate. We can just use + ! the weight of the particle entering the collision as the score - if (scattered) then - ! since scattering has already occured, we do not need to - ! multiply by the scattering cross section + score = last_wgt - select case (macro_bin) - case (MACRO_SCATTER) - score = last_wgt - case (MACRO_NU_SCATTER) - score = wgt - case (MACRO_SCATTER_1) - score = last_wgt * mu - case (MACRO_SCATTER_2) - score = last_wgt * 0.5*(3.0*mu*mu - ONE) - case (MACRO_SCATTER_3) - score = last_wgt * 0.5*(5.0*mu*mu*mu - 3.0*mu) - case (MACRO_N_1N) - if (wgt == last_wgt) then - score = last_wgt - else - cycle - end if - case (MACRO_N_2N) - if (int(wgt/last_wgt) == 2) then - score = last_wgt - else - cycle - end if - case (MACRO_N_3N) - if (int(wgt/last_wgt) == 3) then - score = last_wgt - else - cycle - end if - case (MACRO_N_4N) - if (int(wgt/last_wgt) == 4) then - score = last_wgt - else - cycle - end if - case (MACRO_NU_FISSION) - cycle - case default - message = "Invalid macro reaction on analog tally." - call fatal_error() - end select - end if + case (SCORE_SCATTER) + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle - if (fissioned) then + ! Since only scattering events make it here, again we can use the + ! weight entering the collision as the estimator for the reaction + ! rate - if (macro_bin /= MACRO_NU_FISSION) cycle + score = last_wgt + case (SCORE_NU_SCATTER) + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle + + ! For scattering production, we need to use the post-collision + ! weight as the estimate for the number of neutrons exiting a + ! reaction with neutrons in the exit channel + + score = wgt + + case (SCORE_SCATTER_1) + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle + + ! The first scattering moment can be determined by using the rate + ! of scattering reactions multiplied by the cosine of the change in + ! neutron's angle due to the collision + + score = last_wgt * mu + + case (SCORE_SCATTER_2) + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle + + ! The second scattering moment can be determined in a similar + ! manner to the first scattering moment + + score = last_wgt * 0.5*(3.0*mu*mu - ONE) + + case (SCORE_SCATTER_3) + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle + + ! The first scattering moment can be determined by using the rate + ! of scattering reactions multiplied by the cosine of the change in + ! neutron's angle due to the collision + + score = last_wgt * 0.5*(5.0*mu*mu*mu - 3.0*mu) + + case (SCORE_N_1N) + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle + + ! Skip any events where weight of particle changed + if (wgt /= last_wgt) cycle + + ! All events that reach this point are (n,1n) reactions + score = last_wgt + + case (SCORE_N_2N) + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle + + ! Skip any scattering events where the weight didn't double + if (nint(wgt/last_wgt) /= 2) cycle + + ! All events that reach this point are (n,2n) reactions + score = last_wgt + + case (SCORE_N_3N) + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle + + ! Skip any scattering events where the weight didn't double + if (nint(wgt/last_wgt) /= 3) cycle + + ! All events that reach this point are (n,3n) reactions + score = last_wgt + + case (SCORE_N_4N) + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle + + ! Skip any scattering events where the weight didn't double + if (nint(wgt/last_wgt) /= 4) cycle + + ! All events that reach this point are (n,3n) reactions + score = last_wgt + + case (SCORE_ABSORPTION) + ! Skip any event where the particle wasn't absorbed + if (p % event == EVENT_SCATTER) cycle + + ! All fission and absorption events will contribute here, so we can + ! just use the particle's weight entering the collision + + score = last_wgt + + case (SCORE_FISSION) + ! Skip any non-fission events + if (p % event /= EVENT_FISSION) cycle + + ! All fission events will contribute, so again we can use + ! particle's weight entering the collision as the estimate for the + ! fission reaction rate + + score = last_wgt + + case (SCORE_NU_FISSION) + ! Skip any non-fission events + if (p % event /= EVENT_FISSION) cycle + + if (t % n_filter_bins(FILTER_ENERGYOUT) > 0) then ! Normally, we only need to make contributions to one scoring ! bin. However, in the case of fission, since multiple fission - ! neutrons were emitted with different energies, multiple - ! outgoing energy bins may have been scored to. The following - ! logic treats this special case and scores to multiple bins + ! neutrons were emitted with different energies, multiple outgoing + ! energy bins may have been scored to. The following logic treats + ! this special case and scores to multiple bins - ! save original outgoing energy bin and score index - bin_energyout = bins(T_ENERGYOUT) - score_index0 = score_index - - ! Since the creation of fission sites is weighted such that it - ! is expected to create n_particles sites, we need to multiply - ! the score by keff to get the true nu-fission rate. Otherwise, - ! the sum of all nu-fission rates would be ~1.0. - - score = keff - - ! loop over number of particles banked - do k = 1, p % n_bank - ! determine outgoing energy from fission bank - E_out = fission_bank(n_bank - p % n_bank + k) % E - - ! change outgoing energy bin - bins(T_ENERGYOUT) = binary_search(t % energy_out, n + 1, E_out) - - ! determine scoring index - score_index = sum((bins - 1) * t % stride) + 1 - - ! Add score to tally - call add_to_score(t % scores(score_index, j), score) - end do - - ! reset outgoing energy bin and score index - bins(T_ENERGYOUT) = bin_energyout - score_index = score_index0 + call score_fission_eout(p, t, bins, j) cycle - end if - else - ! For tallies with no outgoing energy filter, the score is - ! calculated normally depending on the quantity specified + else + ! If there is no outgoing energy filter, than we only need to + ! score to one bin. For the score to be 'analog', we need to + ! score the number of particles that were banked in the fission + ! bank. Since this was weighted by 1/keff, we multiply by keff + ! to get the proper score. - select case (macro_bin) - case (MACRO_FLUX) - score = last_wgt / material_xs % total - case (MACRO_TOTAL) - score = last_wgt - case (MACRO_SCATTER) - score = last_wgt * (material_xs % total - material_xs % absorption) & - / material_xs % total - case (MACRO_ABSORPTION) - score = last_wgt * material_xs % absorption / material_xs % total - case (MACRO_FISSION) - score = last_wgt * material_xs % fission / material_xs % total - case (MACRO_NU_FISSION) - score = last_wgt * material_xs % nu_fission / material_xs % total - end select - end if + score = keff * p % n_bank + + end if + + case default + message = "Invalid score type on tally " // to_str(t % id) // "." + call fatal_error() + end select ! Add score to tally call add_to_score(t % scores(score_index, j), score) @@ -507,7 +453,261 @@ contains end do - end subroutine score_tally + end subroutine score_analog_tally + +!=============================================================================== +! SCORE_FISSION_EOUT handles a special case where we need to store neutron +! production rate with an outgoing energy filter (think of a fission matrix). In +! this case, we may need to score to multiple bins if there were multiple +! neutrons produced with different energies. +!=============================================================================== + + subroutine score_fission_eout(p, t, bins, j) + + type(Particle), pointer :: p + type(TallyObject), pointer :: t + integer, intent(inout) :: bins(N_FILTER_TYPES) + integer, intent(in) :: j + + integer :: k + integer :: bin_energyout + integer :: score_index + real(8) :: score + real(8) :: E_out + + ! save original outgoing energy bin and score index + bin_energyout = bins(FILTER_ENERGYOUT) + + ! Since the creation of fission sites is weighted such that it is + ! expected to create n_particles sites, we need to multiply the + ! score by keff to get the true nu-fission rate. Otherwise, the sum + ! of all nu-fission rates would be ~1.0. + + score = keff + + ! loop over number of particles banked + do k = 1, p % n_bank + ! determine outgoing energy from fission bank + E_out = fission_bank(n_bank - p % n_bank + k) % E + + ! change outgoing energy bin + bins(FILTER_ENERGYOUT) = binary_search(t % energy_out, & + size(t % energy_out), E_out) + + ! determine scoring index + score_index = sum((bins - 1) * t % stride) + 1 + + ! Add score to tally + call add_to_score(t % scores(score_index, j), score) + end do + + ! reset outgoing energy bin and score index + bins(FILTER_ENERGYOUT) = bin_energyout + + end subroutine score_fission_eout + +!=============================================================================== +! SCORE_TRACKLENGTH_TALLY calculates fluxes and reaction rates based on the +! track-length estimate of the flux. This is triggered at every event (surface +! crossing, lattice crossing, or collision) and thus cannot be done for tallies +! that require post-collision information. +!=============================================================================== + + subroutine score_tracklength_tally(p, distance) + + type(Particle), pointer :: p + real(8), intent(in) :: distance + + integer :: i + integer :: j + integer :: bins(N_FILTER_TYPES) + integer :: score_index + integer :: score_bin + real(8) :: flux + real(8) :: score + logical :: found_bin + type(TallyObject), pointer :: t => null() + + ! Determine track-length estimate of flux + flux = p % wgt * distance + + ! A loop over all tallies is necessary because we need to simultaneously + ! determine different filter bins for the same tally in order to score to it + + do i = 1, n_tallies + t => tallies(i) + + ! Surface current tallies are treated separately + if (t % type == TALLY_SURFACE_CURRENT) cycle + + ! Skip tallies that are handled by track-length estimators + if (t % estimator == ESTIMATOR_ANALOG) cycle + + ! ======================================================================= + ! DETERMINE SCORING BIN COMBINATION + + call get_scoring_bins(p, i, bins, found_bin) + if (.not. found_bin) cycle + + ! ======================================================================= + ! CALCULATE SCORES AND ACCUMULATE TALLY + + ! If we have made it here, we have a scoring combination of bins for this + ! tally -- now we need to determine where in the scores array we should + ! be accumulating the tally values + + ! Determine scoring index for this filter combination + score_index = sum((bins - 1) * t % stride) + 1 + + ! Determine score for each bin + do j = 1, t % n_score_bins + ! determine what type of score bin + score_bin = t % score_bins(j) % scalar + + ! Determine cross section + select case(score_bin) + case (SCORE_FLUX) + score = flux + case (SCORE_TOTAL) + score = material_xs % total * flux + case (SCORE_SCATTER) + score = (material_xs % total - material_xs % absorption) * flux + case (SCORE_ABSORPTION) + score = material_xs % absorption * flux + case (SCORE_FISSION) + score = material_xs % fission * flux + case (SCORE_NU_FISSION) + score = material_xs % nu_fission * flux + case default + message = "Invalid score type on tally " // to_str(t % id) // "." + call fatal_error() + end select + + ! Add score to tally + call add_to_score(t % scores(score_index, j), score) + + end do + + ! Reset tally map positioning + position = 0 + + end do + + end subroutine score_tracklength_tally + +!=============================================================================== +! GET_SCORING_BINS determines a combination of filter bins that should be scored +! for a tally based on the particle's current attributes. +!=============================================================================== + + subroutine get_scoring_bins(p, index_tally, bins, found_bin) + + type(Particle), pointer :: p + integer, intent(in) :: index_tally + integer, intent(out) :: bins(N_FILTER_TYPES) + logical, intent(out) :: found_bin + + integer :: i + integer :: n + integer :: mesh_bin + type(TallyObject), pointer :: t => null() + type(StructuredMesh), pointer :: m => null() + + found_bin = .true. + t => tallies(index_tally) + bins = 1 + + FILTER_LOOP: do i = 1, t % n_filters + + select case (t % filters(i)) + case (FILTER_MESH) + ! determine mesh bin + m => meshes(t % mesh) + + ! Determine if we're in the mesh first + call get_mesh_bin(m, p % coord0 % xyz, mesh_bin) + if (mesh_bin == NO_BIN_FOUND) then + found_bin = .false. + return + end if + bins(FILTER_MESH) = mesh_bin + + case (FILTER_UNIVERSE) + ! determine next universe bin + ! TODO: Account for multiple universes when performing this filter + bins(FILTER_UNIVERSE) = get_next_bin(FILTER_UNIVERSE, & + p % coord % universe, index_tally) + if (bins(FILTER_UNIVERSE) == NO_BIN_FOUND) then + found_bin = .false. + return + end if + + case (FILTER_MATERIAL) + bins(FILTER_MATERIAL) = get_next_bin(FILTER_MATERIAL, & + p % material, index_tally) + if (bins(FILTER_MATERIAL) == NO_BIN_FOUND) then + found_bin = .false. + return + end if + + case (FILTER_CELL) + ! determine next cell bin + ! TODO: Account for cells in multiple levels when performing this filter + bins(FILTER_CELL) = get_next_bin(FILTER_CELL, & + p % coord % cell, index_tally) + if (bins(FILTER_CELL) == NO_BIN_FOUND) then + found_bin = .false. + return + end if + + case (FILTER_CELLBORN) + ! determine next cellborn bin + bins(FILTER_CELLBORN) = get_next_bin(FILTER_CELLBORN, & + p % cell_born, index_tally) + if (bins(FILTER_CELLBORN) == NO_BIN_FOUND) then + found_bin = .false. + return + end if + + case (FILTER_SURFACE) + ! determine next surface bin + bins(FILTER_SURFACE) = get_next_bin(FILTER_SURFACE, & + p % surface, index_tally) + if (bins(FILTER_SURFACE) == NO_BIN_FOUND) then + found_bin = .false. + return + end if + + case (FILTER_ENERGYIN) + ! determine incoming energy bin + n = t % n_filter_bins(FILTER_ENERGYIN) + ! check if energy of the particle is within energy bins + if (p % last_E < t % energy_in(1) .or. & + p % last_E > t % energy_in(n + 1)) then + found_bin = .false. + return + end if + + ! search to find incoming energy bin + bins(FILTER_ENERGYIN) = binary_search(t % energy_in, n + 1, p % last_E) + + case (FILTER_ENERGYOUT) + ! determine outgoing energy bin + n = t % n_filter_bins(FILTER_ENERGYOUT) + ! check if energy of the particle is within energy bins + if (p % E < t % energy_out(1) .or. p % E > t % energy_out(n + 1)) then + found_bin = .false. + return + end if + + ! search to find incoming energy bin + bins(FILTER_ENERGYOUT) = binary_search(t % energy_out, n + 1, p % E) + + end select + + end do FILTER_LOOP + + end subroutine get_scoring_bins !=============================================================================== ! SCORE_SURFACE_CURRENT tallies surface crossings in a mesh tally by manually @@ -518,26 +718,26 @@ contains type(Particle), pointer :: p - integer :: i ! loop indices - integer :: j ! loop indices - integer :: k ! loop indices - integer :: ijk0(3) ! indices of starting coordinates - integer :: ijk1(3) ! indices of ending coordinates - integer :: n_cross ! number of surface crossings - integer :: n ! number of incoming energy bins - integer :: bins(TALLY_TYPES) ! scoring bin combination - integer :: score_index ! index of scoring bin - real(8) :: uvw(3) ! cosine of angle of particle - real(8) :: xyz0(3) ! starting/intermediate coordinates - real(8) :: xyz1(3) ! ending coordinates of particle - real(8) :: xyz_cross(3) ! coordinates of bounding surfaces - real(8) :: d(3) ! distance to each bounding surface - real(8) :: distance ! actual distance traveled - logical :: start_in_mesh ! particle's starting xyz in mesh? - logical :: end_in_mesh ! particle's ending xyz in mesh? - logical :: x_same ! same starting/ending x index (i) - logical :: y_same ! same starting/ending y index (j) - logical :: z_same ! same starting/ending z index (k) + integer :: i ! loop indices + integer :: j ! loop indices + integer :: k ! loop indices + integer :: ijk0(3) ! indices of starting coordinates + integer :: ijk1(3) ! indices of ending coordinates + integer :: n_cross ! number of surface crossings + integer :: n ! number of incoming energy bins + integer :: bins(N_FILTER_TYPES) ! scoring bin combination + integer :: score_index ! index of scoring bin + real(8) :: uvw(3) ! cosine of angle of particle + real(8) :: xyz0(3) ! starting/intermediate coordinates + real(8) :: xyz1(3) ! ending coordinates of particle + real(8) :: xyz_cross(3) ! coordinates of bounding surfaces + real(8) :: d(3) ! distance to each bounding surface + real(8) :: distance ! actual distance traveled + logical :: start_in_mesh ! particle's starting xyz in mesh? + logical :: end_in_mesh ! particle's ending xyz in mesh? + logical :: x_same ! same starting/ending x index (i) + logical :: y_same ! same starting/ending y index (j) + logical :: z_same ! same starting/ending z index (k) type(TallyObject), pointer :: t => null() type(StructuredMesh), pointer :: m => null() @@ -552,7 +752,7 @@ contains t => tallies(i) ! Skip non-surface-current tallies - if (.not. t % surface_current) cycle + if (t % type /= TALLY_SURFACE_CURRENT) cycle ! Determine indices for starting and ending location m => meshes(t % mesh) @@ -570,16 +770,16 @@ contains uvw = p % coord0 % uvw ! determine incoming energy bin - n = t % n_bins(T_ENERGYIN) + n = t % n_filter_bins(FILTER_ENERGYIN) if (n > 0) then ! check if energy of the particle is within energy bins if (p % last_E < t % energy_in(1) .or. & p % last_E > t % energy_in(n + 1)) cycle ! search to find incoming energy bin - bins(TS_ENERGYIN) = binary_search(t % energy_in, n + 1, p % last_E) + bins(SURF_FILTER_ENERGYIN) = binary_search(t % energy_in, n + 1, p % last_E) else - bins(TS_ENERGYIN) = 1 + bins(SURF_FILTER_ENERGYIN) = 1 end if ! ======================================================================= @@ -595,7 +795,7 @@ contains do j = ijk0(3), ijk1(3) - 1 ijk0(3) = j if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - bins(TS_SURFACE) = OUT_TOP + bins(SURF_FILTER_SURFACE) = OUT_TOP bins(1:3) = ijk0 + 1 score_index = sum((bins - 1) * t % stride) + 1 call add_to_score(t % scores(score_index, 1), p % wgt) @@ -605,7 +805,7 @@ contains do j = ijk0(3) - 1, ijk1(3), -1 ijk0(3) = j if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - bins(TS_SURFACE) = IN_TOP + bins(SURF_FILTER_SURFACE) = IN_TOP bins(1:3) = ijk0 + 1 score_index = sum((bins - 1) * t % stride) + 1 call add_to_score(t % scores(score_index, 1), p % wgt) @@ -619,7 +819,7 @@ contains do j = ijk0(2), ijk1(2) - 1 ijk0(2) = j if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - bins(TS_SURFACE) = OUT_FRONT + bins(SURF_FILTER_SURFACE) = OUT_FRONT bins(1:3) = ijk0 + 1 score_index = sum((bins - 1) * t % stride) + 1 call add_to_score(t % scores(score_index, 1), p % wgt) @@ -629,7 +829,7 @@ contains do j = ijk0(2) - 1, ijk1(2), -1 ijk0(2) = j if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - bins(TS_SURFACE) = IN_FRONT + bins(SURF_FILTER_SURFACE) = IN_FRONT bins(1:3) = ijk0 + 1 score_index = sum((bins - 1) * t % stride) + 1 call add_to_score(t % scores(score_index, 1), p % wgt) @@ -643,7 +843,7 @@ contains do j = ijk0(1), ijk1(1) - 1 ijk0(1) = j if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - bins(TS_SURFACE) = OUT_RIGHT + bins(SURF_FILTER_SURFACE) = OUT_RIGHT bins(1:3) = ijk0 + 1 score_index = sum((bins - 1) * t % stride) + 1 call add_to_score(t % scores(score_index, 1), p % wgt) @@ -653,7 +853,7 @@ contains do j = ijk0(1) - 1, ijk1(1), -1 ijk0(1) = j if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - bins(TS_SURFACE) = IN_RIGHT + bins(SURF_FILTER_SURFACE) = IN_RIGHT bins(1:3) = ijk0 + 1 score_index = sum((bins - 1) * t % stride) + 1 call add_to_score(t % scores(score_index, 1), p % wgt) @@ -677,7 +877,7 @@ contains do k = 1, n_cross ! Reset scoring bin index - bins(TS_SURFACE) = 0 + bins(SURF_FILTER_SURFACE) = 0 ! Calculate distance to each bounding surface. We need to treat ! special case where the cosine of the angle is zero since this would @@ -704,7 +904,7 @@ contains ! Crossing into right mesh cell -- this is treated as outgoing ! current from (i,j,k) if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - bins(TS_SURFACE) = OUT_RIGHT + bins(SURF_FILTER_SURFACE) = OUT_RIGHT bins(1:3) = ijk0 + 1 end if ijk0(1) = ijk0(1) + 1 @@ -715,7 +915,7 @@ contains ijk0(1) = ijk0(1) - 1 xyz_cross(1) = xyz_cross(1) - m % width(1) if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - bins(TS_SURFACE) = IN_RIGHT + bins(SURF_FILTER_SURFACE) = IN_RIGHT bins(1:3) = ijk0 + 1 end if end if @@ -724,7 +924,7 @@ contains ! Crossing into front mesh cell -- this is treated as outgoing ! current in (i,j,k) if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - bins(TS_SURFACE) = OUT_FRONT + bins(SURF_FILTER_SURFACE) = OUT_FRONT bins(1:3) = ijk0 + 1 end if ijk0(2) = ijk0(2) + 1 @@ -735,7 +935,7 @@ contains ijk0(2) = ijk0(2) - 1 xyz_cross(2) = xyz_cross(2) - m % width(2) if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - bins(TS_SURFACE) = IN_FRONT + bins(SURF_FILTER_SURFACE) = IN_FRONT bins(1:3) = ijk0 + 1 end if end if @@ -744,7 +944,7 @@ contains ! Crossing into top mesh cell -- this is treated as outgoing ! current in (i,j,k) if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - bins(TS_SURFACE) = OUT_TOP + bins(SURF_FILTER_SURFACE) = OUT_TOP bins(1:3) = ijk0 + 1 end if ijk0(3) = ijk0(3) + 1 @@ -755,14 +955,14 @@ contains ijk0(3) = ijk0(3) - 1 xyz_cross(3) = xyz_cross(3) - m % width(3) if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - bins(TS_SURFACE) = IN_TOP + bins(SURF_FILTER_SURFACE) = IN_TOP bins(1:3) = ijk0 + 1 end if end if end if ! Determine scoring index - if (bins(TS_SURFACE) > 0) then + if (bins(SURF_FILTER_SURFACE) > 0) then score_index = sum((bins - 1) * t % stride) + 1 ! Check for errors @@ -879,7 +1079,7 @@ contains ! Loop over all filter and scoring bins do j = 1, t % n_total_bins - do k = 1, t % n_macro_bins + do k = 1, t % n_score_bins ! Add the sum and square of the sum of contributions from each ! history within a cycle to the variables val and val_sq. This will ! later allow us to calculate a variance on the tallies @@ -905,49 +1105,49 @@ contains subroutine write_tallies() - integer :: i ! index in tallies array - integer :: j ! level in tally hierarchy - integer :: k ! loop index for scoring bins - integer :: bins(TALLY_TYPES) = 0 ! bins corresponding to each filter - integer :: indent ! number of spaces to preceed output - integer :: io_error ! error in opening/writing file - integer :: last_filter ! lowest level filter type - integer :: score_index ! index in scores array for filters - logical :: file_exists ! does tallies.out file already exists? - logical :: has_filter(TALLY_TYPES) ! does tally have this filter? - character(MAX_FILE_LEN) :: filename ! name of output file - character(15) :: filter_name(TALLY_TYPES) ! names of tally filters - character(27) :: macro_name(N_MACRO_TYPES) ! names of macro scores + integer :: i ! index in tallies array + integer :: j ! level in tally hierarchy + integer :: k ! loop index for scoring bins + integer :: bins(N_FILTER_TYPES) = 0 ! bins corresponding to each filter + integer :: indent ! number of spaces to preceed output + integer :: io_error ! error in opening/writing file + integer :: last_filter ! lowest level filter type + integer :: score_index ! index in scores array for filters + logical :: file_exists ! does tallies.out file already exists? + logical :: has_filter(N_FILTER_TYPES) ! does tally have this filter? + character(MAX_FILE_LEN) :: filename ! name of output file + character(15) :: filter_name(N_FILTER_TYPES) ! names of tally filters + character(27) :: score_name(N_SCORE_TYPES) ! names of scoring function type(TallyObject), pointer :: t ! Skip if there are no tallies if (n_tallies == 0) return ! Initialize names for tally filter types - filter_name(T_UNIVERSE) = "Universe" - filter_name(T_MATERIAL) = "Material" - filter_name(T_CELL) = "Cell" - filter_name(T_CELLBORN) = "Birth Cell" - filter_name(T_SURFACE) = "Surface" - filter_name(T_MESH) = "Mesh" - filter_name(T_ENERGYIN) = "Incoming Energy" - filter_name(T_ENERGYOUT) = "Outgoing Energy" + filter_name(FILTER_UNIVERSE) = "Universe" + filter_name(FILTER_MATERIAL) = "Material" + filter_name(FILTER_CELL) = "Cell" + filter_name(FILTER_CELLBORN) = "Birth Cell" + filter_name(FILTER_SURFACE) = "Surface" + filter_name(FILTER_MESH) = "Mesh" + filter_name(FILTER_ENERGYIN) = "Incoming Energy" + filter_name(FILTER_ENERGYOUT) = "Outgoing Energy" - ! Initialize names for macro scores - macro_name(abs(MACRO_FLUX)) = "Flux" - macro_name(abs(MACRO_TOTAL)) = "Total Reaction Rate" - macro_name(abs(MACRO_SCATTER)) = "Scattering Rate" - macro_name(abs(MACRO_NU_SCATTER)) = "Scattering Production Rate" - macro_name(abs(MACRO_SCATTER_1)) = "First Scattering Moment" - macro_name(abs(MACRO_SCATTER_2)) = "Second Scattering Moment" - macro_name(abs(MACRO_SCATTER_3)) = "Third Scattering Moment" - macro_name(abs(MACRO_N_1N)) = "(n,1n) Rate" - macro_name(abs(MACRO_N_2N)) = "(n,2n) Rate" - macro_name(abs(MACRO_N_3N)) = "(n,3n) Rate" - macro_name(abs(MACRO_N_4N)) = "(n,4n) Rate" - macro_name(abs(MACRO_ABSORPTION)) = "Absorption Rate" - macro_name(abs(MACRO_FISSION)) = "Fission Rate" - macro_name(abs(MACRO_NU_FISSION)) = "Nu-Fission Rate" + ! Initialize names for scores + score_name(abs(SCORE_FLUX)) = "Flux" + score_name(abs(SCORE_TOTAL)) = "Total Reaction Rate" + score_name(abs(SCORE_SCATTER)) = "Scattering Rate" + score_name(abs(SCORE_NU_SCATTER)) = "Scattering Production Rate" + score_name(abs(SCORE_SCATTER_1)) = "First Scattering Moment" + score_name(abs(SCORE_SCATTER_2)) = "Second Scattering Moment" + score_name(abs(SCORE_SCATTER_3)) = "Third Scattering Moment" + score_name(abs(SCORE_N_1N)) = "(n,1n) Rate" + score_name(abs(SCORE_N_2N)) = "(n,2n) Rate" + score_name(abs(SCORE_N_3N)) = "(n,3n) Rate" + score_name(abs(SCORE_N_4N)) = "(n,4n) Rate" + score_name(abs(SCORE_ABSORPTION)) = "Absorption Rate" + score_name(abs(SCORE_FISSION)) = "Fission Rate" + score_name(abs(SCORE_NU_FISSION)) = "Nu-Fission Rate" ! Create filename for tally output filename = trim(path_input) // "tallies.out" @@ -969,14 +1169,14 @@ contains call header("TALLY " // trim(to_str(t % id)), unit=UNIT_TALLY, level=3) ! Handle surface current tallies separately - if (t % surface_current) then + if (t % type == TALLY_SURFACE_CURRENT) then call write_surface_current(t) cycle end if ! First determine which filters this tally has - do j = 1, TALLY_TYPES - if (t % n_bins(j) > 0) then + do j = 1, N_FILTER_TYPES + if (t % n_filter_bins(j) > 0) then has_filter(j) = .true. last_filter = j else @@ -1003,7 +1203,7 @@ contains ! ================================================================= ! REACHED END OF BINS FOR THIS FILTER, MOVE TO NEXT FILTER - if ((has_filter(j) .and. bins(j) > t % n_bins(j)) .or. & + if ((has_filter(j) .and. bins(j) > t % n_filter_bins(j)) .or. & ((.not. has_filter(j)) .and. bins(j) > 1)) then if (j == 1) then ! This means we are done with all bin combinations @@ -1045,9 +1245,9 @@ contains ! Write scores for this filter bin combination indent = indent + 2 - do k = 1, t % n_macro_bins + do k = 1, t % n_score_bins write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & - repeat(" ", indent), macro_name(abs(t % macro_bins(k) % scalar)), & + repeat(" ", indent), score_name(abs(t % score_bins(k) % scalar)), & to_str(t % scores(score_index,k) % val), & trim(to_str(t % scores(score_index,k) % val_sq)) end do @@ -1070,16 +1270,16 @@ contains type(TallyObject), pointer :: t - integer :: i ! mesh index for x - integer :: j ! mesh index for y - integer :: k ! mesh index for z - integer :: l ! mesh index for energy - integer :: bins(TALLY_TYPES) ! bin combination - integer :: n ! number of incoming energy bins - integer :: len1 ! length of string - integer :: len2 ! length of string - integer :: score_index ! index in scores array for filters - logical :: print_ebin ! should incoming energy bin be displayed? + integer :: i ! mesh index for x + integer :: j ! mesh index for y + integer :: k ! mesh index for z + integer :: l ! mesh index for energy + integer :: bins(N_FILTER_TYPES) ! bin combination + integer :: n ! number of incoming energy bins + integer :: len1 ! length of string + integer :: len2 ! length of string + integer :: score_index ! index in scores array for filters + logical :: print_ebin ! should incoming energy bin be displayed? character(MAX_LINE_LEN) :: string type(StructuredMesh), pointer :: m => null() @@ -1090,7 +1290,7 @@ contains bins = 1 ! determine how many energy in bins there are - n = t % n_bins(T_ENERGYIN) + n = t % n_filter_bins(FILTER_ENERGYIN) if (n > 0) then print_ebin = .true. else @@ -1113,22 +1313,22 @@ contains ! Write incoming energy bin if (print_ebin) then write(UNIT=UNIT_TALLY, FMT='(3X,A,1X,A)') & - "Incoming Energy", trim(get_label(t, T_ENERGYIN, l)) + "Incoming Energy", trim(get_label(t, FILTER_ENERGYIN, l)) end if ! Set incoming energy bin - bins(TS_ENERGYIN) = l + bins(SURF_FILTER_ENERGYIN) = l ! Left Surface bins(1:3) = (/ i-1, j, k /) + 1 - bins(TS_SURFACE) = IN_RIGHT + bins(SURF_FILTER_SURFACE) = IN_RIGHT score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Left", & to_str(t % scores(score_index,1) % val), & trim(to_str(t % scores(score_index,1) % val_sq)) - bins(TS_SURFACE) = OUT_RIGHT + bins(SURF_FILTER_SURFACE) = OUT_RIGHT score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Left", & @@ -1137,14 +1337,14 @@ contains ! Right Surface bins(1:3) = (/ i, j, k /) + 1 - bins(TS_SURFACE) = IN_RIGHT + bins(SURF_FILTER_SURFACE) = IN_RIGHT score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Right", & to_str(t % scores(score_index,1) % val), & trim(to_str(t % scores(score_index,1) % val_sq)) - bins(TS_SURFACE) = OUT_RIGHT + bins(SURF_FILTER_SURFACE) = OUT_RIGHT score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Right", & @@ -1153,14 +1353,14 @@ contains ! Back Surface bins(1:3) = (/ i, j-1, k /) + 1 - bins(TS_SURFACE) = IN_FRONT + bins(SURF_FILTER_SURFACE) = IN_FRONT score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Back", & to_str(t % scores(score_index,1) % val), & trim(to_str(t % scores(score_index,1) % val_sq)) - bins(TS_SURFACE) = OUT_FRONT + bins(SURF_FILTER_SURFACE) = OUT_FRONT score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Back", & @@ -1169,14 +1369,14 @@ contains ! Front Surface bins(1:3) = (/ i, j, k /) + 1 - bins(TS_SURFACE) = IN_FRONT + bins(SURF_FILTER_SURFACE) = IN_FRONT score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Front", & to_str(t % scores(score_index,1) % val), & trim(to_str(t % scores(score_index,1) % val_sq)) - bins(TS_SURFACE) = OUT_FRONT + bins(SURF_FILTER_SURFACE) = OUT_FRONT score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Front", & @@ -1185,14 +1385,14 @@ contains ! Bottom Surface bins(1:3) = (/ i, j, k-1 /) + 1 - bins(TS_SURFACE) = IN_TOP + bins(SURF_FILTER_SURFACE) = IN_TOP score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Bottom", & to_str(t % scores(score_index,1) % val), & trim(to_str(t % scores(score_index,1) % val_sq)) - bins(TS_SURFACE) = OUT_TOP + bins(SURF_FILTER_SURFACE) = OUT_TOP score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Bottom", & @@ -1201,14 +1401,14 @@ contains ! Top Surface bins(1:3) = (/ i, j, k /) + 1 - bins(TS_SURFACE) = IN_TOP + bins(SURF_FILTER_SURFACE) = IN_TOP score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Top", & to_str(t % scores(score_index,1) % val), & trim(to_str(t % scores(score_index,1) % val_sq)) - bins(TS_SURFACE) = OUT_TOP + bins(SURF_FILTER_SURFACE) = OUT_TOP score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Top", & @@ -1241,22 +1441,22 @@ contains type(StructuredMesh), pointer :: m select case(filter_type) - case (T_UNIVERSE) + case (FILTER_UNIVERSE) i = t % universe_bins(bin) % scalar label = to_str(universes(i) % id) - case (T_MATERIAL) + case (FILTER_MATERIAL) i = t % material_bins(bin) % scalar label = to_str(materials(i) % id) - case (T_CELL) + case (FILTER_CELL) i = t % cell_bins(bin) % scalar label = to_str(cells(i) % id) - case (T_CELLBORN) + case (FILTER_CELLBORN) i = t % cellborn_bins(bin) % scalar label = to_str(cells(i) % id) - case (T_SURFACE) + case (FILTER_SURFACE) i = t % surface_bins(bin) % scalar label = to_str(surfaces(i) % id) - case (T_MESH) + case (FILTER_MESH) m => meshes(t % mesh) allocate(ijk(m % n_dimension)) call bin_to_mesh_indices(m, bin, ijk) @@ -1267,11 +1467,11 @@ contains label = "Index (" // trim(to_str(ijk(1))) // ", " // & trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")" end if - case (T_ENERGYIN) + case (FILTER_ENERGYIN) E0 = t % energy_in(bin) E1 = t % energy_in(bin + 1) label = "[" // trim(to_str(E0)) // ", " // trim(to_str(E1)) // ")" - case (T_ENERGYOUT) + case (FILTER_ENERGYOUT) E0 = t % energy_out(bin) E1 = t % energy_out(bin + 1) label = "[" // trim(to_str(E0)) // ", " // trim(to_str(E1)) // ")" @@ -1304,7 +1504,7 @@ contains t => tallies(i) do j = 1, t % n_total_bins - do k = 1, t % n_macro_bins + do k = 1, t % n_score_bins ! Copy values from tallies val = t % scores(j,k) % val val2 = t % scores(j,k) % val_sq diff --git a/src/tally_header.F90 b/src/tally_header.F90 index 832b12b70f..b722c2730c 100644 --- a/src/tally_header.F90 +++ b/src/tally_header.F90 @@ -58,12 +58,17 @@ module tally_header type TallyObject ! Basic data - integer :: id - integer :: type - real(8) :: volume - logical :: surface_current = .false. + integer :: id ! user-defined identifier + integer :: type ! volume, surface current + integer :: estimator ! collision, track-length + real(8) :: volume ! volume of region - ! Tally bin specifications + ! Information about what filters should be used + + integer :: n_filters + integer, allocatable :: filters(:) + + ! Filter bin specifications type(TallyFilter), pointer :: universe_bins(:) => null() type(TallyFilter), pointer :: material_bins(:) => null() @@ -74,24 +79,24 @@ module tally_header real(8), allocatable :: energy_in(:) real(8), allocatable :: energy_out(:) - ! Number of bins for each filter - integer :: n_total_bins = 0 + ! Total number of filter bins + 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 + ! for each filter type, e.g. n_filter_bins(FILTER_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 :: n_filter_bins(:) integer, allocatable :: stride(:) ! Macroscopic properties to score - type(TallyFilter), pointer :: macro_bins(:) => null() - integer :: n_macro_bins = 0 + type(TallyFilter), pointer :: score_bins(:) => null() + integer :: n_score_bins = 0 ! Scores for each bin -- the most natural way to have scores would be to ! have a dimension for each different type of bin, but older Fortran diff --git a/src/utils/analyze_source.py b/src/utils/analyze_source.py index 7acd5df45d..b382103283 100755 --- a/src/utils/analyze_source.py +++ b/src/utils/analyze_source.py @@ -34,6 +34,8 @@ for sourceFile in source: commentChar = '!' elif ending == 'f90': commentChar = '!' + elif ending == 'F90': + commentChar = '!' for line in open(sourceFile, 'r'): line = line.strip()