Add finite-difference temperature derivative

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')

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

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

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

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

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 &

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