From d31fd3470c43778ec29911ed9607b48be2ffe6e1 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 8 Nov 2015 22:03:30 -0500 Subject: [PATCH] Add diff tallies to statepoints, PyAPI --- openmc/statepoint.py | 23 +++++++++++ openmc/tallies.py | 92 +++++++++++++++++++++++++++++++++++++++++++- src/initialize.F90 | 10 ++++- src/output.F90 | 17 ++++++++ src/state_point.F90 | 32 ++++++++++++++- src/tally.F90 | 7 ++++ 6 files changed, 178 insertions(+), 3 deletions(-) diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 2edf9badd4..047581d1fc 100644 --- a/openmc/statepoint.py +++ b/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 diff --git a/openmc/tallies.py b/openmc/tallies.py index 0661ab67b1..4ada71519d 100644 --- a/openmc/tallies.py +++ b/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() diff --git a/src/initialize.F90 b/src/initialize.F90 index 6130be0e25..82252952d8 100644 --- a/src/initialize.F90 +++ b/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 diff --git a/src/output.F90 b/src/output.F90 index 0f667c121b..0043404805 100644 --- a/src/output.F90 +++ b/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) diff --git a/src/state_point.F90 b/src/state_point.F90 index 046e7e6669..84554cd947 100644 --- a/src/state_point.F90 +++ b/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 diff --git a/src/tally.F90 b/src/tally.F90 index 7fce90cce5..d5f36bbe74 100644 --- a/src/tally.F90 +++ b/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