From fa5934cb037520c27bb168fa194701703e367651 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Mon, 3 Dec 2012 20:48:14 -0500 Subject: [PATCH] Added the capability to handle an arbitrary number of legendre moments. This required changes to the input, the output, and the tally handling. The input and math routines are not yet complete, as it still will not allow for any value for the scattering order. I do not expect to extend this past n=5 right now, but I want to have the framework in place so that an extension only requires changing an if-statement in input.xml and adding the additional legendre polynomials to math.F90 --- src/cmfd_input.F90 | 20 +- src/constants.F90 | 31 ++- src/input_xml.F90 | 490 +++++++++++++++++++++++-------------------- src/output.F90 | 352 +++++++++++++++++-------------- src/tally.F90 | 188 ++++++++--------- src/tally_header.F90 | 5 + 6 files changed, 578 insertions(+), 508 deletions(-) diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 51e8a1aada..cd1d95fd17 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -369,12 +369,16 @@ contains ! allocate scoring bins allocate(t % score_bins(4)) t % n_score_bins = 4 - + ! allocate scattering order data + allocate(t % scatt_order(4)) + t % scatt_order = 0 + ! set macro_bins - t % score_bins(1) = SCORE_FLUX - t % score_bins(2) = SCORE_TOTAL - t % score_bins(3) = SCORE_SCATTER_1 - t % score_bins(4) = SCORE_DIFFUSION + t % score_bins(1) = SCORE_FLUX + t % score_bins(2) = SCORE_TOTAL + t % score_bins(3) = SCORE_SCATTER_N + t % scatt_order(3) = 1 + t % score_bins(4) = SCORE_DIFFUSION ! Increment the appropriate index and set pointer analog_tallies(n_user_analog_tallies + 1) = i @@ -413,6 +417,9 @@ contains ! allocate macro reactions allocate(t % score_bins(2)) t % n_score_bins = 2 + ! allocate scattering order data + allocate(t % scatt_order(2)) + t % scatt_order = 0 ! set macro_bins t % score_bins(1) = SCORE_NU_SCATTER @@ -454,6 +461,9 @@ contains ! allocate macro reactions allocate(t % score_bins(1)) t % n_score_bins = 1 + ! allocate scattering order data + allocate(t % scatt_order(1)) + t % scatt_order = 0 ! set macro bins t % score_bins(1) = SCORE_CURRENT diff --git a/src/constants.F90 b/src/constants.F90 index e1d8e853a8..eebb90df67 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -251,28 +251,25 @@ module constants EVENT_FISSION = 3 ! Tally score type - integer, parameter :: N_SCORE_TYPES = 20 + integer, parameter :: N_SCORE_TYPES = 17 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_SCATTER_4 = -8, & ! fourth scattering moment - SCORE_SCATTER_5 = -9, & ! fifth scattering moment - SCORE_TRANSPORT = -10, & ! transport reaction rate - SCORE_DIFFUSION = -11, & ! diffusion coefficient - SCORE_N_1N = -12, & ! (n,1n) rate - SCORE_N_2N = -13, & ! (n,2n) rate - SCORE_N_3N = -14, & ! (n,3n) rate - SCORE_N_4N = -15, & ! (n,4n) rate - SCORE_ABSORPTION = -16, & ! absorption rate - SCORE_FISSION = -17, & ! fission rate - SCORE_NU_FISSION = -18, & ! neutron production rate - SCORE_CURRENT = -19, & ! partial current - SCORE_EVENTS = -20 ! number of events + SCORE_SCATTER_N = -5, & ! arbitrary scattering moment + SCORE_SCATTER_PN = -6, & ! system for scoring 0th through nth moment + SCORE_TRANSPORT = -7, & ! transport reaction rate + SCORE_DIFFUSION = -8, & ! diffusion coefficient + SCORE_N_1N = -9, & ! (n,1n) rate + SCORE_N_2N = -10, & ! (n,2n) rate + SCORE_N_3N = -11, & ! (n,3n) rate + SCORE_N_4N = -12, & ! (n,4n) rate + SCORE_ABSORPTION = -13, & ! absorption rate + SCORE_FISSION = -14, & ! fission rate + SCORE_NU_FISSION = -15, & ! neutron production rate + SCORE_CURRENT = -16, & ! partial current + SCORE_EVENTS = -17 ! number of events ! Tally map bin finding integer, parameter :: NO_BIN_FOUND = -1 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index fd5b90ebda..38873162bc 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -19,7 +19,7 @@ module input_xml implicit none type(DictionaryII), pointer :: & ! used to count how many cells each - cells_in_univ_dict => null() ! universe contains + cells_in_univ_dict => null() ! universe contains contains @@ -92,7 +92,7 @@ contains call get_environment_variable("CROSS_SECTIONS", env_variable) if (len_trim(env_variable) == 0) then message = "No cross_sections.xml file was specified in settings.xml & - &or in the CROSS_SECTIONS environment variable." + &or in the CROSS_SECTIONS environment variable." call fatal_error() else path_cross_sections = trim(env_variable) @@ -104,7 +104,7 @@ contains ! Make sure that either criticality or fixed source was specified if (criticality_ % batches == 0 .and. fixed_source_ % batches == 0) then message = "Number of batches on or " & - // "tag was zero." + // "tag was zero." call fatal_error() end if @@ -202,7 +202,7 @@ contains inquire(FILE=path_source, EXIST=file_exists) if (.not. file_exists) then message = "Binary source file '" // trim(path_source) // & - "' does not exist!" + "' does not exist!" call fatal_error() end if @@ -221,7 +221,7 @@ contains coeffs_reqd = 3 case default message = "Invalid spatial distribution for external source: " & - // trim(source_ % space % type) + // trim(source_ % space % type) call fatal_error() end select @@ -235,11 +235,11 @@ contains ! Read parameters for spatial distribution if (n < coeffs_reqd) then message = "Not enough parameters specified for spatial " & - // "distribution of external source." + // "distribution of external source." call fatal_error() elseif (n > coeffs_reqd) then message = "Too many parameters specified for spatial " & - // "distribution of external source." + // "distribution of external source." call fatal_error() elseif (n > 0) then allocate(external_source % params_space(n)) @@ -266,7 +266,7 @@ contains external_source % type_angle = SRC_ANGLE_TABULAR case default message = "Invalid angular distribution for external source: " & - // trim(source_ % angle % type) + // trim(source_ % angle % type) call fatal_error() end select @@ -280,11 +280,11 @@ contains ! Read parameters for angle distribution if (n < coeffs_reqd) then message = "Not enough parameters specified for angle " & - // "distribution of external source." + // "distribution of external source." call fatal_error() elseif (n > coeffs_reqd) then message = "Too many parameters specified for angle " & - // "distribution of external source." + // "distribution of external source." call fatal_error() elseif (n > 0) then allocate(external_source % params_angle(n)) @@ -314,7 +314,7 @@ contains external_source % type_energy = SRC_ENERGY_TABULAR case default message = "Invalid energy distribution for external source: " & - // trim(source_ % energy % type) + // trim(source_ % energy % type) call fatal_error() end select @@ -328,11 +328,11 @@ contains ! Read parameters for energy distribution if (n < coeffs_reqd) then message = "Not enough parameters specified for energy " & - // "distribution of external source." + // "distribution of external source." call fatal_error() elseif (n > coeffs_reqd) then message = "Too many parameters specified for energy " & - // "distribution of external source." + // "distribution of external source." call fatal_error() elseif (n > 0) then allocate(external_source % params_energy(n)) @@ -370,10 +370,10 @@ contains ! Check to make sure enough values were supplied if (size(entropy_(1) % lower_left) /= 3) then message = "Need to specify (x,y,z) coordinates of lower-left corner & - &of Shannon entropy mesh." + &of Shannon entropy mesh." elseif (size(entropy_(1) % upper_right) /= 3) then message = "Need to specify (x,y,z) coordinates of upper-right corner & - &of Shannon entropy mesh." + &of Shannon entropy mesh." end if ! Allocate mesh object and coordinates on mesh @@ -388,7 +388,7 @@ contains ! Check on values provided if (.not. all(entropy_mesh % upper_right > entropy_mesh % lower_left)) then message = "Upper-right coordinate must be greater than lower-left & - &coordinate for Shannon entropy mesh." + &coordinate for Shannon entropy mesh." call fatal_error() end if @@ -399,7 +399,7 @@ contains ! If so, make sure proper number of values were given if (size(entropy_(1) % dimension) /= 3) then message = "Dimension of entropy mesh must be given as three & - &integers." + &integers." call fatal_error() end if @@ -420,13 +420,13 @@ contains ! Check to make sure enough values were supplied if (size(uniform_fs_(1) % lower_left) /= 3) then message = "Need to specify (x,y,z) coordinates of lower-left corner & - &of UFS mesh." + &of UFS mesh." elseif (size(uniform_fs_(1) % upper_right) /= 3) then message = "Need to specify (x,y,z) coordinates of upper-right corner & - &of UFS mesh." + &of UFS mesh." elseif (size(uniform_fs_(1) % dimension) /= 3) then message = "Dimension of UFS mesh must be given as three & - &integers." + &integers." call fatal_error() end if @@ -450,13 +450,13 @@ contains ! Check on values provided if (.not. all(ufs_mesh % upper_right > ufs_mesh % lower_left)) then message = "Upper-right coordinate must be greater than lower-left & - &coordinate for UFS mesh." + &coordinate for UFS mesh." call fatal_error() end if ! Calculate width ufs_mesh % width = (ufs_mesh % upper_right - & - ufs_mesh % lower_left) / ufs_mesh % dimension + ufs_mesh % lower_left) / ufs_mesh % dimension ! Calculate volume fraction of each cell ufs_mesh % volume_frac = ONE/real(product(ufs_mesh % dimension),8) @@ -466,7 +466,7 @@ contains ! Allocate source_frac allocate(source_frac(1, ufs_mesh % dimension(1), & - ufs_mesh % dimension(2), ufs_mesh % dimension(3))) + ufs_mesh % dimension(2), ufs_mesh % dimension(3))) end if ! Check if the user has specified to write state points @@ -488,7 +488,7 @@ contains n_state_points = n_batches / state_point_(1) % interval allocate(statepoint_batch(n_state_points)) statepoint_batch = (/ (state_point_(1) % interval * i, i = 1, & - n_state_points) /) + n_state_points) /) else ! If neither were specified, write state point at last batch n_state_points = 1 @@ -498,7 +498,7 @@ contains ! Check if the user has specified to write binary source file if (trim(state_point_(1) % source_separate) == 'on') & - source_separate = .true. + source_separate = .true. else ! If no tag was present, by default write state point at ! last batch only @@ -629,7 +629,7 @@ contains ! Check to make sure that either material or fill was specified if (c % material == NONE .and. c % fill == NONE) then message = "Neither material nor fill was specified for cell " // & - trim(to_str(c % id)) + trim(to_str(c % id)) call fatal_error() end if @@ -643,7 +643,7 @@ contains ! Check to make sure that surfaces were specified if (.not. associated(cell_(i) % surfaces)) then message = "No surfaces specified for cell " // & - trim(to_str(c % id)) + trim(to_str(c % id)) call fatal_error() end if @@ -659,7 +659,7 @@ contains ! another universe if (c % fill == NONE) then message = "Cannot apply a rotation to cell " // trim(to_str(& - c % id)) // " because it is not filled with another universe" + c % id)) // " because it is not filled with another universe" call fatal_error() end if @@ -667,7 +667,7 @@ contains n = size(cell_(i) % rotation) if (n /= 3) then message = "Incorrect number of rotation parameters on cell " // & - to_str(c % id) + to_str(c % id) call fatal_error() end if @@ -679,13 +679,13 @@ contains ! Calculate rotation matrix based on angles given allocate(c % rotation(3,3)) c % rotation = reshape((/ & - cos(theta)*cos(psi), cos(theta)*sin(psi), -sin(theta), & - -cos(phi)*sin(psi) + sin(phi)*sin(theta)*cos(psi), & - cos(phi)*cos(psi) + sin(phi)*sin(theta)*sin(psi), & - sin(phi)*cos(theta), & - sin(phi)*sin(psi) + cos(phi)*sin(theta)*cos(psi), & - -sin(phi)*cos(psi) + cos(phi)*sin(theta)*sin(psi), & - cos(phi)*cos(theta) /), (/ 3,3 /)) + cos(theta)*cos(psi), cos(theta)*sin(psi), -sin(theta), & + -cos(phi)*sin(psi) + sin(phi)*sin(theta)*cos(psi), & + cos(phi)*cos(psi) + sin(phi)*sin(theta)*sin(psi), & + sin(phi)*cos(theta), & + sin(phi)*sin(psi) + cos(phi)*sin(theta)*cos(psi), & + -sin(phi)*cos(psi) + cos(phi)*sin(theta)*sin(psi), & + cos(phi)*cos(theta) /), (/ 3,3 /)) end if ! Translation vector @@ -694,7 +694,7 @@ contains ! another universe if (c % fill == NONE) then message = "Cannot apply a translation to cell " // trim(to_str(& - c % id)) // " because it is not filled with another universe" + c % id)) // " because it is not filled with another universe" call fatal_error() end if @@ -702,7 +702,7 @@ contains n = size(cell_(i) % translation) if (n /= 3) then message = "Incorrect number of translation parameters on cell " & - // to_str(c % id) + // to_str(c % id) call fatal_error() end if @@ -814,11 +814,11 @@ contains n = size(surface_(i) % coeffs) if (n < coeffs_reqd) then message = "Not enough coefficients specified for surface: " // & - trim(to_str(s % id)) + trim(to_str(s % id)) call fatal_error() elseif (n > coeffs_reqd) then message = "Too many coefficients specified for surface: " // & - trim(to_str(s % id)) + trim(to_str(s % id)) call fatal_error() else allocate(s % coeffs(n)) @@ -839,7 +839,7 @@ contains s % bc = BC_PERIODIC case default message = "Unknown boundary condition '" // trim(word) // & - "' specified on surface " // trim(to_str(s % id)) + "' specified on surface " // trim(to_str(s % id)) call fatal_error() end select @@ -888,7 +888,7 @@ contains ! Read lattice lower-left location if (size(lattice_(i) % dimension) /= size(lattice_(i) % lower_left)) then message = "Number of entries on must be the same as " // & - "the number of entries on ." + "the number of entries on ." call fatal_error() end if l % x0 = lattice_(i) % lower_left(1) @@ -897,7 +897,7 @@ contains ! Read lattice widths if (size(lattice_(i) % width) /= size(lattice_(i) % lower_left)) then message = "Number of entries on must be the same as " // & - "the number of entries on ." + "the number of entries on ." call fatal_error() end if l % width_x = lattice_(i) % width(1) @@ -998,7 +998,7 @@ contains sum_density = .false. if (val <= ZERO) then message = "Need to specify a positive density on material " // & - trim(to_str(mat % id)) // "." + trim(to_str(mat % id)) // "." call fatal_error() end if @@ -1015,7 +1015,7 @@ contains mat % density = 1.0e-24 * val case default message = "Unkwown units '" // trim(material_(i) % density % units) & - // "' specified on material " // trim(to_str(mat % id)) + // "' specified on material " // trim(to_str(mat % id)) call fatal_error() end select end if @@ -1026,7 +1026,7 @@ contains ! Check to ensure material has at least one nuclide if (.not. associated(material_(i) % nuclides)) then message = "No nuclides specified on material " // & - trim(to_str(mat % id)) + trim(to_str(mat % id)) call fatal_error() end if @@ -1044,7 +1044,7 @@ contains ! Check for empty name on nuclide if (len_trim(nuc % name) == 0) then message = "No name specified on nuclide in material " // & - trim(to_str(mat % id)) + trim(to_str(mat % id)) call fatal_error() end if @@ -1052,7 +1052,7 @@ contains if (len_trim(nuc % xs) == 0) then if (default_xs == '') then message = "No cross section specified for nuclide in material " & - // trim(to_str(mat % id)) + // trim(to_str(mat % id)) call fatal_error() else nuc % xs = default_xs @@ -1066,7 +1066,7 @@ contains ! Check that this nuclide is listed in the cross_sections.xml file if (.not. dict_has_key(xs_listing_dict, name)) then message = "Could not find nuclide " // trim(name) // & - " in cross_sections.xml file!" + " in cross_sections.xml file!" call fatal_error() end if @@ -1074,7 +1074,7 @@ contains n = len_trim(name) if (name(n:n) /= 'c') then message = "Cross-section table " // trim(name) // & - " is not a continuous-energy neutron table." + " is not a continuous-energy neutron table." call fatal_error() end if @@ -1099,11 +1099,11 @@ contains ! weight percents were specified if (nuc % ao == ZERO .and. nuc % wo == ZERO) then message = "No atom or weight percent specified for nuclide " // & - trim(name) + trim(name) call fatal_error() elseif (nuc % ao /= ZERO .and. nuc % wo /= ZERO) then message = "Cannot specify both atom and weight percents for a & - &nuclide: " // trim(name) + &nuclide: " // trim(name) call fatal_error() end if @@ -1118,9 +1118,9 @@ contains ! Check to make sure either all atom percents or all weight percents are ! given if (.not. (all(mat % atom_density > ZERO) .or. & - all(mat % atom_density < ZERO))) then + all(mat % atom_density < ZERO))) then message = "Cannot mix atom and weight percents in material " // & - to_str(mat % id) + to_str(mat % id) call fatal_error() end if @@ -1141,7 +1141,7 @@ contains ! Check that this nuclide is listed in the cross_sections.xml file if (.not. dict_has_key(xs_listing_dict, name)) then message = "Could not find S(a,b) table " // trim(name) // & - " in cross_sections.xml file!" + " in cross_sections.xml file!" call fatal_error() end if mat % has_sab_table = .true. @@ -1205,6 +1205,7 @@ contains logical :: file_exists ! does tallies.xml file exist? character(MAX_LINE_LEN) :: filename character(MAX_WORD_LEN) :: word + character(MAX_WORD_LEN) :: scatter_name type(ListKeyValueCI), pointer :: key_list => null() type(TallyObject), pointer :: t => null() type(StructuredMesh), pointer :: m => null() @@ -1301,7 +1302,7 @@ contains ! Check that dimensions are all greater than zero if (any(mesh_(i) % dimension <= 0)) then message = "All entries on the element for a tally mesh & - &must be positive." + &must be positive." call fatal_error() end if @@ -1311,24 +1312,24 @@ contains ! Read mesh lower-left corner location if (m % n_dimension /= size(mesh_(i) % lower_left)) then message = "Number of entries on must be the same as & - &the number of entries on ." + &the number of entries on ." call fatal_error() end if m % lower_left = mesh_(i) % lower_left ! Make sure either upper-right or width was specified if (associated(mesh_(i) % upper_right) .and. & - associated(mesh_(i) % width)) then + associated(mesh_(i) % width)) then message = "Cannot specify both and on a & - &tally mesh." + &tally mesh." call fatal_error() end if ! Make sure either upper-right or width was specified if (.not. associated(mesh_(i) % upper_right) .and. & - .not. associated(mesh_(i) % width)) then + .not. associated(mesh_(i) % width)) then message = "Must specify either and on a & - &tally mesh." + &tally mesh." call fatal_error() end if @@ -1336,7 +1337,7 @@ contains ! Check to ensure width has same dimensions if (size(mesh_(i) % width) /= size(mesh_(i) % lower_left)) then message = "Number of entries on must be the same as the & - &number of entries on ." + &number of entries on ." call fatal_error() end if @@ -1354,14 +1355,14 @@ contains ! Check to ensure width has same dimensions if (size(mesh_(i) % upper_right) /= size(mesh_(i) % lower_left)) then message = "Number of entries on must be the same as & - &the number of entries on ." + &the number of entries on ." call fatal_error() end if ! Check that upper-right is above lower-left if (any(mesh_(i) % upper_right < mesh_(i) % lower_left)) then message = "The coordinates must be greater than the & - & coordinates on a tally mesh." + & coordinates on a tally mesh." call fatal_error() end if @@ -1410,7 +1411,7 @@ contains if (size(tally_(i) % filters) > 0) then message = "Tally filters should be specified with multiple & - &elements. Did you forget to change your element?" + &elements. Did you forget to change your element?" call fatal_error() end if @@ -1442,7 +1443,7 @@ contains allocate(t % filters(j) % int_bins(n_words)) do k = 1, n_words t % filters(j) % int_bins(k) = int(str_to_int(& - tally_(i) % filter(j) % bins(k)),4) + tally_(i) % filter(j) % bins(k)),4) end do case ('cellborn') @@ -1456,7 +1457,7 @@ contains allocate(t % filters(j) % int_bins(n_words)) do k = 1, n_words t % filters(j) % int_bins(k) = int(str_to_int(& - tally_(i) % filter(j) % bins(k)),4) + tally_(i) % filter(j) % bins(k)),4) end do case ('material') @@ -1470,7 +1471,7 @@ contains allocate(t % filters(j) % int_bins(n_words)) do k = 1, n_words t % filters(j) % int_bins(k) = int(str_to_int(& - tally_(i) % filter(j) % bins(k)),4) + tally_(i) % filter(j) % bins(k)),4) end do case ('universe') @@ -1484,7 +1485,7 @@ contains allocate(t % filters(j) % int_bins(n_words)) do k = 1, n_words t % filters(j) % int_bins(k) = int(str_to_int(& - tally_(i) % filter(j) % bins(k)),4) + tally_(i) % filter(j) % bins(k)),4) end do case ('surface') @@ -1498,7 +1499,7 @@ contains allocate(t % filters(j) % int_bins(n_words)) do k = 1, n_words t % filters(j) % int_bins(k) = int(str_to_int(& - tally_(i) % filter(j) % bins(k)),4) + tally_(i) % filter(j) % bins(k)),4) end do case ('mesh') @@ -1520,7 +1521,7 @@ contains m => meshes(i_mesh) else message = "Could not find mesh " // trim(to_str(id)) // & - " specified on tally " // trim(to_str(t % id)) + " specified on tally " // trim(to_str(t % id)) call fatal_error() end if @@ -1544,7 +1545,7 @@ contains allocate(t % filters(j) % real_bins(n_words)) do k = 1, n_words t % filters(j) % real_bins(k) = str_to_real(& - tally_(i) % filter(j) % bins(k)) + tally_(i) % filter(j) % bins(k)) end do case ('energyout') @@ -1558,7 +1559,7 @@ contains allocate(t % filters(j) % real_bins(n_words)) do k = 1, n_words t % filters(j) % real_bins(k) = str_to_real(& - tally_(i) % filter(j) % bins(k)) + tally_(i) % filter(j) % bins(k)) end do end select @@ -1572,9 +1573,9 @@ contains ! Check that both cell and surface weren't specified if (t % find_filter(FILTER_CELL) > 0 .and. & - t % find_filter(FILTER_SURFACE) > 0) then + t % find_filter(FILTER_SURFACE) > 0) then message = "Cannot specify both cell and surface filters for tally " & - // trim(to_str(t % id)) + // trim(to_str(t % id)) call fatal_error() end if @@ -1593,7 +1594,7 @@ contains ! Set bins to 1, 2, 3, ..., n_nuclides_total, -1 t % nuclide_bins(1:n_nuclides_total) = & - (/ (j, j=1, n_nuclides_total) /) + (/ (j, j=1, n_nuclides_total) /) t % nuclide_bins(n_nuclides_total + 1) = -1 ! Set number of nuclide bins @@ -1621,7 +1622,7 @@ contains key_list => dict_keys(nuclide_dict) do while (associated(key_list)) if (starts_with(key_list % data % key, & - tally_(i) % nuclides(j))) then + tally_(i) % nuclides(j))) then word = key_list % data % key exit end if @@ -1633,8 +1634,8 @@ contains ! Check if no nuclide was found if (.not. associated(key_list)) then message = "Could not find the nuclide " // trim(& - tally_(i) % nuclides(j)) // " specified in tally " & - // trim(to_str(t % id)) // " in any material." + tally_(i) % nuclides(j)) // " specified in tally " & + // trim(to_str(t % id)) // " in any material." call fatal_error() end if deallocate(key_list) @@ -1647,7 +1648,7 @@ contains ! Check to make sure nuclide specified is in problem if (.not. dict_has_key(nuclide_dict, word)) then message = "The nuclide " // trim(word) // " from tally " // & - trim(to_str(t % id)) // " is not present in any material." + trim(to_str(t % id)) // " is not present in any material." call fatal_error() end if @@ -1678,26 +1679,28 @@ contains n_words = size(tally_(i) % scores) n_new = 0 do j = 1, n_words - call lower_case(tally_(i) % scores(j)) - select case (tally_(i) % scores(j)) - case ('scatter-p0') - n_new = n_new + 0 - case ('scatter-p1') - n_new = n_new + 1 - case ('scatter-p2') - n_new = n_new + 2 - case ('scatter-p3') - n_new = n_new + 3 - case ('scatter-p4') - n_new = n_new + 4 - case ('scatter-p5') - n_new = n_new + 5 - end select - end do + call lower_case(tally_(i) % scores(j)) + select case (tally_(i) % scores(j)) + case ('scatter-p0') + n_new = n_new + 0 + case ('scatter-p1') + n_new = n_new + 1 + case ('scatter-p2') + n_new = n_new + 2 + case ('scatter-p3') + n_new = n_new + 3 + case ('scatter-p4') + n_new = n_new + 4 + case ('scatter-p5') + n_new = n_new + 5 + end select + end do n_scores = n_words + n_new ! Allocate accordingly allocate(t % score_bins(n_scores)) + allocate(t % scatt_order(n_scores)) + t % scatt_order = 0 j = 0 do l = 1, n_words j = j + 1 @@ -1705,7 +1708,7 @@ contains case ('flux') ! Prohibit user from tallying flux for an individual nuclide if (.not. (t % n_nuclide_bins == 1 .and. & - t % nuclide_bins(1) == -1)) then + t % nuclide_bins(1) == -1)) then message = "Cannot tally flux for an individual nuclide." call fatal_error() end if @@ -1719,7 +1722,7 @@ contains t % score_bins(j) = SCORE_TOTAL if (t % find_filter(FILTER_ENERGYOUT) > 0) then message = "Cannot tally total reaction rate with an & - &outgoing energy filter." + &outgoing energy filter." call fatal_error() end if case ('scatter') @@ -1729,92 +1732,118 @@ contains ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG - case ('scatter-0') !does the same as 'scatter', for convenience + case ('scatter-0') + ! Does the same as 'scatter', for convenience t % score_bins(j) = SCORE_SCATTER case ('scatter-1') - t % score_bins(j) = SCORE_SCATTER_1 + t % score_bins(j) = SCORE_SCATTER_N + t % scatt_order(j) = 1 ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG case ('scatter-2') - t % score_bins(j) = SCORE_SCATTER_2 + t % score_bins(j) = SCORE_SCATTER_N + t % scatt_order(j) = 2 ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG case ('scatter-3') - t % score_bins(j) = SCORE_SCATTER_3 + t % score_bins(j) = SCORE_SCATTER_N + t % scatt_order(j) = 3 ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG case ('scatter-4') - t % score_bins(j) = SCORE_SCATTER_4 + t % score_bins(j) = SCORE_SCATTER_N + t % scatt_order(j) = 4 ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG case ('scatter-5') - t % score_bins(j) = SCORE_SCATTER_5 + t % score_bins(j) = SCORE_SCATTER_N + t % scatt_order(j) = 5 ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG case ('scatter-p0') - ! Same as a scatter-0 or scatter - t % score_bins(j) = SCORE_SCATTER - case ('scatter-p1') - ! Setup each of P0 and P1 - t % estimator = ESTIMATOR_ANALOG - t % score_bins(j) = SCORE_SCATTER - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_1 - - case ('scatter-p2') - ! Setup each of P0:P2 - t % estimator = ESTIMATOR_ANALOG - t % score_bins(j) = SCORE_SCATTER - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_1 - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_2 - - case ('scatter-p3') - ! Setup each of P0:P3 - t % estimator = ESTIMATOR_ANALOG - t % score_bins(j) = SCORE_SCATTER - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_1 - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_2 - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_3 - - case ('scatter-p4') - ! Setup each of P0:P4 - t % estimator = ESTIMATOR_ANALOG - t % score_bins(j) = SCORE_SCATTER - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_1 - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_2 - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_3 - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_4 - - case ('scatter-p5') - ! Setup each of P0:P5 - t % estimator = ESTIMATOR_ANALOG - t % score_bins(j) = SCORE_SCATTER - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_1 - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_2 - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_3 - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_4 - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_5 - + ! Same as a scatter-0 or scatter + t % score_bins(j) = SCORE_SCATTER + case ('scatter-p1') + ! Setup each of P0 and P1 + t % estimator = ESTIMATOR_ANALOG + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 1 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 1 + + case ('scatter-p2') + ! Setup each of P0:P2 + t % estimator = ESTIMATOR_ANALOG + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 2 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 2 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 2 + + case ('scatter-p3') + ! Setup each of P0:P3 + t % estimator = ESTIMATOR_ANALOG + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 3 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 3 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 3 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 3 + + case ('scatter-p4') + ! Setup each of P0:P4 + t % estimator = ESTIMATOR_ANALOG + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 4 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 4 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 4 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 4 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 4 + + case ('scatter-p5') + ! Setup each of P0:P5 + t % estimator = ESTIMATOR_ANALOG + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 5 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 5 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 5 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 5 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 5 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 5 + case('transport') t % score_bins(j) = SCORE_TRANSPORT @@ -1849,14 +1878,14 @@ contains t % score_bins(j) = SCORE_ABSORPTION if (t % find_filter(FILTER_ENERGYOUT) > 0) then message = "Cannot tally absorption rate with an outgoing & - &energy filter." + &energy filter." call fatal_error() end if case ('fission') t % score_bins(j) = SCORE_FISSION if (t % find_filter(FILTER_ENERGYOUT) > 0) then message = "Cannot tally fission rate with an outgoing & - &energy filter." + &energy filter." call fatal_error() end if case ('nu-fission') @@ -1869,8 +1898,8 @@ contains ! for this tally if (n_words > 1) then message = "Cannot tally other scoring functions in the same & - &tally as surface currents. Separate other scoring & - &functions into a distinct tally." + &tally as surface currents. Separate other scoring & + &functions into a distinct tally." call fatal_error() end if @@ -1902,13 +1931,13 @@ contains t % filters(t % n_filters) % type = FILTER_SURFACE t % filters(t % n_filters) % n_bins = 2 * m % n_dimension allocate(t % filters(t % n_filters) % int_bins(& - 2 * m % n_dimension)) + 2 * m % n_dimension)) if (m % n_dimension == 2) then t % filters(t % n_filters) % int_bins = (/ IN_RIGHT, & - OUT_RIGHT, IN_FRONT, OUT_FRONT /) + OUT_RIGHT, IN_FRONT, OUT_FRONT /) elseif (m % n_dimension == 3) then t % filters(t % n_filters) % int_bins = (/ IN_RIGHT, & - OUT_RIGHT, IN_FRONT, OUT_FRONT, IN_TOP, OUT_TOP /) + OUT_RIGHT, IN_FRONT, OUT_FRONT, IN_TOP, OUT_TOP /) end if t % find_filter(FILTER_SURFACE) = t % n_filters @@ -1917,14 +1946,15 @@ contains case default message = "Unknown scoring function: " // & - trim(tally_(i) % scores(l)) + trim(tally_(i) % scores(l)) call fatal_error() end select end do t % n_score_bins = n_scores + t % n_nonPN_score_bins = n_words else message = "No specified on tally " // trim(to_str(t % id)) & - // "." + // "." call fatal_error() end if @@ -1942,7 +1972,7 @@ contains ! tally needs post-collision information if (t % estimator == ESTIMATOR_ANALOG) then message = "Cannot use track-length estimator for tally " & - // to_str(t % id) + // to_str(t % id) call fatal_error() end if @@ -1951,7 +1981,7 @@ contains case default message = "Invalid estimator '" // trim(tally_(i) % estimator) & - // "' on tally " // to_str(t % id) + // "' on tally " // to_str(t % id) call fatal_error() end select end if @@ -2055,7 +2085,7 @@ contains ! Set output file path pl % path_plot = trim(path_input) // trim(to_str(pl % id)) // & - "_" // trim(plot_(i) % filename) // ".ppm" + "_" // trim(plot_(i) % filename) // ".ppm" ! Copy plot pixel size if (size(plot_(i) % pixels) == 2) then @@ -2071,7 +2101,7 @@ contains pl % not_found % rgb = plot_(i) % background else message = "Bad background RGB " & - // "in plot " // trim(to_str(pl % id)) + // "in plot " // trim(to_str(pl % id)) call fatal_error() end if else @@ -2086,7 +2116,7 @@ contains ! pl % type = PLOT_TYPE_POINTS case default message = "Unsupported plot type '" // plot_(i) % type & - // "' in plot " // trim(to_str(pl % id)) + // "' in plot " // trim(to_str(pl % id)) call fatal_error() end select @@ -2100,7 +2130,7 @@ contains pl % basis = PLOT_BASIS_YZ case default message = "Unsupported plot basis '" // plot_(i) % basis & - // "' in plot " // trim(to_str(pl % id)) + // "' in plot " // trim(to_str(pl % id)) call fatal_error() end select @@ -2109,7 +2139,7 @@ contains pl % origin = plot_(i) % origin else message = "Origin must be length 3 " & - // "in plot " // trim(to_str(pl % id)) + // "in plot " // trim(to_str(pl % id)) call fatal_error() end if @@ -2121,7 +2151,7 @@ contains pl % width(2) = plot_(i) % width(2) else message = "Bad plot width " & - // "in plot " // trim(to_str(pl % id)) + // "in plot " // trim(to_str(pl % id)) call fatal_error() end if @@ -2149,7 +2179,7 @@ contains case default message = "Unsupported plot color type '" // plot_(i) % color & - // "' in plot " // trim(to_str(pl % id)) + // "' in plot " // trim(to_str(pl % id)) call fatal_error() end select @@ -2159,7 +2189,7 @@ contains do j = 1, n_cols if (size(plot_(i) % col_spec_(j) % rgb) /= 3) then message = "Bad RGB " & - // "in plot " // trim(to_str(pl % id)) + // "in plot " // trim(to_str(pl % id)) call fatal_error() end if @@ -2171,7 +2201,7 @@ contains pl % colors(col_id) % rgb = plot_(i) % col_spec_(j) % rgb else message = "Could not find cell " // trim(to_str(col_id)) // & - " specified in plot " // trim(to_str(pl % id)) + " specified in plot " // trim(to_str(pl % id)) call fatal_error() end if @@ -2181,7 +2211,7 @@ contains pl % colors(col_id) % rgb = plot_(i) % col_spec_(j) % rgb else message = "Could not find material " // trim(to_str(col_id)) // & - " specified in plot " // trim(to_str(pl % id)) + " specified in plot " // trim(to_str(pl % id)) call fatal_error() end if @@ -2193,7 +2223,7 @@ contains if (associated(plot_(i) % mask_)) then if (size(plot_(i) % mask_) > 1) then message = "Mutliple masks" // & - " specified in plot " // trim(to_str(pl % id)) + " specified in plot " // trim(to_str(pl % id)) call fatal_error() else if (.not. size(plot_(i) % mask_) == 0) then do j=1,size(pl % colors) @@ -2228,10 +2258,10 @@ contains ! Check if cross_sections.xml exists inquire(FILE=path_cross_sections, EXIST=file_exists) if (.not. file_exists) then - ! Could not find cross_sections.xml file - message = "Cross sections XML file '" // trim(path_cross_sections) // & + ! Could not find cross_sections.xml file + message = "Cross sections XML file '" // trim(path_cross_sections) // & "' does not exist!" - call fatal_error() + call fatal_error() end if message = "Reading cross sections XML file..." @@ -2247,25 +2277,25 @@ contains call read_xml_file_cross_sections_t(path_cross_sections) if (len_trim(directory_) > 0) then - ! Copy directory information if present - directory = trim(directory_) + ! Copy directory information if present + directory = trim(directory_) else - ! If no directory is listed in cross_sections.xml, by default select the - ! directory in which the cross_sections.xml file resides - i = index(path_cross_sections, "/", BACK=.true.) - directory = path_cross_sections(1:i) + ! If no directory is listed in cross_sections.xml, by default select the + ! directory in which the cross_sections.xml file resides + i = index(path_cross_sections, "/", BACK=.true.) + directory = path_cross_sections(1:i) end if ! determine whether binary/ascii if (filetype_ == 'ascii') then - filetype = ASCII + filetype = ASCII elseif (filetype_ == 'binary') then - filetype = BINARY + filetype = BINARY elseif (len_trim(filetype_) == 0) then - filetype = ASCII + filetype = ASCII else - message = "Unknown filetype in cross_sections.xml: " // trim(filetype_) - call fatal_error() + message = "Unknown filetype in cross_sections.xml: " // trim(filetype_) + call fatal_error() end if ! copy default record length and entries for binary files @@ -2274,57 +2304,57 @@ contains ! Allocate xs_listings array if (.not. associated(ace_tables_)) then - message = "No ACE table listings present in cross_sections.xml file!" - call fatal_error() + message = "No ACE table listings present in cross_sections.xml file!" + call fatal_error() else - n_listings = size(ace_tables_) - allocate(xs_listings(n_listings)) + n_listings = size(ace_tables_) + allocate(xs_listings(n_listings)) end if do i = 1, n_listings - listing => xs_listings(i) + listing => xs_listings(i) - ! copy a number of attributes - listing % name = trim(ace_tables_(i) % name) - listing % alias = trim(ace_tables_(i) % alias) - listing % zaid = ace_tables_(i) % zaid - listing % awr = ace_tables_(i) % awr - listing % kT = ace_tables_(i) % temperature - listing % location = ace_tables_(i) % location + ! copy a number of attributes + listing % name = trim(ace_tables_(i) % name) + listing % alias = trim(ace_tables_(i) % alias) + listing % zaid = ace_tables_(i) % zaid + listing % awr = ace_tables_(i) % awr + listing % kT = ace_tables_(i) % temperature + listing % location = ace_tables_(i) % location - ! determine type of cross section - if (ends_with(listing % name, 'c')) then + ! determine type of cross section + if (ends_with(listing % name, 'c')) then listing % type = ACE_NEUTRON - elseif (ends_with(listing % name, 't')) then + elseif (ends_with(listing % name, 't')) then listing % type = ACE_THERMAL - end if + end if - ! set filetype, record length, and number of entries - listing % filetype = filetype - listing % recl = recl - listing % entries = entries + ! set filetype, record length, and number of entries + listing % filetype = filetype + listing % recl = recl + listing % entries = entries - ! determine metastable state - if (ace_tables_(i) % metastable == 0) then + ! determine metastable state + if (ace_tables_(i) % metastable == 0) then listing % metastable = .false. - else + else listing % metastable = .true. - end if + end if - ! determine path of cross section table - if (starts_with(ace_tables_(i) % path, '/')) then + ! determine path of cross section table + if (starts_with(ace_tables_(i) % path, '/')) then listing % path = ace_tables_(i) % path - else + else if (ends_with(directory,'/')) then - listing % path = trim(directory) // trim(ace_tables_(i) % path) + listing % path = trim(directory) // trim(ace_tables_(i) % path) else - listing % path = trim(directory) // '/' // trim(ace_tables_(i) % path) + listing % path = trim(directory) // '/' // trim(ace_tables_(i) % path) end if - end if + end if - ! create dictionary entry for both name and alias - call dict_add_key(xs_listing_dict, listing % name, i) - call dict_add_key(xs_listing_dict, listing % alias, i) + ! create dictionary entry for both name and alias + call dict_add_key(xs_listing_dict, listing % name, i) + call dict_add_key(xs_listing_dict, listing % alias, i) end do end subroutine read_cross_sections_xml diff --git a/src/output.F90 b/src/output.F90 index 04cfcb8a9a..2771913a43 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -32,35 +32,35 @@ contains subroutine title() write(UNIT=OUTPUT_UNIT, FMT='(/11(A/))') & - ' .d88888b. 888b d888 .d8888b.', & - ' d88P" "Y88b 8888b d8888 d88P Y88b', & - ' 888 888 88888b.d88888 888 888', & - ' 888 888 88888b. .d88b. 88888b. 888Y88888P888 888 ', & - ' 888 888 888 "88b d8P Y8b 888 "88b 888 Y888P 888 888 ', & - ' 888 888 888 888 88888888 888 888 888 Y8P 888 888 888', & - ' Y88b. .d88P 888 d88P Y8b. 888 888 888 " 888 Y88b d88P', & - ' "Y88888P" 88888P" "Y8888 888 888 888 888 "Y8888P"', & - '__________________888______________________________________________________', & - ' 888', & - ' 888' + ' .d88888b. 888b d888 .d8888b.', & + ' d88P" "Y88b 8888b d8888 d88P Y88b', & + ' 888 888 88888b.d88888 888 888', & + ' 888 888 88888b. .d88b. 88888b. 888Y88888P888 888 ', & + ' 888 888 888 "88b d8P Y8b 888 "88b 888 Y888P 888 888 ', & + ' 888 888 888 888 88888888 888 888 888 Y8P 888 888 888', & + ' Y88b. .d88P 888 d88P Y8b. 888 888 888 " 888 Y88b d88P', & + ' "Y88888P" 88888P" "Y8888 888 888 888 888 "Y8888P"', & + '__________________888______________________________________________________', & + ' 888', & + ' 888' ! Write version information write(UNIT=OUTPUT_UNIT, FMT=*) & - ' Developed At: Massachusetts Institute of Technology' + ' Developed At: Massachusetts Institute of Technology' write(UNIT=OUTPUT_UNIT, FMT='(6X,"Version:",7X,I1,".",I1,".",I1)') & - VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE + VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE #ifdef GIT_SHA1 write(UNIT=OUTPUT_UNIT, FMT='(6X,"Git SHA1:",6X,A)') GIT_SHA1 #endif ! Write the date and time write(UNIT=OUTPUT_UNIT, FMT='(6X,"Date/Time:",5X,A)') & - time_stamp() + time_stamp() #ifdef MPI ! Write number of processors write(UNIT=OUTPUT_UNIT, FMT='(6X,"MPI Processes:",1X,A)') & - trim(to_str(n_procs)) + trim(to_str(n_procs)) #endif end subroutine title @@ -77,7 +77,7 @@ contains call date_and_time(DATE=date_, TIME=time_) current_time = date_(1:4) // "-" // date_(5:6) // "-" // date_(7:8) // & - " " // time_(1:2) // ":" // time_(3:4) // ":" // time_(5:6) + " " // time_(1:2) // ":" // time_(3:4) // ":" // time_(5:6) end function time_stamp @@ -125,13 +125,13 @@ contains select case (header_level) case (1) write(UNIT=unit_, FMT='(/3(1X,A/))') repeat('=', 75), & - repeat('=', n) // '> ' // trim(line) // ' <' // & - repeat('=', m), repeat('=', 75) + repeat('=', n) // '> ' // trim(line) // ' <' // & + repeat('=', m), repeat('=', 75) case (2) write(UNIT=unit_, FMT='(/2(1X,A/))') trim(line), repeat('-', 75) case (3) write(UNIT=unit_, FMT='(/1X,A/)') repeat('=', n) // '> ' // & - trim(line) // ' <' // repeat('=', m) + trim(line) // ' <' // repeat('=', m) end select end subroutine header @@ -145,11 +145,11 @@ contains if (master) then write(UNIT=OUTPUT_UNIT, FMT='(1X,A,1X,I1,".",I1,".",I1)') & - "OpenMC version", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE + "OpenMC version", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE write(UNIT=OUTPUT_UNIT, FMT=*) "Copyright (c) 2011-2012 & - &Massachusetts Institute of Technology" + &Massachusetts Institute of Technology" write(UNIT=OUTPUT_UNIT, FMT=*) "MIT/X license at & - &" + &" end if end subroutine print_version @@ -207,7 +207,7 @@ contains else ! Determine last space in current line i_end = i_start + index(message(i_start+1:i_start+line_wrap), & - ' ', BACK=.true.) + ' ', BACK=.true.) ! Write up to last space write(ou, fmt='(1X,A)') message(i_start+1:i_end-1) @@ -270,8 +270,8 @@ contains l => lattices(coord % lattice) write(ou,*) ' Lattice = ' // trim(to_str(l % id)) write(ou,*) ' Lattice position = (' // trim(to_str(& - p % coord % lattice_x)) // ',' // trim(to_str(& - p % coord % lattice_y)) // ')' + p % coord % lattice_x)) // ',' // trim(to_str(& + p % coord % lattice_y)) // ')' end if ! Print local coordinates @@ -390,7 +390,7 @@ contains case default index_surf = abs(c % surfaces(i)) string = trim(string) // ' ' // to_str(sign(& - surfaces(index_surf) % id, c % surfaces(i))) + surfaces(index_surf) % id, c % surfaces(i))) end select end do write(unit_,*) ' Surface Specification:' // trim(string) @@ -595,7 +595,7 @@ contains ! Write total atom density in atom/b-cm write(unit_,*) ' Atom Density = ' // trim(to_str(mat % density)) & - // ' atom/b-cm' + // ' atom/b-cm' ! Write atom density for each nuclide in material write(unit_,*) ' Nuclides:' @@ -603,7 +603,7 @@ contains nuc => nuclides(mat % nuclide(i)) density = mat % atom_density(i) string = ' ' // trim(nuc % name) // ' = ' // & - trim(to_str(density)) // ' atom/b-cm' + trim(to_str(density)) // ' atom/b-cm' write(unit_,*) trim(string) end do @@ -628,7 +628,9 @@ contains integer :: j ! index in filters array integer :: id ! user-specified id integer :: unit_ ! unit to write to + integer :: n ! scattering order to include in name character(MAX_LINE_LEN) :: string + character(MAX_WORD_LEN) :: pn_string type(Cell), pointer :: c => null() type(Surface), pointer :: s => null() type(Universe), pointer :: u => null() @@ -740,7 +742,7 @@ contains string = "" do i = 1, t % filters(j) % n_bins + 1 string = trim(string) // ' ' // trim(to_str(& - t % filters(j) % real_bins(i))) + t % filters(j) % real_bins(i))) end do write(unit_,*) ' Incoming Energy Bins:' // trim(string) end if @@ -751,7 +753,7 @@ contains string = "" do i = 1, t % filters(j) % n_bins + 1 string = trim(string) // ' ' // trim(to_str(& - t % filters(j) % real_bins(i))) + t % filters(j) % real_bins(i))) end do write(unit_,*) ' Outgoing Energy Bins:' // trim(string) end if @@ -763,18 +765,20 @@ contains write(unit_,fmt='(A)',advance='no') ' total' else write(unit_,fmt='(A)',advance='no') ' ' // trim(adjustl(& - nuclides(t % nuclide_bins(i)) % name)) + nuclides(t % nuclide_bins(i)) % name)) end if if (mod(i,4) == 0 .and. i /= t % n_nuclide_bins) & - write(unit_,'(/18X)',advance='no') + write(unit_,'(/18X)',advance='no') end do write(unit_,*) ! Write score bins - string = "" - do i = 1, t % n_score_bins - select case (t % score_bins(i)) + string = "" + j = 0 + do i = 1, t % n_nonPN_score_bins + j = j + 1 + select case (t % score_bins(j)) case (SCORE_FLUX) string = trim(string) // ' flux' case (SCORE_TOTAL) @@ -783,16 +787,17 @@ contains string = trim(string) // ' scatter' case (SCORE_NU_SCATTER) string = trim(string) // ' nu-scatter' - case (SCORE_SCATTER_1) - string = trim(string) // ' scatter-1' - case (SCORE_SCATTER_2) - string = trim(string) // ' scatter-2' - case (SCORE_SCATTER_3) - string = trim(string) // ' scatter-3' - case (SCORE_SCATTER_4) - string = trim(string) // ' scatter-4' - case (SCORE_SCATTER_5) - string = trim(string) // ' scatter-5' + case (SCORE_SCATTER_N) + pn_string = ' scatter-' // trim(to_str(t % scatt_order(j))) + string = trim(string) // pn_string + case (SCORE_SCATTER_PN) + pn_string = ' scatter' + string = trim(string) // pn_string + do n = 1, t % scatt_order(j) + pn_string = ' scatter-' // trim(to_str(n)) + string = trim(string) // pn_string + end do + j = j + n - 1 case (SCORE_TRANSPORT) string = trim(string) // ' transport' case (SCORE_DIFFUSION) @@ -933,8 +938,8 @@ contains end if write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,A4,1X,I6,1X,I11,1X,I11)') & - reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, & - law(1:4), rxn % threshold, size_angle, size_energy + reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, & + law(1:4), rxn % threshold, size_angle, size_energy ! Accumulate data size size_xs = size_xs + (nuc % n_grid - rxn%threshold + 1) * 8 @@ -971,11 +976,11 @@ contains write(unit_,*) ' Memory Requirements' write(unit_,*) ' Cross sections = ' // trim(to_str(size_xs)) // ' bytes' write(unit_,*) ' Secondary angle distributions = ' // & - trim(to_str(size_angle_total)) // ' bytes' + trim(to_str(size_angle_total)) // ' bytes' write(unit_,*) ' Secondary energy distributions = ' // & - trim(to_str(size_energy_total)) // ' bytes' + trim(to_str(size_energy_total)) // ' bytes' write(unit_,*) ' Probability Tables = ' // & - trim(to_str(size_urr)) // ' bytes' + trim(to_str(size_urr)) // ' bytes' write(unit_,*) ' Total = ' // trim(to_str(size_total)) // ' bytes' ! Blank line at end of nuclide @@ -1011,28 +1016,28 @@ contains ! Inelastic data write(unit_,*) ' # of Incoming Energies (Inelastic) = ' // & - trim(to_str(sab % n_inelastic_e_in)) + trim(to_str(sab % n_inelastic_e_in)) write(unit_,*) ' # of Outgoing Energies (Inelastic) = ' // & - trim(to_str(sab % n_inelastic_e_out)) + trim(to_str(sab % n_inelastic_e_out)) write(unit_,*) ' # of Outgoing Angles (Inelastic) = ' // & - trim(to_str(sab % n_inelastic_mu)) + trim(to_str(sab % n_inelastic_mu)) write(unit_,*) ' Threshold for Inelastic = ' // & - trim(to_str(sab % threshold_inelastic)) + trim(to_str(sab % threshold_inelastic)) ! Elastic data if (sab % n_elastic_e_in > 0) then write(unit_,*) ' # of Incoming Energies (Elastic) = ' // & - trim(to_str(sab % n_elastic_e_in)) + trim(to_str(sab % n_elastic_e_in)) write(unit_,*) ' # of Outgoing Angles (Elastic) = ' // & - trim(to_str(sab % n_elastic_mu)) + trim(to_str(sab % n_elastic_mu)) write(unit_,*) ' Threshold for Elastic = ' // & - trim(to_str(sab % threshold_elastic)) + trim(to_str(sab % threshold_elastic)) end if ! Determine memory used by S(a,b) table and write out size_sab = 8 * (sab % n_inelastic_e_in * (2 + sab % n_inelastic_e_out * & - (1 + sab % n_inelastic_mu)) + sab % n_elastic_e_in * & - (2 + sab % n_elastic_mu)) + (1 + sab % n_inelastic_mu)) + sab % n_elastic_e_in * & + (2 + sab % n_elastic_mu)) write(unit_,*) ' Memory Used = ' // trim(to_str(size_sab)) // ' bytes' ! Blank line at end @@ -1060,19 +1065,19 @@ contains call header("OpenMC Monte Carlo Code", unit=UNIT_SUMMARY, level=1) write(UNIT=UNIT_SUMMARY, FMT=*) & - "Copyright: 2011-2012 Massachusetts Institute of Technology" + "Copyright: 2011-2012 Massachusetts Institute of Technology" write(UNIT=UNIT_SUMMARY, FMT='(1X,A,7X,2(I1,"."),I1)') & - "Version:", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE + "Version:", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE #ifdef GIT_SHA1 write(UNIT=UNIT_SUMMARY, FMT='(1X,"Git SHA1:",6X,A)') GIT_SHA1 #endif write(UNIT=UNIT_SUMMARY, FMT='(1X,"Date/Time:",5X,A)') & - time_stamp() + time_stamp() ! Write information on number of processors #ifdef MPI write(UNIT=UNIT_SUMMARY, FMT='(1X,"MPI Processes:",1X,A)') & - trim(to_str(n_procs)) + trim(to_str(n_procs)) #endif ! Display problem summary @@ -1117,7 +1122,7 @@ contains call header("UNIONIZED ENERGY GRID", unit=UNIT_SUMMARY) write(UNIT_SUMMARY,*) "Points on energy grid: " // trim(to_str(n_grid)) write(UNIT_SUMMARY,*) "Extra storage required: " // trim(to_str(& - n_grid*n_nuclides_total*4)) // " bytes" + n_grid*n_nuclides_total*4)) // " bytes" ! print summary of variance reduction call header("VARIANCE REDUCTION", unit=UNIT_SUMMARY) @@ -1227,16 +1232,16 @@ contains ! write out information batch and option independent output write(UNIT=OUTPUT_UNIT, FMT='(2X,I5)', ADVANCE='NO') current_batch write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') & - k_batch(current_batch) + k_batch(current_batch) ! write out entropy info if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & - entropy(current_batch) + entropy(current_batch) ! write out accumulated k-effective if after first active batch if (current_batch > n_inactive + 1) then write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5," +/-",F8.5)', ADVANCE='NO') & - keff, keff_std + keff, keff_std else write(UNIT=OUTPUT_UNIT, FMT='(23X)', ADVANCE='NO') end if @@ -1244,11 +1249,11 @@ contains ! write out cmfd keff if it is active if (cmfd_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & - cmfd % keff + cmfd % keff ! write out cmfd entopy if (cmfd_on .and. entropy_on) write(UNIT=OUTPUT_UNIT, & - FMT='(3X, F8.5)', ADVANCE='NO') cmfd % entropy + FMT='(3X, F8.5)', ADVANCE='NO') cmfd % entropy ! next line write(UNIT=OUTPUT_UNIT, FMT=*) @@ -1275,18 +1280,18 @@ contains ! Write plotting origin write(ou,100) "Origin:", trim(to_str(pl % origin(1))) // & - " " // trim(to_str(pl % origin(2))) // " " // & - trim(to_str(pl % origin(3))) + " " // trim(to_str(pl % origin(2))) // " " // & + trim(to_str(pl % origin(3))) ! Write plotting width if (pl % type == PLOT_TYPE_SLICE) then write(ou,100) "Width:", trim(to_str(pl % width(1))) // & - " " // trim(to_str(pl % width(2))) + " " // trim(to_str(pl % width(2))) write(ou,100) "Coloring:", trim(to_str(pl % color_by)) write(ou,100) "Basis:", trim(to_str(pl % basis)) write(ou,100) "Pixels:", trim(to_str(pl % pixels(1))) // " " // & - trim(to_str(pl % pixels(2))) + trim(to_str(pl % pixels(2))) end if write(ou,*) @@ -1319,7 +1324,7 @@ contains write(ou,100) " Reading cross sections", time_read_xs % elapsed write(ou,100) " Unionizing energy grid", time_unionize % elapsed write(ou,100) "Total time in simulation", time_inactive % elapsed + & - time_active % elapsed + time_active % elapsed write(ou,100) " Time in transport only", time_transport % elapsed if(cmfd_run) write(ou,100) "Total CMFD time", time_cmfd % elapsed write(ou,100) " Time in inactive batches", time_inactive % elapsed @@ -1333,14 +1338,14 @@ contains if (restart_run) then total_particles = n_particles * (n_batches - & - restart_batch) * gen_per_batch + restart_batch) * gen_per_batch else total_particles = n_particles * n_batches * gen_per_batch end if ! display calculate rate speed = real(total_particles) / (time_inactive % elapsed + & - time_active % elapsed) + time_active % elapsed) string = to_str(speed) write(ou,101) "Calculation Rate", trim(string) @@ -1358,13 +1363,13 @@ contains ! write global tallies write(ou,102) "k-effective (Analog)", global_tallies(K_ANALOG) % sum, & - global_tallies(K_ANALOG) % sum_sq + global_tallies(K_ANALOG) % sum_sq write(ou,102) "k-effective (Collision)", global_tallies(K_COLLISION) % sum, & - global_tallies(K_COLLISION) % sum_sq + global_tallies(K_COLLISION) % sum_sq write(ou,102) "k-effective (Track-length)", global_tallies(K_TRACKLENGTH) % sum, & - global_tallies(K_TRACKLENGTH) % sum_sq + global_tallies(K_TRACKLENGTH) % sum_sq write(ou,102) "Leakage Fraction", global_tallies(LEAKAGE) % sum, & - global_tallies(LEAKAGE) % sum_sq + global_tallies(LEAKAGE) % sum_sq write(ou,*) ! format for write statements @@ -1385,17 +1390,21 @@ contains integer :: j ! level in tally hierarchy integer :: k ! loop index for scoring bins integer :: n ! loop index for nuclides + integer :: l ! loop index for nonPN scores integer :: type ! type of tally filter integer :: indent ! number of spaces to preceed output integer :: filter_index ! index in scores array for filters integer :: score_index ! scoring bin index integer :: i_nuclide ! index in nuclides array integer :: i_listing ! index in xs_listings array + integer :: nOrder ! loop index for scattering orders real(8) :: t_value ! t-values for confidence intervals real(8) :: alpha ! significance level for CI 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 + character(27) :: score_names(N_SCORE_TYPES) ! names of scoring function + character(27) :: score_name ! names of scoring function + ! to be applied at write-time type(TallyObject), pointer :: t ! Skip if there are no tallies @@ -1412,25 +1421,22 @@ contains filter_name(FILTER_ENERGYOUT) = "Outgoing Energy" ! 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_SCATTER_4)) = "Fourth Scattering Moment" - score_name(abs(SCORE_SCATTER_5)) = "Fifth Scattering Moment" - score_name(abs(SCORE_TRANSPORT)) = "Transport Rate" - score_name(abs(SCORE_DIFFUSION)) = "Diffusion Coefficient" - 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" - score_name(abs(SCORE_EVENTS)) = "Events" + score_names(abs(SCORE_FLUX)) = "Flux" + score_names(abs(SCORE_TOTAL)) = "Total Reaction Rate" + score_names(abs(SCORE_SCATTER)) = "Scattering Rate" + score_names(abs(SCORE_NU_SCATTER)) = "Scattering Production Rate" + score_names(abs(SCORE_SCATTER_N)) = "" + score_names(abs(SCORE_SCATTER_PN)) = "" + score_names(abs(SCORE_TRANSPORT)) = "Transport Rate" + score_names(abs(SCORE_DIFFUSION)) = "Diffusion Coefficient" + score_names(abs(SCORE_N_1N)) = "(n,1n) Rate" + score_names(abs(SCORE_N_2N)) = "(n,2n) Rate" + score_names(abs(SCORE_N_3N)) = "(n,3n) Rate" + score_names(abs(SCORE_N_4N)) = "(n,4n) Rate" + score_names(abs(SCORE_ABSORPTION)) = "Absorption Rate" + score_names(abs(SCORE_FISSION)) = "Fission Rate" + score_names(abs(SCORE_NU_FISSION)) = "Nu-Fission Rate" + score_names(abs(SCORE_EVENTS)) = "Events" ! Create filename for tally output if (run_mode == MODE_TALLIES) then @@ -1468,10 +1474,10 @@ contains ! Write header block if (t % label == "") then call header("TALLY " // trim(to_str(t % id)), unit=UNIT_TALLY, & - level=3) + level=3) else call header("TALLY " // trim(to_str(t % id)) // ": " & - // trim(t % label), unit=UNIT_TALLY, level=3) + // trim(t % label), unit=UNIT_TALLY, level=3) endif ! Handle surface current tallies separately @@ -1520,7 +1526,7 @@ contains ! Print current filter information type = t % filters(j) % type write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A)') repeat(" ", indent), & - trim(filter_name(type)), trim(get_label(t, j)) + trim(filter_name(type)), trim(get_label(t, j)) indent = indent + 2 j = j + 1 end if @@ -1531,7 +1537,7 @@ contains if (t % n_filters > 0) then type = t % filters(j) % type write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A)') repeat(" ", indent), & - trim(filter_name(type)), trim(get_label(t, j)) + trim(filter_name(type)), trim(get_label(t, j)) end if ! Determine scoring index for this bin combination -- note that unlike @@ -1543,7 +1549,7 @@ contains else filter_index = 1 end if - + ! Write scores for this filter bin combination score_index = 0 if (t % n_filters > 0) indent = indent + 2 @@ -1552,20 +1558,52 @@ contains i_nuclide = t % nuclide_bins(n) if (i_nuclide == -1) then write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A)') repeat(" ", indent), & - "Total Material" + "Total Material" else i_listing = nuclides(i_nuclide) % listing write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A)') repeat(" ", indent), & - trim(xs_listings(i_listing) % alias) + trim(xs_listings(i_listing) % alias) end if indent = indent + 2 - do k = 1, t % n_score_bins + k = 0 + do l = 1, t % n_nonPN_score_bins + k = k + 1 score_index = score_index + 1 - write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & - repeat(" ", indent), score_name(abs(t % score_bins(k))), & - to_str(t % scores(score_index,filter_index) % sum), & - trim(to_str(t % scores(score_index,filter_index) % sum_sq)) + if (t % score_bins(k) == SCORE_SCATTER_N) then + if (t % scatt_order(k) == 0) then + score_name = "Scattering Rate" + else + score_name = 'Order-' // trim(to_str(t % scatt_order(k))) // & + ' Scattering Moment' + end if + write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & + repeat(" ", indent), score_name, & + to_str(t % scores(score_index,filter_index) % sum), & + trim(to_str(t % scores(score_index,filter_index) % sum_sq)) + else if (t % score_bins(k) == SCORE_SCATTER_PN) then + score_name = "Scattering Rate" + write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & + repeat(" ", indent), score_name, & + to_str(t % scores(score_index,filter_index) % sum), & + trim(to_str(t % scores(score_index,filter_index) % sum_sq)) + do nOrder = 1, t % scatt_order(k) + score_index = score_index + 1 + score_name = 'Order-' // trim(to_str(nOrder)) // & + ' Scattering Moment' + write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & + repeat(" ", indent), score_name, & + to_str(t % scores(score_index,filter_index) % sum), & + trim(to_str(t % scores(score_index,filter_index) % sum_sq)) + end do + k = k + nOrder - 1 + else + score_name = score_names(abs(t % score_bins(k))) + write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & + repeat(" ", indent), score_name, & + to_str(t % scores(score_index,filter_index) % sum), & + trim(to_str(t % scores(score_index,filter_index) % sum_sq)) + end if end do indent = indent - 2 @@ -1642,110 +1680,110 @@ contains ! Write incoming energy bin write(UNIT=UNIT_TALLY, FMT='(3X,A,1X,A)') & - "Incoming Energy", trim(get_label(t, i_filter_ein)) + "Incoming Energy", trim(get_label(t, i_filter_ein)) end if ! Left Surface t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i-1, j, k /) + 1, .true.) + mesh_indices_to_bin(m, (/ i-1, j, k /) + 1, .true.) t % matching_bins(i_filter_surf) = IN_RIGHT filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Outgoing Current to Left", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Outgoing Current to Left", & + to_str(t % scores(1,filter_index) % sum), & + trim(to_str(t % scores(1,filter_index) % sum_sq)) t % matching_bins(i_filter_surf) = OUT_RIGHT filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Incoming Current from Left", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Incoming Current from Left", & + to_str(t % scores(1,filter_index) % sum), & + trim(to_str(t % scores(1,filter_index) % sum_sq)) ! Right Surface t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) + mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) t % matching_bins(i_filter_surf) = IN_RIGHT filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Incoming Current from Right", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Incoming Current from Right", & + to_str(t % scores(1,filter_index) % sum), & + trim(to_str(t % scores(1,filter_index) % sum_sq)) t % matching_bins(i_filter_surf) = OUT_RIGHT filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Outgoing Current to Right", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Outgoing Current to Right", & + to_str(t % scores(1,filter_index) % sum), & + trim(to_str(t % scores(1,filter_index) % sum_sq)) ! Back Surface t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j-1, k /) + 1, .true.) + mesh_indices_to_bin(m, (/ i, j-1, k /) + 1, .true.) t % matching_bins(i_filter_surf) = IN_FRONT filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Outgoing Current to Back", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Outgoing Current to Back", & + to_str(t % scores(1,filter_index) % sum), & + trim(to_str(t % scores(1,filter_index) % sum_sq)) t % matching_bins(i_filter_surf) = OUT_FRONT filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Incoming Current from Back", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Incoming Current from Back", & + to_str(t % scores(1,filter_index) % sum), & + trim(to_str(t % scores(1,filter_index) % sum_sq)) ! Front Surface t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) + mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) t % matching_bins(i_filter_surf) = IN_FRONT filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Incoming Current from Front", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Incoming Current from Front", & + to_str(t % scores(1,filter_index) % sum), & + trim(to_str(t % scores(1,filter_index) % sum_sq)) t % matching_bins(i_filter_surf) = OUT_FRONT filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Outgoing Current to Front", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Outgoing Current to Front", & + to_str(t % scores(1,filter_index) % sum), & + trim(to_str(t % scores(1,filter_index) % sum_sq)) ! Bottom Surface t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k-1 /) + 1, .true.) + mesh_indices_to_bin(m, (/ i, j, k-1 /) + 1, .true.) t % matching_bins(i_filter_surf) = IN_TOP filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Outgoing Current to Bottom", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Outgoing Current to Bottom", & + to_str(t % scores(1,filter_index) % sum), & + trim(to_str(t % scores(1,filter_index) % sum_sq)) t % matching_bins(i_filter_surf) = OUT_TOP filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Incoming Current from Bottom", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Incoming Current from Bottom", & + to_str(t % scores(1,filter_index) % sum), & + trim(to_str(t % scores(1,filter_index) % sum_sq)) ! Top Surface t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) + mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) t % matching_bins(i_filter_surf) = IN_TOP filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Incoming Current from Top", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Incoming Current from Top", & + to_str(t % scores(1,filter_index) % sum), & + trim(to_str(t % scores(1,filter_index) % sum_sq)) t % matching_bins(i_filter_surf) = OUT_TOP filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Outgoing Current to Top", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Outgoing Current to Top", & + to_str(t % scores(1,filter_index) % sum), & + trim(to_str(t % scores(1,filter_index) % sum_sq)) end do end do @@ -1793,10 +1831,10 @@ contains call bin_to_mesh_indices(m, bin, ijk) if (m % n_dimension == 2) then label = "Index (" // trim(to_str(ijk(1))) // ", " // & - trim(to_str(ijk(2))) // ")" + trim(to_str(ijk(2))) // ")" elseif (m % n_dimension == 3) then label = "Index (" // trim(to_str(ijk(1))) // ", " // & - trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")" + trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")" end if case (FILTER_ENERGYIN, FILTER_ENERGYOUT) E0 = t % filters(i_filter) % real_bins(bin) diff --git a/src/tally.F90 b/src/tally.F90 index daf3973e21..d7f05a5fb2 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -6,7 +6,7 @@ module tally use global use math, only: t_percentile, calc_pn use mesh, only: get_mesh_bin, bin_to_mesh_indices, get_mesh_indices, & - mesh_indices_to_bin, mesh_intersects + mesh_indices_to_bin, mesh_intersects use mesh_header, only: StructuredMesh use output, only: header use search, only: binary_search @@ -34,6 +34,9 @@ contains integer :: j ! loop index for scoring bins integer :: k ! loop index for nuclide bins + integer :: n ! loop index for scattering order + integer :: l ! scoring bin loop index, allowing for changing + ! position during the loop integer :: filter_index ! single index for single bin integer :: score_bin ! scoring bin, e.g. SCORE_FLUX integer :: i_nuclide ! index in nuclides array @@ -119,7 +122,9 @@ contains end if ! Determine score for each bin - SCORE_LOOP: do j = 1, t % n_score_bins + j = 0 + SCORE_LOOP: do l = 1, t % n_nonPN_score_bins + j = j + 1 ! determine what type of score bin score_bin = t % score_bins(j) @@ -160,52 +165,36 @@ contains ! reaction with neutrons in the exit channel score = wgt - - case (SCORE_SCATTER_1) + + case (SCORE_SCATTER_N) ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP - ! 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 + ! Find the scattering order for a singly requested moment, and + ! store its moment contribution. - score = last_wgt * mu + if (t % scatt_order(j) == 1) then + score = last_wgt * mu ! avoid function call overhead + else + score = last_wgt * calc_pn(t % scatt_order(j), mu) + endif - case (SCORE_SCATTER_2) + case (SCORE_SCATTER_PN) ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP + score_index = score_index - 1 + ! Find the scattering order for a collection of requested moments + ! and store the moment contribution of each + do n = 0, t % scatt_order(j) + ! determine scoring bin index + score_index = score_index + 1 + ! get the score and tally it + score = last_wgt * calc_pn(n, mu) + if (n < t % scatt_order(j)) & + call add_to_score(t % scores(score_index, filter_index), score) + end do + j = j + t % scatt_order(j) - ! The second scattering moment can be determined in a similar - ! manner to the first scattering moment - - score = last_wgt * calc_pn(2, mu) - - case (SCORE_SCATTER_3) - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP - - ! The third scattering moment can be determined in a similar - ! manner to the first scattering moment - - score = last_wgt * calc_pn(3, mu) - - case (SCORE_SCATTER_4) - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP - - ! The fourth scattering moment can be determined in a similar - ! manner to the first scattering moment - - score = last_wgt * calc_pn(4, mu) - - case (SCORE_SCATTER_5) - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP - - ! The fifth scattering moment can be determined in a similar - ! manner to the first scattering moment - - score = last_wgt * calc_pn(5, mu) case (SCORE_TRANSPORT) ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP @@ -233,7 +222,7 @@ contains ! mu*Sigma_s)). score = last_wgt / (3.0_8 * score * (material_xs % total - & - mu * score)) + mu * score)) case (SCORE_N_1N) ! Skip any event where the particle didn't scatter @@ -326,6 +315,7 @@ contains case default message = "Invalid score type on tally " // to_str(t % id) // "." + write(*,*) score_index, score_bin, l,j, t % n_score_bins call fatal_error() end select @@ -393,7 +383,7 @@ contains ! change outgoing energy bin t % matching_bins(i) = binary_search(t % filters(i) % real_bins, & - size(t % filters(i) % real_bins), E_out) + size(t % filters(i) % real_bins), E_out) ! determine scoring index i_filter = sum((t % matching_bins - 1) * t % stride) + 1 @@ -474,7 +464,7 @@ contains if (t % all_nuclides) then call score_all_nuclides(tracklength_tallies(curr_ptr % data), flux, & - filter_index) + filter_index) else NUCLIDE_BIN_LOOP: do k = 1, t % n_nuclide_bins @@ -511,20 +501,20 @@ contains select case(score_bin) case (SCORE_TOTAL) score = micro_xs(i_nuclide) % total * & - atom_density * flux + atom_density * flux case (SCORE_SCATTER) score = (micro_xs(i_nuclide) % total - & - micro_xs(i_nuclide) % absorption) * & - atom_density * flux + micro_xs(i_nuclide) % absorption) * & + atom_density * flux case (SCORE_ABSORPTION) score = micro_xs(i_nuclide) % absorption * & - atom_density * flux + atom_density * flux case (SCORE_FISSION) score = micro_xs(i_nuclide) % fission * & - atom_density * flux + atom_density * flux case (SCORE_NU_FISSION) score = micro_xs(i_nuclide) % nu_fission * & - atom_density * flux + atom_density * flux case (SCORE_EVENTS) score = ONE case default @@ -634,7 +624,7 @@ contains score = micro_xs(i_nuclide) % total * atom_density * flux case (SCORE_SCATTER) score = (micro_xs(i_nuclide) % total - & - micro_xs(i_nuclide) % absorption) * atom_density * flux + micro_xs(i_nuclide) % absorption) * atom_density * flux case (SCORE_ABSORPTION) score = micro_xs(i_nuclide) % absorption * atom_density * flux case (SCORE_FISSION) @@ -777,27 +767,27 @@ contains ! determine next universe bin ! TODO: Account for multiple universes when performing this filter t % matching_bins(i) = get_next_bin(FILTER_UNIVERSE, & - p % coord % universe, i_tally) + p % coord % universe, i_tally) case (FILTER_MATERIAL) t % matching_bins(i) = get_next_bin(FILTER_MATERIAL, & - p % material, i_tally) + p % material, i_tally) case (FILTER_CELL) ! determine next cell bin ! TODO: Account for cells in multiple levels when performing this filter t % matching_bins(i) = get_next_bin(FILTER_CELL, & - p % coord % cell, i_tally) + p % coord % cell, i_tally) case (FILTER_CELLBORN) ! determine next cellborn bin t % matching_bins(i) = get_next_bin(FILTER_CELLBORN, & - p % cell_born, i_tally) + p % cell_born, i_tally) case (FILTER_SURFACE) ! determine next surface bin t % matching_bins(i) = get_next_bin(FILTER_SURFACE, & - p % surface, i_tally) + p % surface, i_tally) case (FILTER_ENERGYIN) ! determine incoming energy bin @@ -805,12 +795,12 @@ contains ! check if energy of the particle is within energy bins if (p % E < t % filters(i) % real_bins(1) .or. & - p % E > t % filters(i) % real_bins(k + 1)) then + p % E > t % filters(i) % real_bins(k + 1)) then t % matching_bins(i) = NO_BIN_FOUND else ! search to find incoming energy bin t % matching_bins(i) = binary_search(t % filters(i) % real_bins, & - k + 1, p % E) + k + 1, p % E) end if end select @@ -931,25 +921,25 @@ contains select case(score_bin) case (SCORE_TOTAL) score = micro_xs(i_nuclide) % total * & - atom_density * flux + atom_density * flux case (SCORE_SCATTER) score = (micro_xs(i_nuclide) % total - & - micro_xs(i_nuclide) % absorption) * & - atom_density * flux + micro_xs(i_nuclide) % absorption) * & + atom_density * flux case (SCORE_ABSORPTION) score = micro_xs(i_nuclide) % absorption * & - atom_density * flux + atom_density * flux case (SCORE_FISSION) score = micro_xs(i_nuclide) % fission * & - atom_density * flux + atom_density * flux case (SCORE_NU_FISSION) score = micro_xs(i_nuclide) % nu_fission * & - atom_density * flux + atom_density * flux case (SCORE_EVENTS) score = ONE case default message = "Invalid score type on tally " // & - to_str(t % id) // "." + to_str(t % id) // "." call fatal_error() end select @@ -972,7 +962,7 @@ contains score = ONE case default message = "Invalid score type on tally " // & - to_str(t % id) // "." + to_str(t % id) // "." call fatal_error() end select end if @@ -1030,27 +1020,27 @@ contains ! determine next universe bin ! TODO: Account for multiple universes when performing this filter t % matching_bins(i) = get_next_bin(FILTER_UNIVERSE, & - p % coord % universe, i_tally) + p % coord % universe, i_tally) case (FILTER_MATERIAL) t % matching_bins(i) = get_next_bin(FILTER_MATERIAL, & - p % material, i_tally) + p % material, i_tally) case (FILTER_CELL) ! determine next cell bin ! TODO: Account for cells in multiple levels when performing this filter t % matching_bins(i) = get_next_bin(FILTER_CELL, & - p % coord % cell, i_tally) + p % coord % cell, i_tally) case (FILTER_CELLBORN) ! determine next cellborn bin t % matching_bins(i) = get_next_bin(FILTER_CELLBORN, & - p % cell_born, i_tally) + p % cell_born, i_tally) case (FILTER_SURFACE) ! determine next surface bin t % matching_bins(i) = get_next_bin(FILTER_SURFACE, & - p % surface, i_tally) + p % surface, i_tally) case (FILTER_ENERGYIN) ! determine incoming energy bin @@ -1065,12 +1055,12 @@ contains ! check if energy of the particle is within energy bins if (E < t % filters(i) % real_bins(1) .or. & - E > t % filters(i) % real_bins(n + 1)) then + E > t % filters(i) % real_bins(n + 1)) then t % matching_bins(i) = NO_BIN_FOUND else ! search to find incoming energy bin t % matching_bins(i) = binary_search(t % filters(i) % real_bins, & - n + 1, E) + n + 1, E) end if case (FILTER_ENERGYOUT) @@ -1079,12 +1069,12 @@ contains ! check if energy of the particle is within energy bins if (E < t % filters(i) % real_bins(1) .or. & - E > t % filters(i) % real_bins(n + 1)) then + E > t % filters(i) % real_bins(n + 1)) then t % matching_bins(i) = NO_BIN_FOUND else ! search to find incoming energy bin t % matching_bins(i) = binary_search(t % filters(i) % real_bins, & - n + 1, p % E) + n + 1, p % E) end if end select @@ -1175,14 +1165,14 @@ contains n = t % filters(j) % n_bins ! check if energy of the particle is within energy bins if (p % E < t % filters(j) % real_bins(1) .or. & - p % E > t % filters(j) % real_bins(n + 1)) then + p % E > t % filters(j) % real_bins(n + 1)) then curr_ptr => curr_ptr % next ! select next tally cycle end if ! search to find incoming energy bin t % matching_bins(j) = binary_search(t % filters(j) % real_bins, & - n + 1, p % E) + n + 1, p % E) end if ! ======================================================================= @@ -1200,7 +1190,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = OUT_TOP t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 call add_to_score(t % scores(1, filter_index), p % wgt) end if @@ -1211,7 +1201,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = IN_TOP t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 call add_to_score(t % scores(1, filter_index), p % wgt) end if @@ -1227,7 +1217,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = OUT_FRONT t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 call add_to_score(t % scores(1, filter_index), p % wgt) end if @@ -1238,7 +1228,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = IN_FRONT t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 call add_to_score(t % scores(1, filter_index), p % wgt) end if @@ -1254,7 +1244,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = OUT_RIGHT t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 call add_to_score(t % scores(1, filter_index), p % wgt) end if @@ -1265,7 +1255,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = IN_RIGHT t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 call add_to_score(t % scores(1, filter_index), p % wgt) end if @@ -1318,7 +1308,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = OUT_RIGHT t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) end if ijk0(1) = ijk0(1) + 1 xyz_cross(1) = xyz_cross(1) + m % width(1) @@ -1330,7 +1320,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = IN_RIGHT t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) end if end if elseif (distance == d(2)) then @@ -1340,7 +1330,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = OUT_FRONT t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) end if ijk0(2) = ijk0(2) + 1 xyz_cross(2) = xyz_cross(2) + m % width(2) @@ -1352,7 +1342,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = IN_FRONT t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) end if end if else if (distance == d(3)) then @@ -1362,7 +1352,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = OUT_TOP t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) end if ijk0(3) = ijk0(3) + 1 xyz_cross(3) = xyz_cross(3) + m % width(3) @@ -1374,7 +1364,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = IN_TOP t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) end if end if end if @@ -1385,7 +1375,7 @@ contains ! Check for errors if (filter_index <= 0 .or. filter_index > & - t % total_filter_bins) then + t % total_filter_bins) then message = "Score index outside range." call fatal_error() end if @@ -1442,7 +1432,7 @@ contains end if i_tally_check = tally_maps(filter_type) % items(filter_value) % & - elements(position(filter_type)) % index_tally + elements(position(filter_type)) % index_tally if (i_tally_check > i_tally) then ! Since the index being checked against is greater than the index we @@ -1454,7 +1444,7 @@ contains elseif (i_tally_check == i_tally) then ! Found a match bin = tally_maps(filter_type) % items(filter_value) % & - elements(position(filter_type)) % index_bin + elements(position(filter_type)) % index_bin return end if @@ -1539,14 +1529,14 @@ contains ! The MPI_IN_PLACE specifier allows the master to copy values into ! a receive buffer without having a temporary variable call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, & - MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) + MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) ! Transfer values to value on master t % scores(:,:) % value = tally_temp else ! Receive buffer not significant at other processors call MPI_REDUCE(tally_temp, dummy, n_bins, MPI_REAL8, & - MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) + MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) ! Reset value on other processors t % scores(:,:) % value = 0 @@ -1562,14 +1552,14 @@ contains if (master) then call MPI_REDUCE(MPI_IN_PLACE, global_temp, N_GLOBAL_TALLIES, & - MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) + MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) ! Transfer values back to global_tallies on master global_tallies(:) % value = global_temp else ! Receive buffer not significant at other processors call MPI_REDUCE(global_temp, dummy, N_GLOBAL_TALLIES, & - MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) + MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) ! Reset value on other processors global_tallies(:) % value = ZERO @@ -1580,11 +1570,11 @@ contains ! last realization if (master) then call MPI_REDUCE(MPI_IN_PLACE, total_weight, 1, MPI_REAL8, MPI_SUM, & - 0, MPI_COMM_WORLD, mpi_err) + 0, MPI_COMM_WORLD, mpi_err) else ! Receive buffer not significant at other processors call MPI_REDUCE(total_weight, dummy, 1, MPI_REAL8, MPI_SUM, & - 0, MPI_COMM_WORLD, mpi_err) + 0, MPI_COMM_WORLD, mpi_err) end if if (associated(curr_ptr)) nullify(curr_ptr) @@ -1690,7 +1680,7 @@ contains score % sum = score % sum/n score % sum_sq = sqrt((score % sum_sq/n - score % sum * & - score % sum) / (n - 1)) + score % sum) / (n - 1)) end subroutine statistics_score @@ -1911,7 +1901,7 @@ contains end if do i = n_cmfd_current_tallies + n_user_current_tallies, & - n_user_current_tallies + 1, -1 + n_user_current_tallies + 1, -1 ! allocate node allocate(curr_ptr) diff --git a/src/tally_header.F90 b/src/tally_header.F90 index f47f0ab1d9..47d992f4dd 100644 --- a/src/tally_header.F90 +++ b/src/tally_header.F90 @@ -97,8 +97,13 @@ module tally_header logical :: all_nuclides = .false. ! Values to score, e.g. flux, absorption, etc. + ! scat_order is the scattering order for each score. + ! It is to be 0 if the scattering order is 0, or if the score is not a + ! scattering response. integer :: n_score_bins = 0 integer, allocatable :: score_bins(:) + integer, allocatable :: scatt_order(:) + integer :: n_nonPN_score_bins = 0 ! Scores for each bin -- the first dimenion of the array is for scores ! (e.g. flux, total reaction rate, fission reaction rate, etc.) and the