Add diff tallies to statepoints, PyAPI

openmc/statepoint.py

@@ -349,6 +349,29 @@ class StatePoint(object):
                     base, tally_key)].value.decode()
                 tally.num_realizations = n_realizations
 
+                # Read derivative information.
+                derivatives_present = self._f[
+                     '{0}{1}/derivative present'.format(base, tally_key)].value
+                assert derivatives_present in (0, 1)
+                if derivatives_present == 1:
+                    tally.diff_variable = self._f[
+                         '{0}{1}/derivative/dependent variable'.format(
+                         base, tally_key)].value.decode()
+                    if tally.diff_variable == 'density':
+                        tally.diff_material = self._f[
+                             '{0}{1}/derivative/material'.format(
+                             base, tally_key)].value
+                    elif tally.diff_variable == 'nuclide_density':
+                        tally.diff_material = self._f[
+                             '{0}{1}/derivative/material'.format(
+                             base, tally_key)].value
+                        tally.diff_nuclide = self._f[
+                             '{0}{1}/derivative/nuclide'.format(
+                             base, tally_key)].value.decode()
+                    else:
+                        raise RuntimeError('Unrecognized tally differential'
+                             'variable')
+
                 # Read the number of Filters
                 n_filters = self._f['{0}{1}/n_filters'.format(base, tally_key)].value
 

openmc/tallies.py

@@ -93,6 +93,9 @@ class Tally(object):
         self._scores = []
         self._estimator = None
         self._triggers = []
+        self._diff_variable = None
+        self._diff_material = None
+        self._diff_nuclide = None
 
         self._num_score_bins = 0
         self._num_realizations = 0
@@ -117,6 +120,9 @@ class Tally(object):
             clone.id = self.id
             clone.name = self.name
             clone.estimator = self.estimator
+            clone._diff_variable = self.diff_variable
+            clone._diff_material = self.diff_material
+            clone._diff_nuclide = self.diff_nuclide
             clone.num_score_bins = self.num_score_bins
             clone.num_realizations = self.num_realizations
             clone._sum = copy.deepcopy(self._sum, memo)
@@ -173,6 +179,14 @@ class Tally(object):
             if nuclide not in other.nuclides:
                 return False
 
+        # Check derivatives
+        if self.diff_variable != other.diff_variable:
+            return False
+        if self.diff_material != other.diff_material:
+            return False
+        if self.diff_nuclide != other.diff_nuclide:
+            return False
+
         # Check all scores
         if len(self.scores) != len(other.scores):
             return False
@@ -203,6 +217,9 @@ class Tally(object):
 
         hashable.append(self.estimator)
         hashable.append(self.name)
+        hashable.append(self._diff_variable)
+        hashable.append(self._diff_material)
+        hashable.append(self._diff_nuclide)
 
         return hash(tuple(hashable))
 
@@ -211,6 +228,21 @@ class Tally(object):
         string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self.id)
         string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self.name)
 
+        if self.diff_variable is not None:
+            string += '{0: <16}{1}{2}\n'.format('\tDerivative', '=\t',
+                                                self.diff_variable)
+            if self.diff_variable == 'density':
+                string += '{0: <16}{1}{2}\n'.format('\tDiff_material', '=\t',
+                                                    self.diff_material)
+            elif self.diff_variable == 'nuclide_density':
+                string += '{0: <16}{1}{2}\n'.format('\tDiff_material', '=\t',
+                                                    self.diff_material)
+                string += '{0: <16}{1}{2}\n'.format('\tDiff_nuclide', '=\t',
+                                                    self.diff_nuclide)
+            else:
+                raise RuntimeError("Encountered unrecognized differential "
+                                   "variable in a tally.")
+
         string += '{0: <16}{1}\n'.format('\tFilters', '=\t')
 
         for filter in self.filters:
@@ -379,6 +411,18 @@ class Tally(object):
     def derived(self):
         return self._derived
 
