diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 0192f30df..1ad4b5f9c 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -798,13 +798,19 @@ The following responses can be tallied. Total reaction rate :scatter: - Total scattering rate + Total scattering rate. Can also be identified with the ``scatter-0`` + response type. :nu-scatter: Total production of neutrons due to scattering. This accounts for multiplicity from (n,2n), (n,3n), and (n,4n) reactions and should be slightly higher than the scattering rate. + :scatter-0: + Zeroth scattering moment. Can also be identified with the ``scatter`` + response type, but is provided for consistency with the higher order + scattering moments. + :scatter-1: First scattering moment @@ -813,6 +819,12 @@ The following responses can be tallied. :scatter-3: Third scattering moment + + :scatter-4: + Fourth scattering moment + + :scatter-5: + Fifth scattering moment :absorption: Total absorption rate. This accounts for all reactions which do not produce diff --git a/src/constants.F90 b/src/constants.F90 index eb4e5e700..e1d8e853a 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -251,7 +251,7 @@ module constants EVENT_FISSION = 3 ! Tally score type - integer, parameter :: N_SCORE_TYPES = 18 + integer, parameter :: N_SCORE_TYPES = 20 integer, parameter :: & SCORE_FLUX = -1, & ! flux SCORE_TOTAL = -2, & ! total reaction rate @@ -260,17 +260,19 @@ module constants SCORE_SCATTER_1 = -5, & ! first scattering moment SCORE_SCATTER_2 = -6, & ! second scattering moment SCORE_SCATTER_3 = -7, & ! third scattering moment - SCORE_TRANSPORT = -8, & ! transport reaction rate - SCORE_DIFFUSION = -9, & ! diffusion coefficient - SCORE_N_1N = -10, & ! (n,1n) rate - SCORE_N_2N = -11, & ! (n,2n) rate - SCORE_N_3N = -12, & ! (n,3n) rate - SCORE_N_4N = -13, & ! (n,4n) rate - SCORE_ABSORPTION = -14, & ! absorption rate - SCORE_FISSION = -15, & ! fission rate - SCORE_NU_FISSION = -16, & ! neutron production rate - SCORE_CURRENT = -17, & ! partial current - SCORE_EVENTS = -18 ! number of events + 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 ! Tally map bin finding integer, parameter :: NO_BIN_FOUND = -1 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index c93aa4a33..621856834 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1700,6 +1700,8 @@ contains ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG + 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 @@ -1713,6 +1715,16 @@ contains case ('scatter-3') t % score_bins(j) = SCORE_SCATTER_3 + ! Set tally estimator to analog + t % estimator = ESTIMATOR_ANALOG + case ('scatter-4') + t % score_bins(j) = SCORE_SCATTER_4 + + ! Set tally estimator to analog + t % estimator = ESTIMATOR_ANALOG + case ('scatter-5') + t % score_bins(j) = SCORE_SCATTER_5 + ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG case('transport') diff --git a/src/math.F90 b/src/math.F90 index 7e36fb9ed..47a5e01f9 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -112,4 +112,39 @@ contains end function t_percentile +!=============================================================================== +! CALC_PN calculates the n-th order Legendre polynomial at the value of x. +! Since this function is called repeatedly during the neutron transport process, +! neither n or x is checked to see if they are in the applicable range. +! This is left to the client developer to use where applicable. x is to be in +! the domain of [-1,1], and 0<=n<=5. If x is outside of the range, the return +! value will be outside the expected range; if n is outside the stated range, +! the return value will be 1.0. +!=============================================================================== + + pure function calc_pn(n,x) result(pnx) + integer, intent(in) :: n ! Legendre order requested + real(8), intent(in) :: x ! Independent variable the Legendre is to be + ! evaluated at; x must be in the domain [-1,1] + real(8) :: pnx ! The Legendre poly of order n evaluated at x + + select case(n) + case(1) + pnx = x + case(2) + pnx = 1.5_8 * x * x - 0.5_8 + case(3) + pnx = 2.5_8 * x * x * x - 1.5_8 * x + case(4) + pnx = 4.375_8 * x * x * x * x - 3.75_8 * x * x + 0.375_8 + case(5) + pnx = 7.875_8 * x * x * x * x * x - 8.75_8 * x * x * x + 1.875 * x + case default + pnx = ONE + end select + + + end function calc_pn + + end module math diff --git a/src/output.F90 b/src/output.F90 index b72c859b7..04cfcb8a9 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -789,6 +789,10 @@ contains 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_TRANSPORT) string = trim(string) // ' transport' case (SCORE_DIFFUSION) @@ -1415,6 +1419,8 @@ contains 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" diff --git a/src/tally.F90 b/src/tally.F90 index c9596af39..daf3973e2 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -4,7 +4,7 @@ module tally use datatypes_header, only: ListInt use error, only: fatal_error use global - use math, only: t_percentile + 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 use mesh_header, only: StructuredMesh @@ -178,18 +178,34 @@ contains ! 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) + 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 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 + ! The third scattering moment can be determined in a similar + ! manner to the first scattering moment - score = last_wgt * 0.5*(5.0*mu*mu*mu - 3.0*mu) + 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