diff --git a/src/constants.f90 b/src/constants.f90 index 225116c6d7..cc67b297a7 100644 --- a/src/constants.f90 +++ b/src/constants.f90 @@ -153,7 +153,7 @@ module constants N_DA = 117 ! Tally macro reactions - integer, parameter :: N_MACRO_TYPES = 14 + integer, parameter :: N_MACRO_TYPES = 15 integer, parameter :: & MACRO_FLUX = -1, & ! flux MACRO_TOTAL = -2, & ! total reaction rate @@ -168,7 +168,8 @@ module constants 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_NU_FISSION = -14, & ! neutron production rate + MACRO_CURRENT = -15 ! partial current ! Tally map bin finding integer, parameter :: NO_BIN_FOUND = -1 diff --git a/src/input_xml.f90 b/src/input_xml.f90 index 2d3eb7e789..b97d7c5a1e 100644 --- a/src/input_xml.f90 +++ b/src/input_xml.f90 @@ -790,6 +790,46 @@ contains end if case ('nu-fission') t % macro_bins(j) % scalar = MACRO_NU_FISSION + case ('current') + t % macro_bins(j) % scalar = MACRO_CURRENT + t % surface_current = .true. + + ! Check to make sure that current is the only desired response + ! for this tally + if (n_words > 1) then + msg = "Cannot tally other macro reactions in the same " & + // "tally as surface currents. Separate other macro " & + // "reactions into a distinct tally." + call fatal_error(msg) + 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 + msg = "Surface currents must be used with a mesh filter only." + call fatal_error(msg) + 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) + + ! Get pointer to mesh + uid = t % mesh + index = dict_get_key(mesh_dict, uid) + m => meshes(index) + + ! 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 + elseif (size(m % dimension) == 3) then + t % n_bins(T_MESH) = t % n_bins(T_MESH) + & + product(m % dimension + 1) * 6 + end if + case default msg = "Unknown macro reaction: " // trim(words(j)) call fatal_error(msg) diff --git a/src/tally.f90 b/src/tally.f90 index 363a1a878d..db767a533c 100644 --- a/src/tally.f90 +++ b/src/tally.f90 @@ -274,6 +274,7 @@ contains integer :: score_index integer :: score_index0 integer :: macro_bin + integer :: mesh_bin real(8) :: score real(8) :: last_wgt real(8) :: wgt @@ -343,8 +344,17 @@ contains ! determine mesh bin if (t % n_bins(T_MESH) > 0) then m => meshes(t % mesh) - call get_mesh_bin(m, p % xyz, bins(T_MESH), in_mesh) + + ! Determine if we're in the mesh first + call get_mesh_bin(m, p % xyz, mesh_bin, in_mesh) if (.not. in_mesh) cycle + + if (t % surface_current) then + msg = "Surface current mesh tally not yet implemented." + call fatal_error(msg) + else + bins(T_MESH) = mesh_bin + end if else bins(T_MESH) = 1 end if diff --git a/src/tally_header.f90 b/src/tally_header.f90 index 9d745f0bbf..a525461795 100644 --- a/src/tally_header.f90 +++ b/src/tally_header.f90 @@ -61,6 +61,7 @@ module tally_header integer :: uid integer :: type real(8) :: volume + logical :: surface_current = .false. ! Tally bin specifications