diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 8cd11d7359..a1f19ecd8b 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -475,6 +475,8 @@ class StatePoint(object): deriv.material = self._f[base + '/material'].value deriv.nuclide = \ self._f[base + '/nuclide'].value.decode() + elif deriv.variable == 'temperature': + deriv.material = self._f[base + '/material'].value else: raise RuntimeError('Unrecognized tally differential ' 'variable') diff --git a/openmc/tally_derivative.py b/openmc/tally_derivative.py index 9fbad187dc..0b16ad46fa 100644 --- a/openmc/tally_derivative.py +++ b/openmc/tally_derivative.py @@ -33,7 +33,7 @@ class TallyDerivative(object): id : Integral Unique identifier for the tally derivative variable : str - Either 'density' or 'nuclide_density' + Accepted values are 'density', 'nuclide_density', and 'temperature' material : Integral The perutrubed material nuclide : str @@ -97,6 +97,9 @@ class TallyDerivative(object): self.material) string += '{0: <16}{1}{2}\n'.format('\tNuclide', '=\t', self.nuclide) + elif self.variable == 'temperature': + string += '{0: <16}{1}{2}\n'.format('\tMaterial', '=\t', + self.material) return string @@ -132,9 +135,9 @@ class TallyDerivative(object): def variable(self, var): if var is not None: cv.check_type('derivative variable', var, basestring) - if var not in ('density', 'nuclide_density'): - raise ValueError("A tally differential variable must be either " - "'density' or 'nuclide_density'") + if var not in ('density', 'nuclide_density', 'temperature'): + raise ValueError("A tally differential variable must be " + "'density', 'nuclide_density', or 'temperature'") self._variable = var @material.setter diff --git a/src/constants.F90 b/src/constants.F90 index 37f85b8303..bdd8c041d0 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -357,7 +357,8 @@ module constants ! Differential tally independent variables integer, parameter :: & DIFF_DENSITY = 1, & - DIFF_NUCLIDE_DENSITY = 2 + DIFF_NUCLIDE_DENSITY = 2, & + DIFF_TEMPERATURE = 3 ! ============================================================================ ! RANDOM NUMBER STREAM CONSTANTS diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 033e80b174..c739695001 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -729,6 +729,40 @@ contains end if end subroutine +!=============================================================================== +! MULTIPOLE_DERIV_EVAL evaluates the windowed multipole equations for the +! derivative of cross sections in the resolved resonance regions with respect to +! temperature. +!=============================================================================== + + subroutine multipole_deriv_eval(multipole, Emev, sqrtkT, sigT, sigA, sigF) + type(MultipoleArray), intent(in) :: multipole ! The windowed multipole + ! object to process. + real(8), intent(in) :: Emev ! The energy at which to + ! evaluate the cross section + ! derivative in MeV + real(8), intent(in) :: sqrtkT ! The temperature in the form + ! sqrt(kT (in MeV)), at which + ! to evaluate the XS. + real(8), intent(out) :: sigT ! Total cross section + real(8), intent(out) :: sigA ! Absorption cross section + real(8), intent(out) :: sigF ! Fission cross section + real(8), parameter :: T_DIFF = 10.0_8 + real(8) :: kT, sqrtkT_hi, sqrtkT_lo + real(8) :: sigT_hi, sigT_lo, sigA_hi, sigA_lo, sigF_hi, sigF_lo + + kT = sqrtkT**2 + sqrtkT_hi = sqrt(kT + K_BOLTZMANN*T_DIFF/2.0) + sqrtkT_lo = sqrt(kT - K_BOLTZMANN*T_DIFF/2.0) + + call multipole_eval(multipole, Emev, sqrtkT_hi, sigT_hi, sigA_hi, sigF_hi) + call multipole_eval(multipole, Emev, sqrtkT_lo, sigT_lo, sigA_lo, sigF_lo) + + sigT = (sigT_hi - sigT_lo) / T_DIFF + sigA = (sigA_hi - sigA_lo) / T_DIFF + sigF = (sigF_hi - sigF_lo) / T_DIFF + end subroutine multipole_deriv_eval + !=============================================================================== ! FILL_FACTORS calculates the value of phi, the hardsphere phase shift factor, ! and sigT_factor, a factor inside of the sigT equation not present in the diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 2d29734b72..6939923fb5 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1901,7 +1901,6 @@ contains type(Node), pointer :: doc => null() type(Node), pointer :: node_mat => null() type(Node), pointer :: node_dens => null() - type(Node), pointer :: node_temp => null() type(Node), pointer :: node_nuc => null() type(Node), pointer :: node_ele => null() type(Node), pointer :: node_sab => null() @@ -2642,6 +2641,10 @@ contains deallocate(pair_list) deriv % diff_nuclide = nuclide_dict % get_key(word) + + case("temperature") + deriv % variable = DIFF_TEMPERATURE + call get_node_value(node_deriv, "material", deriv % diff_material) end select end associate end do diff --git a/src/state_point.F90 b/src/state_point.F90 index b7a913dd01..ee78730342 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -212,6 +212,10 @@ contains i_list = nuclides(deriv % diff_nuclide) % listing call write_dataset(deriv_group, "nuclide", & xs_listings(i_list) % alias) + case (DIFF_TEMPERATURE) + call write_dataset(deriv_group, "independent variable", & + "temperature") + call write_dataset(deriv_group, "material", deriv % diff_material) case default call fatal_error("Independent variable for derivative " & // trim(to_str(deriv % id)) // " not defined in & diff --git a/src/tally.F90 b/src/tally.F90 index 4dea2e249f..a74e8c52e0 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -2,6 +2,7 @@ module tally use ace_header, only: Reaction use constants + use cross_section, only: multipole_deriv_eval use error, only: fatal_error use geometry_header use global @@ -2409,6 +2410,21 @@ contains call fatal_error("Differential tallies are only implemented for & &analog and collision estimators.") end select + + case (DIFF_TEMPERATURE) + select case (t % estimator) + + case (ESTIMATOR_ANALOG) + + case (ESTIMATOR_COLLISION) + + select case (score_bin) + + case (SCORE_FLUX) + score = score * deriv % flux_deriv + + end select + end select end select end associate end subroutine apply_derivative_to_score @@ -2422,7 +2438,8 @@ contains type(particle), intent(in) :: p real(8), intent(in) :: distance ! Neutron flight distance - integer :: i + integer :: i, l + real(8) :: dsigT, dsigA, dsigF ! A void material cannot be perturbed so it will not affect flux derivatives if (p % material == MATERIAL_VOID) return @@ -2452,6 +2469,24 @@ contains - distance * micro_xs(deriv % diff_nuclide) % total end if end associate + + case (DIFF_TEMPERATURE) + associate (mat => materials(p % material)) + if (mat % id == deriv % diff_material) then + do l=1, mat % n_nuclides + associate (nuc => nuclides(mat % nuclide(l))) + if (nuc % mp_present .and. & + p % E >= nuc % multipole % start_E/1.0e6_8 .and. & + p % E <= nuc % multipole % end_E/1.0e6_8) then + call multipole_deriv_eval(nuc % multipole, p % E, & + p % sqrtkT, dsigT, dsigA, dsigF) + deriv % flux_deriv = deriv % flux_deriv & + - distance * dsigT * mat % atom_density(l) + end if + end associate + end do + end if + end associate end select end associate end do