+    @property
+    def diff_variable(self):
+        return self._diff_variable
+
+    @property
+    def diff_material(self):
+        return self._diff_material
+
+    @property
+    def diff_nuclide(self):
+        return self._diff_nuclide
+
     @estimator.setter
     def estimator(self, estimator):
         cv.check_value('estimator', estimator,
@@ -422,6 +466,27 @@ class Tally(object):
         else:
             self._name = ''
 
+    @diff_variable.setter
+    def diff_variable(self, var):
+        if var is not None:
+            cv.check_type('differential variable', var, basestring)
+            if var not in ('density', 'nuclide_density'):
+                raise ValueError("A tally differential variable must be either "
+                     "'density' or 'nuclide_density'")
+            self._diff_variable = var
+
+    @diff_material.setter
+    def diff_material(self, mat):
+        if mat is not None:
+            cv.check_type('differential material', mat, Integral)
+            self._diff_material = mat
+
+    @diff_nuclide.setter
+    def diff_nuclide(self, nuc):
+        if nuc is not None:
+            cv.check_type('differential nuclide', nuc, basestring)
+            self._diff_nuclide = nuc
+
     def add_filter(self, filter):
         """Add a filter to the tally
 
@@ -703,6 +768,21 @@ class Tally(object):
             subelement = ET.SubElement(element, "nuclides")
             subelement.text = nuclides.rstrip(' ')
 
+        # Optional derivative
+        if self.diff_variable is not None:
+            if self.diff_variable == 'density':
+                subelement = ET.SubElement(element, "derivative")
+                subelement.set("variable", self.diff_variable)
+                subelement.set("material", str(self.diff_material))
+            elif self.diff_variable == 'nuclide_density':
+                subelement = ET.SubElement(element, "derivative")
+                subelement.set("variable", self.diff_variable)
+                subelement.set("material", str(self.diff_material))
+                subelement.set("nuclide", self.diff_nuclide)
+            else:
+                raise RuntimeError("Encountered unrecognized differential "
+                                   "variable in a tally.")
+
         # Scores
         if len(self.scores) == 0:
             msg = 'Unable to get XML for Tally ID="{0}" since it does not ' \
@@ -1122,7 +1202,7 @@ class Tally(object):
         return data
 
     def get_pandas_dataframe(self, filters=True, nuclides=True,
-                             scores=True, summary=None):
+                             scores=True, derivative=True, summary=None):
         """Build a Pandas DataFrame for the Tally data.
 
         This method constructs a Pandas DataFrame object for the Tally data
@@ -1140,6 +1220,8 @@ class Tally(object):
             Include columns with nuclide bin information (default is True).
         scores : bool
             Include columns with score bin information (default is True).
+        derivative : bool
+            Include columns with differential tally info (default is True).
         summary : None or Summary
             An optional Summary object to be used to construct columns for
             distribcell tally filters (default is None). The geometric
@@ -1221,6 +1303,14 @@ class Tally(object):
             tile_factor = data_size / len(self.scores)
             df['score'] = np.tile(self.scores, tile_factor)
 
+        # Include columns for derivatives if user requested it
+        if derivative and (self.diff_variable is not None):
+            df['d_variable'] = self.diff_variable
+            if self.diff_material is not None:
+                df['d_material'] = self.diff_material
+            if self.diff_nuclide is not None:
+                df['d_nuclide'] = self.diff_nuclide
+
         # Append columns with mean, std. dev. for each tally bin
         df['mean'] = self.mean.ravel()
         df['std. dev.'] = self.std_dev.ravel()

src/initialize.F90

@@ -831,7 +831,6 @@ contains
       ! Change density in g/cm^3 to atom/b-cm. Since all values are now in atom
       ! percent, the sum needs to be re-evaluated as 1/sum(x*awr)
       if (.not. density_in_atom) then
-        mat % density_gpcc = -mat % density
         sum_percent = ZERO
         do j = 1, mat%n_nuclides
           index_list = xs_listing_dict%get_key(mat%names(j))
@@ -846,6 +845,15 @@ contains
 
       ! Calculate nuclide atom densities
       mat%atom_density = mat%density * mat%atom_density
+
+      ! Calculate density in g/cm^3.
+      mat%density_gpcc = ZERO
+      do j = 1, mat%n_nuclides
+        index_list = xs_listing_dict%get_key(mat%names(j))
+        awr = xs_listings(index_list)%awr
+        mat%density_gpcc = mat%density_gpcc &
+             + mat%atom_density(j) * awr * MASS_NEUTRON / N_AVOGADRO
+      end do
     end do
 
   end subroutine normalize_ao

src/output.F90

@@ -1023,6 +1023,23 @@ contains
              // trim(t % name), unit=unit_tally, level=3)
       endif
 
+      ! Write derivative information.
+      if (allocated(t % deriv)) then
+        select case (t % deriv % dep_var)
+        case (DIFF_DENSITY)
+          write(unit=unit_tally, fmt="(' Density derivative  Material ',A)") &
+               to_str(t % deriv % diff_material)
+        case (DIFF_NUCLIDE_DENSITY)
+          i_listing = nuclides(t % deriv % diff_nuclide) % listing
+          write(unit=unit_tally, fmt="(' Nuclide density derivative  Material '&
+               &,A,'  Nuclide ',A)") trim(to_str(t % deriv % diff_material)), &
+               trim(xs_listings(i_listing) % alias)
+        case default
+          call fatal_error("Differential tally dependent variable for tally " &
+               // trim(to_str(t % id)) // " not defined in output.F90.")
+        end select
+      end if
+
       ! Handle surface current tallies separately
       if (t % type == TALLY_SURFACE_CURRENT) then
         call write_surface_current(t, unit_tally)

src/state_point.F90

@@ -51,7 +51,7 @@ contains
     integer(HID_T) :: cmfd_group
     integer(HID_T) :: tallies_group, tally_group
     integer(HID_T) :: meshes_group, mesh_group
-    integer(HID_T) :: filter_group
+    integer(HID_T) :: filter_group, deriv_group
     character(20), allocatable :: str_array(:)
     character(MAX_FILE_LEN)    :: filename
     type(RegularMesh), pointer :: meshp
@@ -307,6 +307,34 @@ contains
           call write_dataset(tally_group, "nuclides", str_array)
           deallocate(str_array)
 
+          ! Write derivative information.
+          if (allocated(tally % deriv)) then
+            call write_dataset(tally_group, "derivative present", 1)
+            deriv_group = create_group(tally_group, "derivative")
+            select case (tally % deriv % dep_var)
+            case (DIFF_DENSITY)
+              call write_dataset(deriv_group, "dependent variable", "density")
+              call write_dataset(deriv_group, "material", &
+                   tally % deriv % diff_material)
+            case (DIFF_NUCLIDE_DENSITY)
+              call write_dataset(deriv_group, "dependent variable", &
+                   "nuclide_density")
+              call write_dataset(deriv_group, "material", &
+                   tally % deriv % diff_material)
+              i_list = nuclides(tally % deriv % diff_nuclide) % listing
+              call write_dataset(deriv_group, "nuclide", &
+                   xs_listings(i_list) % alias)
+            case default
+              call fatal_error("Differential tally dependent variable for &
+                   &tally " // trim(to_str(tally % id)) // " not defined in &
+                   &state_point.F90.")
+            end select
+            call close_group(deriv_group)
+          else
+            call write_dataset(tally_group, "derivative present", 0)
+          end if
+
+          ! Write scores.
           call write_dataset(tally_group, "n_score_bins", tally%n_score_bins)
           allocate(str_array(size(tally%score_bins)))
           do j = 1, size(tally%score_bins)
@@ -355,6 +383,8 @@ contains
               str_array(j) = "events"
             case (SCORE_INVERSE_VELOCITY)
               str_array(j) = "inverse-velocity"
+            case (SCORE_KEFF)
+              str_array(j) = "keff"
             case default
               str_array(j) = reaction_name(tally%score_bins(j))
             end select

src/tally.F90

@@ -745,6 +745,9 @@ contains
       if (allocated(t % deriv) .and. t % estimator == ESTIMATOR_COLLISION) then
         select case (score_bin)
 
+        case (SCORE_FLUX)
+          score = score * t % deriv % accumulator
+
         case (SCORE_KEFF)
           select case (t % deriv % dep_var)
 
@@ -765,6 +768,10 @@ contains
             end if
 
           end select
+
+        case default
+          call fatal_error('Tally derivative not defined for a score on tally '&
+               // trim(to_str(t % id)))
         end select
       end if