diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index 7551f4293f..5f34d4700c 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -162,6 +162,18 @@ if run_mode == 'k-eigenvalue': Width of each mesh cell in each dimension. +**/tallies/derivatives/derivative /independent variable** (*char[]*) + + Independent variable of tally derivative + +**/tallies/derivatives/derivative /material** (*int*) + + ID of the perturbed material + +**/tallies/derivatives/derivative /nuclide** (*char[]*) + + Alias of the perturbed nuclide + **/tallies/n_tallies** (*int*) Number of user-defined tallies. @@ -204,6 +216,10 @@ if run_mode == 'k-eigenvalue': Array of nuclides to tally. Note that if no nuclide is specified in the user input, a single 'total' nuclide appears here. +**/tallies/tally /derivative** (*int*) + + ID of the derivative applied to the tally. + **/tallies/tally /n_score_bins** (*int*) Number of scoring bins for a single nuclide. In general, this can be greater diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 92f53e6bf9..4a81615a10 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1916,6 +1916,13 @@ The ```` element accepts the following sub-elements: *Default*: "all" + :derivative: + The id of a ``derivative`` element. This derivative will be applied to all + scores in the tally. Differential tallies are currently only implemented + for collision and analog estimators. + + *Default*: None + ```` Element ------------------ @@ -1944,6 +1951,39 @@ attributes/sub-elements: One of ```` or ```` must be specified, but not both (even if they are consistent with one another). +```` Element +------------------------ + +OpenMC can take the first-order derivative of many tallies with respect to +material perturbations. It works by propagating a derivative through the +transport equation. Essentially, OpenMC keeps track of how each particle's +weight would change as materials are perturbed, and then accounts for that +weight change in the tallies. Note that this assumes material perturbations are +small enough not to change the distribution of fission sites. This element has +the following attributes/sub-elements: + + :id: + A unique integer that can be used to identify the derivative. + + :variable: + The independent variable of the derivative. Accepted options are "density", + "nuclide_density", and "temperature". A "density" derivative will give the + derivative with respect to the density of the material in [g / cm^3]. A + "nuclide_density" derivative will give the derivative with respect to the + density of a particular nuclide in units of [atom / b / cm]. A + "temperature" derivative is with respect to a material temperature in units + of [K]. The temperature derivative requires windowed multipole to be + turned on. Note also that the temperature derivative only accounts for + resolved resonance Doppler broadening. It does not account for thermal + expansion, S(a, b) scattering, resonance scattering, or unresolved Doppler + broadening. + + :material: + The perturbed material. (Necessary for all derivative types) + + :nuclide: + The perturbed nuclide. (Necessary only for "nuclide_density") + ```` Element ----------------------------- diff --git a/openmc/__init__.py b/openmc/__init__.py index 48f43be91d..b5530ed84b 100644 --- a/openmc/__init__.py +++ b/openmc/__init__.py @@ -16,6 +16,7 @@ from openmc.universe import * from openmc.mesh import * from openmc.filter import * from openmc.trigger import * +from openmc.tally_derivative import * from openmc.tallies import * from openmc.mgxs_library import * from openmc.cmfd import * diff --git a/openmc/statepoint.py b/openmc/statepoint.py index e05b18ebed..144146ba3a 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -97,6 +97,9 @@ class StatePoint(object): Dictionary whose keys are tally IDs and whose values are Tally objects tallies_present : bool Indicate whether user-defined tallies are present + tally_derivatives : dict + Dictionary whose keys are tally derivative IDs and whose values are + TallyDerivative objects version: tuple of Integral Version of OpenMC summary : None or openmc.Summary @@ -129,6 +132,7 @@ class StatePoint(object): self._summary = None self._global_tallies = None self._sparse = False + self._derivs_read = False # Automatically link in a summary file if one exists if autolink: @@ -391,6 +395,12 @@ class StatePoint(object): base, tally_key)].value.decode() tally.num_realizations = n_realizations + # Read derivative information. + if 'derivative' in self._f['{0}{1}'.format(base, tally_key)]: + deriv_id = self._f['{0}{1}/derivative'.format( + base, tally_key)].value + tally.derivative = self.tally_derivatives[deriv_id] + # Read the number of Filters n_filters = \ self._f['{0}{1}/n_filters'.format(base, tally_key)].value @@ -452,6 +462,43 @@ class StatePoint(object): def tallies_present(self): return self._f['tallies/tallies_present'].value + @property + def tally_derivatives(self): + if not self._derivs_read: + # Initialize dictionaries for the Meshes + # Keys - Derivative IDs + # Values - TallyDerivative objects + self._derivs = {} + + # Populate the dictionary if any derivatives are present. + if 'derivatives' in self._f['tallies']: + # Read the derivative ids. + base = 'tallies/derivatives' + deriv_ids = [int(k.split(' ')[1]) for k in self._f[base]] + + # Create each derivative object and add it to the dictionary. + for d_id in deriv_ids: + base = 'tallies/derivatives/derivative {:d}'.format(d_id) + deriv = openmc.TallyDerivative(derivative_id=d_id) + deriv.variable = \ + self._f[base + '/independent variable'].value.decode() + if deriv.variable == 'density': + deriv.material = self._f[base + '/material'].value + elif deriv.variable == 'nuclide_density': + 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') + self._derivs[d_id] = deriv + + self._derivs_read = True + + return self._derivs + @property def version(self): return (self._f['version_major'].value, diff --git a/openmc/tallies.py b/openmc/tallies.py index 28da08badf..d5d6df03fd 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -104,6 +104,8 @@ class Tally(object): sparse : bool Whether or not the tally uses SciPy's LIL sparse matrix format for compressed data storage + derivative : openmc.TallyDerivative + A material perturbation derivative to apply to all scores in the tally. """ @@ -116,6 +118,7 @@ class Tally(object): self._scores = cv.CheckedList(_SCORE_CLASSES, 'tally scores') self._estimator = None self._triggers = cv.CheckedList(openmc.Trigger, 'tally triggers') + self._derivative = None self._num_realizations = 0 self._with_summary = False @@ -151,6 +154,10 @@ class Tally(object): if nuclide not in other.nuclides: return False + # Check derivatives + if self.derivative != other.derivative: + return False + # Check all scores if len(self.scores) != len(other.scores): return False @@ -172,27 +179,31 @@ class Tally(object): def __repr__(self): string = 'Tally\n' - string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self.id) - string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self.name) + string += '{: <16}=\t{}\n'.format('\tID', self.id) + string += '{: <16}=\t{}\n'.format('\tName', self.name) - string += '{0: <16}{1}\n'.format('\tFilters', '=\t') + if self.derivative is not None: + string += '{: <16}=\t{}\n'.format('\tDerivative ID', + str(self.derivative.id)) + + string += '{: <16}=\n'.format('\tFilters') for self_filter in self.filters: - string += '{0: <16}\t\t{1}\t{2}\n'.format('', + string += '{: <16}\t\t{}\t{}\n'.format('', type(self_filter).__name__, self_filter.bins) - string += '{0: <16}{1}'.format('\tNuclides', '=\t') + string += '{: <16}=\t'.format('\tNuclides') for nuclide in self.nuclides: if isinstance(nuclide, openmc.Nuclide): - string += '{0} '.format(nuclide.name) + string += nuclide.name + ' ' else: - string += '{0} '.format(nuclide) + string += nuclide + ' ' string += '\n' - string += '{0: <16}{1}{2}\n'.format('\tScores', '=\t', self.scores) - string += '{0: <16}{1}{2}\n'.format('\tEstimator', '=\t', self.estimator) + string += '{: <16}=\t{}\n'.format('\tScores', self.scores) + string += '{: <16}=\t{}\n'.format('\tEstimator', self.estimator) return string @@ -375,6 +386,10 @@ class Tally(object): def derived(self): return self._derived + @property + def derivative(self): + return self._derivative + @property def sparse(self): return self._sparse @@ -429,6 +444,12 @@ class Tally(object): else: self._name = '' + @derivative.setter + def derivative(self, deriv): + if deriv is not None: + cv.check_type('tally derivative', deriv, openmc.TallyDerivative) + self._derivative = deriv + @filters.setter def filters(self, filters): cv.check_type('tally filters', filters, MutableSequence) @@ -1055,6 +1076,11 @@ class Tally(object): for trigger in self.triggers: trigger.get_trigger_xml(element) + # Optional derivatives + if self.derivative is not None: + subelement = ET.SubElement(element, "derivative") + subelement.text = str(self.derivative.id) + return element def contains_filter(self, filter_type): @@ -1489,7 +1515,8 @@ class Tally(object): return data def get_pandas_dataframe(self, filters=True, nuclides=True, scores=True, - distribcell_paths=True, float_format='{:.2e}'): + derivative=True, distribcell_paths=True, + float_format='{:.2e}'): """Build a Pandas DataFrame for the Tally data. This method constructs a Pandas DataFrame object for the Tally data @@ -1507,6 +1534,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). distribcell_paths : bool, optional Construct columns for distribcell tally filters (default is True). The geometric information in the Summary object is embedded into a @@ -1587,6 +1616,14 @@ class Tally(object): tile_factor = data_size / len(self.scores) df[column_name] = np.tile(scores, int(tile_factor)) + # Include columns for derivatives if user requested it + if derivative and (self.derivative is not None): + df['d_variable'] = self.derivative.variable + if self.derivative.material is not None: + df['d_material'] = self.derivative.material + if self.derivative.nuclide is not None: + df['d_nuclide'] = self.derivative.nuclide + # Append columns with mean, std. dev. for each tally bin df['mean'] = self.mean.ravel() df['std. dev.'] = self.std_dev.ravel() @@ -3554,6 +3591,18 @@ class Tallies(cv.CheckedList): self._tallies_file.append(xml_element) already_written.add(f.mesh) + def _create_derivative_subelements(self): + # Get a list of all derivatives referenced in a tally. + derivs = [] + for tally in self: + deriv = tally.derivative + if deriv is not None and deriv not in derivs: + derivs.append(deriv) + + # Add the derivatives to the XML tree. + for d in derivs: + self._tallies_file.append(d.to_xml_element()) + def export_to_xml(self): """Create a tallies.xml file that can be used for a simulation. @@ -3564,6 +3613,7 @@ class Tallies(cv.CheckedList): self._create_mesh_subelements() self._create_tally_subelements() + self._create_derivative_subelements() # Clean the indentation in the file to be user-readable clean_xml_indentation(self._tallies_file) diff --git a/openmc/tally_derivative.py b/openmc/tally_derivative.py new file mode 100644 index 0000000000..1c8a316cf5 --- /dev/null +++ b/openmc/tally_derivative.py @@ -0,0 +1,141 @@ +from __future__ import division + +import sys +from numbers import Integral +from xml.etree import ElementTree as ET + +from six import string_types + +import openmc.checkvalue as cv +from openmc.mixin import EqualityMixin + + +# "Static" variable for auto-generated TallyDerivative IDs +AUTO_TALLY_DERIV_ID = 10000 + +def reset_auto_tally_deriv_id(): + global AUTO_TALLY_ID + AUTO_TALLY_DERIV_ID = 10000 + + +class TallyDerivative(EqualityMixin): + """A material perturbation derivative to apply to a tally. + + Parameters + ---------- + derivative_id : Integral, optional + Unique identifier for the tally derivative. If none is specified, an + identifier will automatically be assigned + variable : str, optional + Accepted values are 'density', 'nuclide_density', and 'temperature' + material : Integral, optional + The perturubed material ID + nuclide : str, optional + The perturbed nuclide. Only needed for 'nuclide_density' derivatives. + Ex: 'Xe135' + + Attributes + ---------- + id : Integral + Unique identifier for the tally derivative + variable : str + Accepted values are 'density', 'nuclide_density', and 'temperature' + material : Integral + The perturubed material ID + nuclide : str + The perturbed nuclide. Only needed for 'nuclide_density' derivatives. + Ex: 'Xe135' + + """ + + def __init__(self, derivative_id=None, variable=None, material=None, + nuclide=None): + # Initialize Tally class attributes + self.id = derivative_id + self.variable = variable + self.material = material + self.nuclide = nuclide + + def __hash__(self): + return hash(repr(self)) + + def __repr__(self): + string = 'Tally Derivative\n' + string += '{: <16}=\t{}\n'.format('\tID', self.id) + string += '{: <16}=\t{}\n'.format('\tVariable', self.variable) + + if self.variable == 'density': + string += '{: <16}=\t{}\n'.format('\tMaterial', self.material) + elif self.variable == 'nuclide_density': + string += '{: <16}=\t{}\n'.format('\tMaterial', self.material) + string += '{: <16}=\t{}\n'.format('\tNuclide', self.nuclide) + elif self.variable == 'temperature': + string += '{: <16}=\t{}\n'.format('\tMaterial', self.material) + + return string + + @property + def id(self): + return self._id + + @property + def variable(self): + return self._variable + + @property + def material(self): + return self._material + + @property + def nuclide(self): + return self._nuclide + + @id.setter + def id(self, deriv_id): + if deriv_id is None: + global AUTO_TALLY_DERIV_ID + self._id = AUTO_TALLY_DERIV_ID + AUTO_TALLY_DERIV_ID += 1 + else: + cv.check_type('tally derivative ID', deriv_id, Integral) + cv.check_greater_than('tally derivative ID', deriv_id, 0, + equality=True) + self._id = deriv_id + + @variable.setter + def variable(self, var): + if var is not None: + cv.check_type('derivative variable', var, string_types) + cv.check_value('derivative variable', var, + ('density', 'nuclide_density', 'temperature')) + self._variable = var + + @material.setter + def material(self, mat): + if mat is not None: + cv.check_type('derivative material', mat, Integral) + self._material = mat + + @nuclide.setter + def nuclide(self, nuc): + if nuc is not None: + cv.check_type('derivative nuclide', nuc, string_types) + self._nuclide = nuc + + def to_xml_element(self): + """Return XML representation of the tally derivative + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing derivative data + + """ + + element = ET.Element("derivative") + element.set("id", str(self.id)) + element.set("variable", self.variable) + element.set("material", str(self.material)) + if self.variable == 'nuclide_density': + element.set("nuclide", self.nuclide) + return element diff --git a/src/constants.F90 b/src/constants.F90 index f2a9fb6f1a..044ea1c6bc 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -391,6 +391,12 @@ module constants K_TRACKLENGTH = 3, & LEAKAGE = 4 + ! Differential tally independent variables + integer, parameter :: & + DIFF_DENSITY = 1, & + DIFF_NUCLIDE_DENSITY = 2, & + DIFF_TEMPERATURE = 3 + ! ============================================================================ ! RANDOM NUMBER STREAM CONSTANTS diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 79ac84a54e..39552a7169 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -7,7 +7,7 @@ module cross_section use global use list_header, only: ListElemInt use material_header, only: Material - use math, only: faddeeva, broaden_wmp_polynomials + use math, only: faddeeva, w_derivative, broaden_wmp_polynomials use multipole_header, only: FORM_RM, FORM_MLBW, MP_EA, RM_RT, RM_RA, RM_RF, & MLBW_RT, MLBW_RX, MLBW_RA, MLBW_RF, FIT_T, FIT_A,& FIT_F, MultipoleArray @@ -706,6 +706,94 @@ contains end if end subroutine multipole_eval +!=============================================================================== +! 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, E, sqrtkT, sigT, sigA, sigF) + type(MultipoleArray), intent(in) :: multipole ! The windowed multipole + ! object to process. + real(8), intent(in) :: E ! The energy at which to + ! evaluate the cross section + real(8), intent(in) :: sqrtkT ! The temperature in the form + ! sqrt(kT), 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 + complex(8) :: w_val ! The faddeeva function evaluated at Z + complex(8) :: Z ! sqrt(atomic weight ratio / kT) * (sqrt(E) - pole) + complex(8) :: sigT_factor(multipole % num_l) + real(8) :: sqrtE ! sqrt(E), eV + real(8) :: invE ! 1/E, eV + real(8) :: dopp ! sqrt(atomic weight ratio / kT) + integer :: i_pole ! index of pole + integer :: i_window ! index of window + integer :: startw ! window start pointer (for poles) + integer :: endw ! window end pointer + real(8) :: T + + ! ========================================================================== + ! Bookkeeping + + ! Define some frequently used variables. + sqrtE = sqrt(E) + invE = ONE / E + dopp = multipole % sqrtAWR / sqrtkT + T = sqrtkT**2 / K_BOLTZMANN + + if (sqrtkT == ZERO) call fatal_error("Windowed multipole temperature & + &derivatives are not implemented for 0 Kelvin cross sections.") + + ! Locate us + i_window = floor((sqrtE - sqrt(multipole % start_E)) / multipole % spacing & + + ONE) + startw = multipole % w_start(i_window) + endw = multipole % w_end(i_window) + + ! Fill in factors. + if (startw <= endw) then + call compute_sigT_factor(multipole, sqrtE, sigT_factor) + end if + + ! Initialize the ouptut cross sections. + sigT = ZERO + sigA = ZERO + sigF = ZERO + + ! TODO Polynomials: Some of the curvefit polynomials Doppler broaden so + ! rigorously we should be computing the derivative of those. But in + ! practice, those derivatives are only large at very low energy and they + ! have no effect on reactor calculations. + + ! ========================================================================== + ! Add the contribution from the poles in this window. + + if (endw >= startw) then + do i_pole = startw, endw + Z = (sqrtE - multipole % data(MP_EA, i_pole)) * dopp + w_val = -invE * SQRT_PI * HALF * w_derivative(Z, 2) + if (multipole % formalism == FORM_MLBW) then + sigT = sigT + real((multipole % data(MLBW_RT, i_pole) * & + sigT_factor(multipole%l_value(i_pole)) + & + multipole % data(MLBW_RX, i_pole)) * w_val) + sigA = sigA + real(multipole % data(MLBW_RA, i_pole) * w_val) + sigF = sigF + real(multipole % data(MLBW_RF, i_pole) * w_val) + else if (multipole % formalism == FORM_RM) then + sigT = sigT + real(multipole % data(RM_RT, i_pole) * w_val * & + sigT_factor(multipole % l_value(i_pole))) + sigA = sigA + real(multipole % data(RM_RA, i_pole) * w_val) + sigF = sigF + real(multipole % data(RM_RF, i_pole) * w_val) + end if + end do + sigT = -HALF*multipole % sqrtAWR / sqrt(K_BOLTZMANN) * T**(-1.5) * sigT + sigA = -HALF*multipole % sqrtAWR / sqrt(K_BOLTZMANN) * T**(-1.5) * sigA + sigF = -HALF*multipole % sqrtAWR / sqrt(K_BOLTZMANN) * T**(-1.5) * sigF + end if + end subroutine multipole_deriv_eval + !=============================================================================== ! COMPUTE_SIGT_FACTOR calculates the sigT_factor, a factor inside of the sigT ! equation not present in the sigA and sigF equations. diff --git a/src/global.F90 b/src/global.F90 index eb4b33d7d6..02f6870189 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -20,7 +20,7 @@ module global use set_header, only: SetInt use surface_header, only: SurfaceContainer use source_header, only: SourceDistribution - use tally_header, only: TallyObject + use tally_header, only: TallyObject, TallyDerivative use trigger_header, only: KTrigger use timer_header, only: Timer use volume_header, only: VolumeCalculation @@ -185,6 +185,10 @@ module global integer :: n_tallies = 0 ! # of tallies integer :: n_user_tallies = 0 ! # of user tallies + ! Tally derivatives + type(TallyDerivative), allocatable :: tally_derivs(:) +!$omp threadprivate(tally_derivs) + ! Normalization for statistics integer :: n_realizations = 0 ! # of independent realizations real(8) :: total_weight ! total starting particle weight in realization diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 53cb11f985..1bbf68f48b 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2628,11 +2628,13 @@ contains type(Node), pointer :: node_mesh => null() type(Node), pointer :: node_tal => null() type(Node), pointer :: node_filt => null() - type(Node), pointer :: node_trigger=>null() + type(Node), pointer :: node_trigger => null() + type(Node), pointer :: node_deriv => null() type(NodeList), pointer :: node_mesh_list => null() type(NodeList), pointer :: node_tal_list => null() type(NodeList), pointer :: node_filt_list => null() type(NodeList), pointer :: node_trigger_list => null() + type(NodeList), pointer :: node_deriv_list => null() type(ElemKeyValueCI), pointer :: scores type(ElemKeyValueCI), pointer :: next @@ -2640,6 +2642,12 @@ contains filename = trim(path_input) // "tallies.xml" inquire(FILE=filename, EXIST=file_exists) if (.not. file_exists) then + ! We need to allocate tally_derivs to avoid segfaults. Also needs to be + ! done in parallel because tally derivs are threadprivate. +!$omp parallel + allocate(tally_derivs(0)) +!$omp end parallel + ! Since a tallies.xml file is optional, no error is issued here return end if @@ -2818,6 +2826,94 @@ contains ! We only need the mesh info for plotting if (run_mode == MODE_PLOTTING) return + ! ========================================================================== + ! READ DATA FOR DERIVATIVES + + ! Get pointer list to XML nodes and allocate global array. + ! The array is threadprivate so it must be allocated in parallel. + call get_node_list(doc, "derivative", node_deriv_list) +!$omp parallel + allocate(tally_derivs(get_list_size(node_deriv_list))) +!$omp end parallel + + ! Make sure this is not an MG run. + if (.not. run_CE .and. get_list_size(node_deriv_list) > 0) then + call fatal_error("Differential tallies not supported in multi-group mode") + end if + + ! Read derivative attributes. + do i = 1, get_list_size(node_deriv_list) + associate(deriv => tally_derivs(i)) + ! Get pointer to derivative node. + call get_list_item(node_deriv_list, i, node_deriv) + + ! Copy the derivative id. + if (check_for_node(node_deriv, "id")) then + call get_node_value(node_deriv, "id", deriv % id) + else + call fatal_error("Must specify an ID for elements in the& + & tally XML file") + end if + + ! Make sure the id is > 0. + if (deriv % id <= 0) then + call fatal_error(" IDs must be an integer greater than & + &zero") + end if + + ! Make sure this id has not already been used. + do j = 1, i-1 + if (tally_derivs(j) % id == deriv % id) then + call fatal_error("Two or more 's use the same unique & + &ID: " // trim(to_str(deriv % id))) + end if + end do + + ! Read the independent variable name. + temp_str = "" + call get_node_value(node_deriv, "variable", temp_str) + temp_str = to_lower(temp_str) + + select case(temp_str) + + case("density") + deriv % variable = DIFF_DENSITY + call get_node_value(node_deriv, "material", deriv % diff_material) + + case("nuclide_density") + deriv % variable = DIFF_NUCLIDE_DENSITY + call get_node_value(node_deriv, "material", deriv % diff_material) + + call get_node_value(node_deriv, "nuclide", word) + word = trim(to_lower(word)) + pair_list => nuclide_dict % keys() + do while (associated(pair_list)) + if (starts_with(pair_list % key, word)) then + word = pair_list % key(1:150) + exit + end if + + ! Advance to next + pair_list => pair_list % next + end do + + ! Check if no nuclide was found + if (.not. associated(pair_list)) then + call fatal_error("Could not find the nuclide " & + // trim(word) // " specified in derivative " & + // trim(to_str(deriv % id)) // " in any material.") + end if + 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 + ! ========================================================================== ! READ TALLY DATA @@ -2839,7 +2935,7 @@ contains t % estimator = ESTIMATOR_TRACKLENGTH - ! Copy material id + ! Copy tally id if (check_for_node(node_tal, "id")) then call get_node_value(node_tal, "id", t % id) else @@ -3429,6 +3525,7 @@ contains call fatal_error("Cannot tally flux with an outgoing energy & &filter.") end if + case ('flux-yn') ! Prohibit user from tallying flux for an individual nuclide if (.not. (t % n_nuclide_bins == 1 .and. & @@ -3766,6 +3863,48 @@ contains // trim(to_str(t % id)) // ".") end if + ! Check for a tally derivative. + if (check_for_node(node_tal, "derivative")) then + ! Temporarily store the derivative id. + call get_node_value(node_tal, "derivative", t % deriv) + + ! Find the derivative with the given id, and store it's index. + do j = 1, size(tally_derivs) + if (tally_derivs(j) % id == t % deriv) then + t % deriv = j + ! Only analog or collision estimators are supported for differential + ! tallies. + if (t % estimator == ESTIMATOR_TRACKLENGTH) then + t % estimator = ESTIMATOR_COLLISION + end if + ! We found the derivative we were looking for; exit the do loop. + exit + end if + if (j == size(tally_derivs)) then + call fatal_error("Could not find derivative " & + // trim(to_str(t % deriv)) // " specified on tally " & + // trim(to_str(t % id))) + end if + end do + + if (tally_derivs(t % deriv) % variable == DIFF_NUCLIDE_DENSITY & + .or. tally_derivs(t % deriv) % variable == DIFF_TEMPERATURE) then + if (any(t % nuclide_bins == -1)) then + if (t % find_filter(FILTER_ENERGYOUT) > 0) then + call fatal_error("Error on tally " // trim(to_str(t % id)) & + // ": Cannot use a 'nuclide_density' or 'temperature' & + &derivative on a tally with an outgoing energy filter and & + &'total' nuclide rate. Instead, tally each nuclide in the & + &material individually.") + ! Note that diff tallies with these characteristics would work + ! correctly if no tally events occur in the perturbed material + ! (e.g. pertrubing moderator but only tallying fuel), but this + ! case would be hard to check for by only reading inputs. + end if + end if + end if + end if + ! If settings.xml trigger is turned on, create tally triggers if (trigger_on) then @@ -4787,8 +4926,8 @@ contains sum_percent = sum(mat % atom_density) mat % atom_density = mat % atom_density / sum_percent - ! 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) + ! 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 sum_percent = ZERO do j = 1, mat % n_nuclides @@ -4807,6 +4946,18 @@ 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 + if (run_CE) then + awr = nuclides(mat % nuclide(j)) % awr + else + awr = ONE + end if + mat % density_gpcc = mat % density_gpcc & + + mat % atom_density(j) * awr * MASS_NEUTRON / N_AVOGADRO + end do end associate end do diff --git a/src/material_header.F90 b/src/material_header.F90 index 772e3a415d..f7c1b5db09 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -13,6 +13,7 @@ module material_header integer, allocatable :: nuclide(:) ! index in nuclides array real(8) :: density ! total atom density in atom/b-cm real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm + real(8) :: density_gpcc ! total density in g/cm^3 ! Energy grid information integer :: n_grid ! # of union material grid points diff --git a/src/math.F90 b/src/math.F90 index 11ace6976a..eac4ea94cc 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -753,6 +753,22 @@ contains end function faddeeva + recursive function w_derivative(z, order) result(wv) + complex(C_DOUBLE_COMPLEX), intent(in) :: z ! The point to evaluate Z at + integer, intent(in) :: order + complex(8) :: wv ! The resulting w(z) value + + select case(order) + case (0) + wv = faddeeva(z) + case (1) + wv = -TWO * z * faddeeva(z) + TWO * ONEI / SQRT_PI + case default + wv = -TWO * z * w_derivative(z, order-1) & + - TWO * (order-1) * w_derivative(z, order-2) + end select + end function w_derivative + !=============================================================================== ! BROADEN_WMP_POLYNOMIALS Doppler broadens the windowed multipole curvefit. The ! curvefit is a polynomial of the form diff --git a/src/output.F90 b/src/output.F90 index 6a8ba23b58..f5766187ee 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -777,6 +777,28 @@ contains // trim(t % name), unit=unit_tally, level=3) endif + ! Write derivative information. + if (t % deriv /= NONE) then + associate(deriv => tally_derivs(t % deriv)) + select case (deriv % variable) + case (DIFF_DENSITY) + write(unit=unit_tally, fmt="(' Density derivative Material ',A)") & + to_str(deriv % diff_material) + case (DIFF_NUCLIDE_DENSITY) + write(unit=unit_tally, fmt="(' Nuclide density derivative & + &Material ',A,' Nuclide ',A)") & + trim(to_str(deriv % diff_material)), & + trim(nuclides(deriv % diff_nuclide) % name) + case (DIFF_TEMPERATURE) + write(unit=unit_tally, fmt="(' Temperature derivative Material ',& + &A)") to_str(deriv % diff_material) + case default + call fatal_error("Differential tally dependent variable for tally "& + // trim(to_str(t % id)) // " not defined in output.F90.") + end select + end associate + 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/relaxng/tallies.rnc b/src/relaxng/tallies.rnc index cc327f8048..2143c8e472 100644 --- a/src/relaxng/tallies.rnc +++ b/src/relaxng/tallies.rnc @@ -15,6 +15,23 @@ element tallies { ) }* & + element derivative { + (element id { xsd:int } | attribute id { xsd:int }) & + (element material { xsd:int } | attribute material { xsd:int }) & + ( (element variable { ( "density") } + | attribute variable { ( "density" ) } ) | + ( + (element variable { ( "nuclide_density" ) } + | attribute variable { ( "nuclide_density" ) } ) + & + (element nuclide { xsd:string { maxLength = "12" } } + | attribute nuclide { xsd:string { maxLength = "12" } } ) | + ) + (element variable { ( "temperature") } + | attribute variable { ( "temperature" ) } ) + ) + }* & + element tally { (element id { xsd:int } | attribute id { xsd:int }) & (element name { xsd:string { maxLength="52" } } | @@ -41,7 +58,8 @@ element tallies { (element type { xsd:string } | attribute type { xsd:string }) & (element threshold { xsd:double} | attribute threshold { xsd:double }) & (element scores { list { xsd:string { maxLength = "20" }+ } } | attribute scores { list { xsd:string { maxLength = "20"}+ } } )? - }* + }* & + (element derivative { xsd:int } | attribute derivative { xsd:int } )? }* & element assume_separate { xsd:boolean }? diff --git a/src/relaxng/tallies.rng b/src/relaxng/tallies.rng index 755e9e90dc..fc60b8373a 100644 --- a/src/relaxng/tallies.rng +++ b/src/relaxng/tallies.rng @@ -89,6 +89,68 @@ + + + + + + + + + + + + + + + + + + + + + + + density + + + density + + + + + + nuclide_density + + + nuclide_density + + + + + + 12 + + + + + 12 + + + + + + + temperature + + + temperature + + + + + + @@ -254,6 +316,16 @@ + + + + + + + + + + diff --git a/src/simulation.F90 b/src/simulation.F90 index ef055a5493..e6b20ddf8f 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -83,7 +83,7 @@ contains ! ==================================================================== ! LOOP OVER PARTICLES -!$omp parallel do schedule(static) firstprivate(p) +!$omp parallel do schedule(static) firstprivate(p) copyin(tally_derivs) PARTICLE_LOOP: do i_work = 1, work current_work = i_work diff --git a/src/state_point.F90 b/src/state_point.F90 index f9616f2285..23543775cb 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -50,7 +50,8 @@ contains integer, allocatable :: key_array(:) integer(HID_T) :: file_id integer(HID_T) :: cmfd_group, tallies_group, tally_group, meshes_group, & - mesh_group, filter_group, runtime_group + mesh_group, filter_group, derivs_group, deriv_group, & + runtime_group character(MAX_WORD_LEN), allocatable :: str_array(:) character(MAX_FILE_LEN) :: filename type(RegularMesh), pointer :: meshp @@ -198,6 +199,38 @@ contains call close_group(meshes_group) + ! Write information for derivatives. + if (size(tally_derivs) > 0) then + derivs_group = create_group(tallies_group, "derivatives") + do i = 1, size(tally_derivs) + associate(deriv => tally_derivs(i)) + deriv_group = create_group(derivs_group, "derivative " & + // trim(to_str(deriv % id))) + select case (deriv % variable) + case (DIFF_DENSITY) + call write_dataset(deriv_group, "independent variable", "density") + call write_dataset(deriv_group, "material", deriv % diff_material) + case (DIFF_NUCLIDE_DENSITY) + call write_dataset(deriv_group, "independent variable", & + "nuclide_density") + call write_dataset(deriv_group, "material", deriv % diff_material) + call write_dataset(deriv_group, "nuclide", & + nuclides(deriv % diff_nuclide) % name) + 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 & + &state_point.F90.") + end select + call close_group(deriv_group) + end associate + end do + call close_group(derivs_group) + end if + ! Write number of tallies call write_dataset(tallies_group, "n_tallies", n_tallies) @@ -273,6 +306,13 @@ contains call write_dataset(tally_group, "nuclides", str_array) deallocate(str_array) + ! Write derivative information. + if (tally % deriv /= NONE) then + call write_dataset(tally_group, "derivative", & + tally_derivs(tally % deriv) % id) + 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) diff --git a/src/tally.F90 b/src/tally.F90 index d4968d9d68..79f22610e0 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -8,6 +8,7 @@ module tally use algorithm, only: binary_search use constants + use cross_section, only: multipole_deriv_eval use error, only: fatal_error use geometry_header use global @@ -1087,6 +1088,14 @@ contains end select + !######################################################################### + ! Add derivative information on score for differential tallies. + + if (t % deriv /= NONE) then + call apply_derivative_to_score(p, t, i_nuclide, atom_density, & + score_bin, score) + end if + !######################################################################### ! Expand score if necessary and add to tally results. call expand_and_score(p, t, score_index, filter_index, score_bin, & @@ -2420,6 +2429,13 @@ contains ! determine score based on bank site weight and keff score = keff * fission_bank(n_bank - p % n_bank + k) % wgt + ! Add derivative information for differential tallies. Note that the + ! i_nuclide and atom_density arguments do not matter since this is an + ! analog estimator. + if (t % deriv /= NONE) then + call apply_derivative_to_score(p, t, 0, ZERO, SCORE_NU_FISSION, score) + end if + if (.not. run_CE .and. eo_filt % matches_transport_groups) then ! determine outgoing energy from fission bank @@ -3107,6 +3123,778 @@ contains end subroutine score_surface_current +!=============================================================================== +! APPLY_DERIVATIVE_TO_SCORE multiply the given score by its relative derivative +!=============================================================================== + + subroutine apply_derivative_to_score(p, t, i_nuclide, atom_density, & + score_bin, score) + type(Particle), intent(in) :: p + type(TallyObject), intent(in) :: t + integer, intent(in) :: i_nuclide + real(8), intent(in) :: atom_density ! atom/b-cm + integer, intent(in) :: score_bin + real(8), intent(inout) :: score + + integer :: l + logical :: scoring_diff_nuclide + real(8) :: flux_deriv + real(8) :: dsigT, dsigA, dsigF, cum_dsig + + if (score == ZERO) return + + ! If our score was previously c then the new score is + ! c * (1/f * d_f/d_p + 1/c * d_c/d_p) + ! where (1/f * d_f/d_p) is the (logarithmic) flux derivative and p is the + ! perturbated variable. + + associate(deriv => tally_derivs(t % deriv)) + flux_deriv = deriv % flux_deriv + + select case (tally_derivs(t % deriv) % variable) + + !========================================================================= + ! Density derivative: + ! c = Sigma_MT + ! c = sigma_MT * N + ! c = sigma_MT * rho * const + ! d_c / d_rho = sigma_MT * const + ! (1 / c) * (d_c / d_rho) = 1 / rho + + case (DIFF_DENSITY) + select case (t % estimator) + + case (ESTIMATOR_ANALOG) + + select case (score_bin) + + case (SCORE_FLUX) + score = score * flux_deriv + + case (SCORE_TOTAL, SCORE_SCATTER, SCORE_ABSORPTION, SCORE_FISSION, & + SCORE_NU_FISSION) + if (materials(p % material) % id == deriv % diff_material) then + score = score * (flux_deriv + ONE & + / materials(p % material) % density_gpcc) + else + score = score * flux_deriv + end if + + case default + call fatal_error('Tally derivative not defined for a score on & + &tally ' // trim(to_str(t % id))) + end select + + case (ESTIMATOR_COLLISION) + + select case (score_bin) + + case (SCORE_FLUX) + score = score * flux_deriv + + case (SCORE_TOTAL, SCORE_SCATTER, SCORE_ABSORPTION, SCORE_FISSION, & + SCORE_NU_FISSION) + if (materials(p % material) % id == deriv % diff_material) then + score = score * (flux_deriv + ONE & + / materials(p % material) % density_gpcc) + else + score = score * flux_deriv + end if + + case default + call fatal_error('Tally derivative not defined for a score on & + &tally ' // trim(to_str(t % id))) + end select + + case default + call fatal_error("Differential tallies are only implemented for & + &analog and collision estimators.") + end select + + !========================================================================= + ! Nuclide density derivative: + ! If we are scoring a reaction rate for a single nuclide then + ! c = Sigma_MT_i + ! c = sigma_MT_i * N_i + ! d_c / d_N_i = sigma_MT_i + ! (1 / c) * (d_c / d_N_i) = 1 / N_i + ! If the score is for the total material (i_nuclide = -1) + ! c = Sum_i(Sigma_MT_i) + ! d_c / d_N_i = sigma_MT_i + ! (1 / c) * (d_c / d_N) = sigma_MT_i / Sigma_MT + ! where i is the perturbed nuclide. + + case (DIFF_NUCLIDE_DENSITY) + select case (t % estimator) + + case (ESTIMATOR_ANALOG) + + select case (score_bin) + + case (SCORE_FLUX) + score = score * flux_deriv + + case (SCORE_TOTAL, SCORE_SCATTER, SCORE_ABSORPTION, SCORE_FISSION, & + SCORE_NU_FISSION) + if (materials(p % material) % id == deriv % diff_material & + .and. p % event_nuclide == deriv % diff_nuclide) then + associate(mat => materials(p % material)) + ! Search for the index of the perturbed nuclide. + do l = 1, mat % n_nuclides + if (mat % nuclide(l) == deriv % diff_nuclide) exit + end do + + score = score * (flux_deriv & + + ONE / mat % atom_density(l)) + end associate + else + score = score * flux_deriv + end if + + case default + call fatal_error('Tally derivative not defined for a score on & + &tally ' // trim(to_str(t % id))) + end select + + case (ESTIMATOR_COLLISION) + scoring_diff_nuclide = & + (materials(p % material) % id == deriv % diff_material) & + .and. (i_nuclide == deriv % diff_nuclide) + + select case (score_bin) + + case (SCORE_FLUX) + score = score * flux_deriv + + case (SCORE_TOTAL) + if (i_nuclide == -1 .and. & + materials(p % material) % id == deriv % diff_material .and. & + material_xs % total /= ZERO) then + score = score * (flux_deriv & + + micro_xs(deriv % diff_nuclide) % total & + / material_xs % total) + else if (scoring_diff_nuclide .and. & + micro_xs(deriv % diff_nuclide) % total /= ZERO) then + score = score * (flux_deriv + ONE / atom_density) + else + score = score * flux_deriv + end if + + case (SCORE_SCATTER) + if (i_nuclide == -1 .and. & + materials(p % material) % id == deriv % diff_material .and. & + material_xs % total - material_xs % absorption /= ZERO) then + score = score * (flux_deriv & + + (micro_xs(deriv % diff_nuclide) % total & + - micro_xs(deriv % diff_nuclide) % absorption) & + / (material_xs % total - material_xs % absorption)) + else if (scoring_diff_nuclide .and. & + (micro_xs(deriv % diff_nuclide) % total & + - micro_xs(deriv % diff_nuclide) % absorption) /= ZERO) then + score = score * (flux_deriv + ONE / atom_density) + else + score = score * flux_deriv + end if + + case (SCORE_ABSORPTION) + if (i_nuclide == -1 .and. & + materials(p % material) % id == deriv % diff_material .and. & + material_xs % absorption /= ZERO) then + score = score * (flux_deriv & + + micro_xs(deriv % diff_nuclide) % absorption & + / material_xs % absorption ) + else if (scoring_diff_nuclide .and. & + micro_xs(deriv % diff_nuclide) % absorption /= ZERO) then + score = score * (flux_deriv + ONE / atom_density) + else + score = score * flux_deriv + end if + + case (SCORE_FISSION) + if (i_nuclide == -1 .and. & + materials(p % material) % id == deriv % diff_material .and. & + material_xs % fission /= ZERO) then + score = score * (flux_deriv & + + micro_xs(deriv % diff_nuclide) % fission & + / material_xs % fission) + else if (scoring_diff_nuclide .and. & + micro_xs(deriv % diff_nuclide) % fission /= ZERO) then + score = score * (flux_deriv + ONE / atom_density) + else + score = score * flux_deriv + end if + + case (SCORE_NU_FISSION) + if (i_nuclide == -1 .and. & + materials(p % material) % id == deriv % diff_material .and. & + material_xs % nu_fission /= ZERO) then + score = score * (flux_deriv & + + micro_xs(deriv % diff_nuclide) % nu_fission & + / material_xs % nu_fission) + else if (scoring_diff_nuclide .and. & + micro_xs(deriv % diff_nuclide) % nu_fission /= ZERO) then + score = score * (flux_deriv + ONE / atom_density) + else + score = score * flux_deriv + end if + + case default + call fatal_error('Tally derivative not defined for a score on & + &tally ' // trim(to_str(t % id))) + end select + + case default + call fatal_error("Differential tallies are only implemented for & + &analog and collision estimators.") + end select + + !========================================================================= + ! Temperature derivative: + ! If we are scoring a reaction rate for a single nuclide then + ! c = Sigma_MT_i + ! c = sigma_MT_i * N_i + ! d_c / d_T = (d_sigma_Mt_i / d_T) * N_i + ! (1 / c) * (d_c / d_T) = (d_sigma_MT_i / d_T) / sigma_MT_i + ! If the score is for the total material (i_nuclide = -1) + ! (1 / c) * (d_c / d_T) = Sum_i((d_sigma_MT_i / d_T) * N_i) / Sigma_MT_i + ! where i is the perturbed nuclide. The d_sigma_MT_i / d_T term is + ! computed by multipole_deriv_eval. It only works for the resolved + ! resonance range and requires multipole data. + + case (DIFF_TEMPERATURE) + select case (t % estimator) + + case (ESTIMATOR_ANALOG) + + select case (score_bin) + + case (SCORE_FLUX) + score = score * flux_deriv + + case (SCORE_TOTAL) + if (materials(p % material) % id == deriv % diff_material .and. & + micro_xs(p % event_nuclide) % total > ZERO) then + associate(mat => materials(p % material)) + ! Search for the index of the perturbed nuclide. + do l = 1, mat % n_nuclides + if (mat % nuclide(l) == p % event_nuclide) exit + end do + + dsigT = ZERO + associate (nuc => nuclides(p % event_nuclide)) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + end if + end associate + score = score * (flux_deriv & + + dsigT * mat % atom_density(l) / material_xs % total) + end associate + else + score = score * flux_deriv + end if + + case (SCORE_SCATTER) + if (materials(p % material) % id == deriv % diff_material .and. & + (micro_xs(p % event_nuclide) % total & + - micro_xs(p % event_nuclide) % absorption) > ZERO) then + associate(mat => materials(p % material)) + ! Search for the index of the perturbed nuclide. + do l = 1, mat % n_nuclides + if (mat % nuclide(l) == p % event_nuclide) exit + end do + + dsigT = ZERO + dsigA = ZERO + associate (nuc => nuclides(p % event_nuclide)) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + end if + end associate + score = score * (flux_deriv + (dsigT - dsigA) & + * mat % atom_density(l) / & + (material_xs % total - material_xs % absorption)) + end associate + else + score = score * flux_deriv + end if + + case (SCORE_ABSORPTION) + if (materials(p % material) % id == deriv % diff_material .and. & + micro_xs(p % event_nuclide) % absorption > ZERO) then + associate(mat => materials(p % material)) + ! Search for the index of the perturbed nuclide. + do l = 1, mat % n_nuclides + if (mat % nuclide(l) == p % event_nuclide) exit + end do + + dsigA = ZERO + associate (nuc => nuclides(p % event_nuclide)) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + end if + end associate + score = score * (flux_deriv & + + dsigA * mat % atom_density(l) / material_xs % absorption) + end associate + else + score = score * flux_deriv + end if + + case (SCORE_FISSION) + if (materials(p % material) % id == deriv % diff_material .and. & + micro_xs(p % event_nuclide) % fission > ZERO) then + associate(mat => materials(p % material)) + ! Search for the index of the perturbed nuclide. + do l = 1, mat % n_nuclides + if (mat % nuclide(l) == p % event_nuclide) exit + end do + + dsigF = ZERO + associate (nuc => nuclides(p % event_nuclide)) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + end if + end associate + score = score * (flux_deriv & + + dsigF * mat % atom_density(l) / material_xs % fission) + end associate + else + score = score * flux_deriv + end if + + case (SCORE_NU_FISSION) + if (materials(p % material) % id == deriv % diff_material .and. & + micro_xs(p % event_nuclide) % nu_fission > ZERO) then + associate(mat => materials(p % material)) + ! Search for the index of the perturbed nuclide. + do l = 1, mat % n_nuclides + if (mat % nuclide(l) == p % event_nuclide) exit + end do + + dsigF = ZERO + associate (nuc => nuclides(p % event_nuclide)) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + end if + end associate + score = score * (flux_deriv & + + dsigF * mat % atom_density(l) / material_xs % nu_fission& + * micro_xs(p % event_nuclide) % nu_fission & + / micro_xs(p % event_nuclide) % fission) + end associate + else + score = score * flux_deriv + end if + + case default + call fatal_error('Tally derivative not defined for a score on & + &tally ' // trim(to_str(t % id))) + end select + + case (ESTIMATOR_COLLISION) + + select case (score_bin) + + case (SCORE_FLUX) + score = score * flux_deriv + + case (SCORE_TOTAL) + if (i_nuclide == -1 .and. & + materials(p % material) % id == deriv % diff_material .and. & + material_xs % total > ZERO) then + cum_dsig = ZERO + associate(mat => materials(p % material)) + do l = 1, mat % n_nuclides + associate (nuc => nuclides(mat % nuclide(l))) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E .and. & + micro_xs(mat % nuclide(l)) % total > ZERO) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + cum_dsig = cum_dsig + dsigT * mat % atom_density(l) + end if + end associate + end do + end associate + score = score * (flux_deriv & + + cum_dsig / material_xs % total) + else if (materials(p % material) % id == deriv % diff_material & + .and. material_xs % total > ZERO) then + dsigT = ZERO + associate (nuc => nuclides(i_nuclide)) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + end if + end associate + score = score * (flux_deriv & + + dsigT / micro_xs(i_nuclide) % total) + else + score = score * flux_deriv + end if + + case (SCORE_SCATTER) + if (i_nuclide == -1 .and. & + materials(p % material) % id == deriv % diff_material .and. & + (material_xs % total - material_xs % absorption) > ZERO) then + cum_dsig = ZERO + associate(mat => materials(p % material)) + do l = 1, mat % n_nuclides + associate (nuc => nuclides(mat % nuclide(l))) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E .and. & + (micro_xs(mat % nuclide(l)) % total & + - micro_xs(mat % nuclide(l)) % absorption) > ZERO) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + cum_dsig = cum_dsig & + + (dsigT - dsigA) * mat % atom_density(l) + end if + end associate + end do + end associate + score = score * (flux_deriv + cum_dsig & + / (material_xs % total - material_xs % absorption)) + else if ( materials(p % material) % id == deriv % diff_material & + .and. (material_xs % total - material_xs % absorption) > ZERO)& + then + dsigT = ZERO + dsigA = ZERO + associate (nuc => nuclides(i_nuclide)) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + end if + end associate + score = score * (flux_deriv + (dsigT - dsigA) & + / (micro_xs(i_nuclide) % total & + - micro_xs(i_nuclide) % absorption)) + else + score = score * flux_deriv + end if + + case (SCORE_ABSORPTION) + if (i_nuclide == -1 .and. & + materials(p % material) % id == deriv % diff_material .and. & + material_xs % absorption > ZERO) then + cum_dsig = ZERO + associate(mat => materials(p % material)) + do l = 1, mat % n_nuclides + associate (nuc => nuclides(mat % nuclide(l))) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E .and. & + micro_xs(mat % nuclide(l)) % absorption > ZERO) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + cum_dsig = cum_dsig + dsigA * mat % atom_density(l) + end if + end associate + end do + end associate + score = score * (flux_deriv & + + cum_dsig / material_xs % absorption) + else if (materials(p % material) % id == deriv % diff_material & + .and. material_xs % absorption > ZERO) then + dsigA = ZERO + associate (nuc => nuclides(i_nuclide)) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + end if + end associate + score = score * (flux_deriv & + + dsigA / micro_xs(i_nuclide) % absorption) + else + score = score * flux_deriv + end if + + case (SCORE_FISSION) + if (i_nuclide == -1 .and. & + materials(p % material) % id == deriv % diff_material .and. & + material_xs % fission > ZERO) then + cum_dsig = ZERO + associate(mat => materials(p % material)) + do l = 1, mat % n_nuclides + associate (nuc => nuclides(mat % nuclide(l))) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E .and. & + micro_xs(mat % nuclide(l)) % fission > ZERO) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + cum_dsig = cum_dsig + dsigF * mat % atom_density(l) + end if + end associate + end do + end associate + score = score * (flux_deriv & + + cum_dsig / material_xs % fission) + else if (materials(p % material) % id == deriv % diff_material & + .and. material_xs % fission > ZERO) then + dsigF = ZERO + associate (nuc => nuclides(i_nuclide)) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + end if + end associate + score = score * (flux_deriv & + + dsigF / micro_xs(i_nuclide) % fission) + else + score = score * flux_deriv + end if + + case (SCORE_NU_FISSION) + if (i_nuclide == -1 .and. & + materials(p % material) % id == deriv % diff_material .and. & + material_xs % nu_fission > ZERO) then + cum_dsig = ZERO + associate(mat => materials(p % material)) + do l = 1, mat % n_nuclides + associate (nuc => nuclides(mat % nuclide(l))) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E .and. & + micro_xs(mat % nuclide(l)) % nu_fission > ZERO) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + cum_dsig = cum_dsig + dsigF * mat % atom_density(l) & + * micro_xs(mat % nuclide(l)) % nu_fission & + / micro_xs(mat % nuclide(l)) % fission + end if + end associate + end do + end associate + score = score * (flux_deriv & + + cum_dsig / material_xs % nu_fission) + else if (materials(p % material) % id == deriv % diff_material & + .and. material_xs % nu_fission > ZERO) then + dsigF = ZERO + associate (nuc => nuclides(i_nuclide)) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + end if + end associate + score = score * (flux_deriv & + + dsigF / micro_xs(i_nuclide) % fission) + else + score = score * flux_deriv + end if + + case default + call fatal_error('Tally derivative not defined for a score on & + &tally ' // trim(to_str(t % id))) + end select + + case default + call fatal_error("Differential tallies are only implemented for & + &analog and collision estimators.") + end select + end select + end associate + end subroutine apply_derivative_to_score + +!=============================================================================== +! SCORE_TRACK_DERIVATIVE Adjust flux derivatives on differential tallies to +! account for a neutron travelling through a perturbed material. +!=============================================================================== + + subroutine score_track_derivative(p, distance) + type(Particle), intent(in) :: p + real(8), intent(in) :: distance ! Neutron flight distance + + 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 + + do i = 1, size(tally_derivs) + associate(deriv => tally_derivs(i)) + select case (deriv % variable) + + case (DIFF_DENSITY) + associate (mat => materials(p % material)) + if (mat % id == deriv % diff_material) then + ! phi is proportional to e^(-Sigma_tot * dist) + ! (1 / phi) * (d_phi / d_rho) = - (d_Sigma_tot / d_rho) * dist + ! (1 / phi) * (d_phi / d_rho) = - Sigma_tot / rho * dist + deriv % flux_deriv = deriv % flux_deriv & + - distance * material_xs % total / mat % density_gpcc + end if + end associate + + case (DIFF_NUCLIDE_DENSITY) + associate (mat => materials(p % material)) + if (mat % id == deriv % diff_material) then + ! phi is proportional to e^(-Sigma_tot * dist) + ! (1 / phi) * (d_phi / d_N) = - (d_Sigma_tot / d_N) * dist + ! (1 / phi) * (d_phi / d_N) = - sigma_tot * dist + deriv % flux_deriv = deriv % flux_deriv & + - 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 .and. & + p % E <= nuc % multipole % end_E) then + ! phi is proportional to e^(-Sigma_tot * dist) + ! (1 / phi) * (d_phi / d_T) = - (d_Sigma_tot / d_T) * dist + ! (1 / phi) * (d_phi / d_T) = - N (d_sigma_tot / d_T) * dist + 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 + end subroutine score_track_derivative + +!=============================================================================== +! SCORE_COLLISION_DERIVATIVE Adjust flux derivatives on differential tallies to +! account for a neutron scattering in the perturbed material. Note that this +! subroutine will be called after absorption events in addition to scattering +! events, but any flux derivatives scored after an absorption will never be +! tallied. This is because the order of operations is +! 1. Particle is moved. +! 2. score_track_derivative is called. +! 3. Collision physics are computed, and the particle is labeled absorbed. +! 4. Analog- and collision-estimated tallies are scored. +! 5. This subroutine is called. +! 6. Particle is killed and no more tallies are scored. +! Hence, it is safe to assume that only derivative of the scattering cross +! section need to be computed here. +!=============================================================================== + + subroutine score_collision_derivative(p) + type(Particle), intent(in) :: p + + integer :: i, j, 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 + + do i = 1, size(tally_derivs) + associate(deriv => tally_derivs(i)) + select case (deriv % variable) + + case (DIFF_DENSITY) + associate (mat => materials(p % material)) + if (mat % id == deriv % diff_material) then + ! phi is proportional to Sigma_s + ! (1 / phi) * (d_phi / d_rho) = (d_Sigma_s / d_rho) / Sigma_s + ! (1 / phi) * (d_phi / d_rho) = 1 / rho + deriv % flux_deriv = deriv % flux_deriv & + + ONE / mat % density_gpcc + end if + end associate + + case (DIFF_NUCLIDE_DENSITY) + associate (mat => materials(p % material)) + if (mat % id == deriv % diff_material & + .and. p % event_nuclide == deriv % diff_nuclide) then + ! Find the index in this material for the diff_nuclide. + do j = 1, mat % n_nuclides + if (mat % nuclide(j) == deriv % diff_nuclide) exit + end do + ! Make sure we found the nuclide. + if (mat % nuclide(j) /= deriv % diff_nuclide) then + call fatal_error("Couldn't find the right nuclide.") + end if + ! phi is proportional to Sigma_s + ! (1 / phi) * (d_phi / d_N) = (d_Sigma_s / d_N) / Sigma_s + ! (1 / phi) * (d_phi / d_N) = sigma_s / Sigma_s + ! (1 / phi) * (d_phi / d_N) = 1 / N + deriv % flux_deriv = deriv % flux_deriv & + + ONE / mat % atom_density(j) + 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 (mat % nuclide(l) == p % event_nuclide .and. & + nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E) then + ! phi is proportional to Sigma_s + ! (1 / phi) * (d_phi / d_T) = (d_Sigma_s / d_T) / Sigma_s + ! (1 / phi) * (d_phi / d_T) = (d_sigma_s / d_T) / sigma_s + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + deriv % flux_deriv = deriv % flux_deriv + (dsigT - dsigA)& + / (micro_xs(mat % nuclide(l)) % total & + - micro_xs(mat % nuclide(l)) % absorption) + ! Note that this is an approximation! The real scattering + ! cross section is Sigma_s(E'->E, uvw'->uvw) = + ! Sigma_s(E') * P(E'->E, uvw'->uvw). We are assuming that + ! d_P(E'->E, uvw'->uvw) / d_T = 0 and only computing + ! d_S(E') / d_T. Using this approximation in the vicinity + ! of low-energy resonances causes errors (~2-5% for PWR + ! pincell eigenvalue derivatives). + end if + end associate + end do + end if + end associate + end select + end associate + end do + end subroutine score_collision_derivative + +!=============================================================================== +! ZERO_FLUX_DERIVS Set the flux derivatives on differential tallies to zero. +!=============================================================================== + + subroutine zero_flux_derivs() + integer :: i + do i = 1, size(tally_derivs) + tally_derivs(i) % flux_deriv = ZERO + end do + end subroutine zero_flux_derivs + !=============================================================================== ! SYNCHRONIZE_TALLIES accumulates the sum of the contributions from each history ! within the batch to a new random variable diff --git a/src/tally_header.F90 b/src/tally_header.F90 index eeda03e762..1d34779aec 100644 --- a/src/tally_header.F90 +++ b/src/tally_header.F90 @@ -10,6 +10,19 @@ module tally_header implicit none +!=============================================================================== +! TALLYDERIVATIVE describes a first-order derivative that can be applied to +! tallies. +!=============================================================================== + + type TallyDerivative + integer :: id + integer :: variable + integer :: diff_material + integer :: diff_nuclide + real(8) :: flux_deriv + end type TallyDerivative + !=============================================================================== ! TALLYOBJECT describes a user-specified tally. The region of phase space to ! tally in is given by the TallyFilters and the results are stored in a @@ -72,6 +85,9 @@ module tally_header integer :: n_triggers = 0 ! # of triggers type(TriggerObject), allocatable :: triggers(:) ! Array of triggers + ! Index for the TallyDerivative for differential tallies. + integer :: deriv = NONE + contains procedure :: write_results_hdf5 procedure :: read_results_hdf5 diff --git a/src/tracking.F90 b/src/tracking.F90 index b62cfc85a0..a45f8e448f 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -14,7 +14,9 @@ module tracking use random_lcg, only: prn use string, only: to_str use tally, only: score_analog_tally, score_tracklength_tally, & - score_collision_tally, score_surface_current + score_collision_tally, score_surface_current, & + score_track_derivative, & + score_collision_derivative, zero_flux_derivs use track_output, only: initialize_particle_track, write_particle_track, & add_particle_track, finalize_particle_track @@ -63,6 +65,9 @@ contains call initialize_particle_track() endif + ! Every particle starts with no accumulated flux derivative. + if (active_tallies % size() > 0) call zero_flux_derivs() + EVENT_LOOP: do ! If the cell hasn't been determined based on the particle's location, ! initiate a search for the current cell. This generally happens at the @@ -137,6 +142,9 @@ contains distance * material_xs % nu_fission end if + ! Score flux derivative accumulators for differential tallies. + if (active_tallies % size() > 0) call score_track_derivative(p, distance) + if (d_collision > d_boundary) then ! ==================================================================== ! PARTICLE CROSSES SURFACE @@ -213,6 +221,9 @@ contains p % coord(j + 1) % uvw = p % coord(j) % uvw end if end do + + ! Score flux derivative accumulators for differential tallies. + if (active_tallies % size() > 0) call score_collision_derivative(p) end if ! Save coordinates for tallying purposes diff --git a/tests/test_diff_tally/inputs_true.dat b/tests/test_diff_tally/inputs_true.dat new file mode 100644 index 0000000000..9a67475791 --- /dev/null +++ b/tests/test_diff_tally/inputs_true.dat @@ -0,0 +1 @@ +df5a0ee7d353fe25344496058cea42e3244a7ab32c13c34c03b4b2f6adf88579f08a71ec0d22d8d264e50e3663717068ff54308e2d16aed8e085c222a8c2fdbd \ No newline at end of file diff --git a/tests/test_diff_tally/results_true.dat b/tests/test_diff_tally/results_true.dat new file mode 100644 index 0000000000..7c4f7f7764 --- /dev/null +++ b/tests/test_diff_tally/results_true.dat @@ -0,0 +1,177 @@ +d_material,d_nuclide,d_variable,score,mean,std. dev. +1,,density,nu-fission,0.0000000e+00,0.0000000e+00 +1,,density,scatter,3.9902949e-02,9.1258428e-03 +1,,density,nu-fission,0.0000000e+00,0.0000000e+00 +1,,density,scatter,9.8213745e-04,9.8213745e-04 +1,,density,nu-fission,2.7945389e-02,2.0999368e-02 +1,,density,scatter,4.0626155e-01,1.7017674e-02 +1,,density,nu-fission,2.4595279e-02,2.1224443e-02 +1,,density,scatter,2.6505962e-03,2.3698970e-03 +1,,density,nu-fission,0.0000000e+00,0.0000000e+00 +1,,density,scatter,-1.2721728e-01,1.8810986e-01 +1,,density,nu-fission,0.0000000e+00,0.0000000e+00 +1,,density,scatter,0.0000000e+00,0.0000000e+00 +1,,density,nu-fission,0.0000000e+00,0.0000000e+00 +1,,density,scatter,4.1950592e-02,7.3680059e-02 +1,,density,nu-fission,0.0000000e+00,0.0000000e+00 +1,,density,scatter,0.0000000e+00,0.0000000e+00 +1,O16,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00 +1,O16,nuclide_density,scatter,-1.9798576e-02,1.9798576e-02 +1,O16,nuclide_density,nu-fission,9.3632921e-01,2.0322286e+00 +1,O16,nuclide_density,scatter,-2.2042881e-01,1.6781767e-01 +1,O16,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00 +1,O16,nuclide_density,scatter,0.0000000e+00,0.0000000e+00 +1,O16,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00 +1,O16,nuclide_density,scatter,0.0000000e+00,0.0000000e+00 +1,U235,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00 +1,U235,nuclide_density,scatter,2.3761276e+00,2.3761276e+00 +1,U235,nuclide_density,nu-fission,7.7441690e+02,8.6460933e+01 +1,U235,nuclide_density,scatter,9.6917103e+01,1.0750582e+01 +1,U235,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00 +1,U235,nuclide_density,scatter,0.0000000e+00,0.0000000e+00 +1,U235,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00 +1,U235,nuclide_density,scatter,0.0000000e+00,0.0000000e+00 +1,,temperature,nu-fission,0.0000000e+00,0.0000000e+00 +1,,temperature,scatter,-2.0607998e-06,2.0607998e-06 +1,,temperature,nu-fission,3.8845881e-06,3.5405083e-05 +1,,temperature,scatter,-6.6772431e-07,2.0750826e-07 +1,,temperature,nu-fission,0.0000000e+00,0.0000000e+00 +1,,temperature,scatter,0.0000000e+00,0.0000000e+00 +1,,temperature,nu-fission,0.0000000e+00,0.0000000e+00 +1,,temperature,scatter,0.0000000e+00,0.0000000e+00 +3,,density,flux,-7.6179774e+00,5.0326447e+00 +3,,density,flux,-1.6242019e+01,6.6231774e+00 +1,,density,flux,-2.2394746e-01,5.4534297e-02 +1,,density,flux,-5.7132854e-02,1.6898134e-01 +1,O16,nuclide_density,flux,-1.3403658e+01,1.3017127e+01 +1,O16,nuclide_density,flux,-1.0499715e+01,2.1684953e+01 +1,U235,nuclide_density,flux,-2.4741168e+03,5.7986812e+01 +1,U235,nuclide_density,flux,-1.9955533e+03,4.3254820e+02 +1,,temperature,flux,1.0601815e-04,3.3076550e-04 +1,,temperature,flux,1.6664450e-04,2.8761112e-04 +3,,density,total,-3.3161585e+00,2.5615932e+00 +3,,density,absorption,-4.3491248e-01,5.1559021e-01 +3,,density,scatter,-2.8812460e+00,2.0540305e+00 +3,,density,fission,-2.3060366e-01,3.1191109e-01 +3,,density,nu-fission,-5.6817321e-01,7.6143389e-01 +3,,density,total,-2.8908007e-01,3.9383184e-01 +3,,density,absorption,-2.5237485e-01,3.6565053e-01 +3,,density,scatter,-3.6705221e-02,2.9275371e-02 +3,,density,fission,-2.1271120e-01,3.0873698e-01 +3,,density,nu-fission,-5.1866339e-01,7.5235533e-01 +3,,density,total,2.7320569e+00,8.7709465e+00 +3,,density,absorption,8.4322466e-02,1.2807898e-01 +3,,density,scatter,2.6477344e+00,8.6432567e+00 +3,,density,fission,0.0000000e+00,0.0000000e+00 +3,,density,nu-fission,0.0000000e+00,0.0000000e+00 +3,,density,total,0.0000000e+00,0.0000000e+00 +3,,density,absorption,0.0000000e+00,0.0000000e+00 +3,,density,scatter,0.0000000e+00,0.0000000e+00 +3,,density,fission,0.0000000e+00,0.0000000e+00 +3,,density,nu-fission,0.0000000e+00,0.0000000e+00 +1,,density,total,4.7523439e-01,2.5652004e-02 +1,,density,absorption,3.1301234e-02,7.6793142e-03 +1,,density,scatter,4.4393315e-01,1.8023252e-02 +1,,density,fission,1.5559355e-02,3.2310294e-03 +1,,density,nu-fission,3.8658109e-02,7.8012233e-03 +1,,density,total,2.2062939e-02,4.3689061e-03 +1,,density,absorption,1.6845626e-02,4.0750601e-03 +1,,density,scatter,5.2173132e-03,2.9663452e-04 +1,,density,fission,1.3503649e-02,3.2642875e-03 +1,,density,nu-fission,3.2945474e-02,7.9507685e-03 +1,,density,total,-9.4735691e-02,2.5584665e-01 +1,,density,absorption,-4.0246397e-03,5.2111770e-03 +1,,density,scatter,-9.0711052e-02,2.5073762e-01 +1,,density,fission,0.0000000e+00,0.0000000e+00 +1,,density,nu-fission,0.0000000e+00,0.0000000e+00 +1,,density,total,0.0000000e+00,0.0000000e+00 +1,,density,absorption,0.0000000e+00,0.0000000e+00 +1,,density,scatter,0.0000000e+00,0.0000000e+00 +1,,density,fission,0.0000000e+00,0.0000000e+00 +1,,density,nu-fission,0.0000000e+00,0.0000000e+00 +1,O16,nuclide_density,total,4.4488215e+01,5.3981603e+00 +1,O16,nuclide_density,absorption,-2.9372611e-01,9.3282813e-01 +1,O16,nuclide_density,scatter,4.4781941e+01,4.5257555e+00 +1,O16,nuclide_density,fission,9.9827872e-02,4.0698927e-01 +1,O16,nuclide_density,nu-fission,2.3342880e-01,9.8749212e-01 +1,O16,nuclide_density,total,4.0914051e-02,6.2301812e-01 +1,O16,nuclide_density,absorption,8.8289354e-02,5.5570922e-01 +1,O16,nuclide_density,scatter,-4.7375303e-02,6.7714071e-02 +1,O16,nuclide_density,fission,1.3678733e-01,4.4006413e-01 +1,O16,nuclide_density,nu-fission,3.3257520e-01,1.0719313e+00 +1,O16,nuclide_density,total,-1.6055918e+01,2.3439335e+01 +1,O16,nuclide_density,absorption,-3.9601231e-01,3.0131622e-01 +1,O16,nuclide_density,scatter,-1.5659906e+01,2.3140407e+01 +1,O16,nuclide_density,fission,0.0000000e+00,0.0000000e+00 +1,O16,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00 +1,O16,nuclide_density,total,0.0000000e+00,0.0000000e+00 +1,O16,nuclide_density,absorption,0.0000000e+00,0.0000000e+00 +1,O16,nuclide_density,scatter,0.0000000e+00,0.0000000e+00 +1,O16,nuclide_density,fission,0.0000000e+00,0.0000000e+00 +1,O16,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00 +1,U235,nuclide_density,total,-6.1444902e+02,6.0111692e+01 +1,U235,nuclide_density,absorption,2.3977134e+02,4.5868405e+01 +1,U235,nuclide_density,scatter,-8.5422036e+02,2.4877920e+01 +1,U235,nuclide_density,fission,3.1319132e+02,3.1465126e+01 +1,U235,nuclide_density,nu-fission,7.6415465e+02,7.6501604e+01 +1,U235,nuclide_density,total,5.0916179e+02,3.5499319e+01 +1,U235,nuclide_density,absorption,3.9477360e+02,3.4326165e+01 +1,U235,nuclide_density,scatter,1.1438818e+02,2.7378217e+00 +1,U235,nuclide_density,fission,3.1360739e+02,3.2243605e+01 +1,U235,nuclide_density,nu-fission,7.6527228e+02,7.8683697e+01 +1,U235,nuclide_density,total,-4.1815711e+03,8.7796611e+02 +1,U235,nuclide_density,absorption,-1.0563226e+02,2.5308627e+01 +1,U235,nuclide_density,scatter,-4.0759388e+03,8.5332185e+02 +1,U235,nuclide_density,fission,0.0000000e+00,0.0000000e+00 +1,U235,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00 +1,U235,nuclide_density,total,0.0000000e+00,0.0000000e+00 +1,U235,nuclide_density,absorption,0.0000000e+00,0.0000000e+00 +1,U235,nuclide_density,scatter,0.0000000e+00,0.0000000e+00 +1,U235,nuclide_density,fission,0.0000000e+00,0.0000000e+00 +1,U235,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00 +1,,temperature,total,2.1367130e-04,1.8632816e-04 +1,,temperature,absorption,6.9456249e-05,3.2199390e-05 +1,,temperature,scatter,1.4421505e-04,1.5754180e-04 +1,,temperature,fission,3.0674981e-06,1.8048377e-05 +1,,temperature,nu-fission,7.4768877e-06,4.3976766e-05 +1,,temperature,total,5.8223116e-06,2.3786474e-05 +1,,temperature,absorption,5.2142698e-06,2.1901911e-05 +1,,temperature,scatter,6.0804185e-07,1.9834890e-06 +1,,temperature,fission,3.0674200e-06,1.8048542e-05 +1,,temperature,nu-fission,7.4767028e-06,4.3977153e-05 +1,,temperature,total,2.1510156e-04,4.6388707e-04 +1,,temperature,absorption,2.2409527e-06,8.2901060e-06 +1,,temperature,scatter,2.1286061e-04,4.5560034e-04 +1,,temperature,fission,0.0000000e+00,0.0000000e+00 +1,,temperature,nu-fission,0.0000000e+00,0.0000000e+00 +1,,temperature,total,0.0000000e+00,0.0000000e+00 +1,,temperature,absorption,0.0000000e+00,0.0000000e+00 +1,,temperature,scatter,0.0000000e+00,0.0000000e+00 +1,,temperature,fission,0.0000000e+00,0.0000000e+00 +1,,temperature,nu-fission,0.0000000e+00,0.0000000e+00 +3,,density,absorption,-1.6517201e-01,3.0984738e-01 +3,,density,absorption,8.0402344e-03,1.4689335e-01 +1,,density,absorption,2.9069882e-02,2.2139609e-03 +1,,density,absorption,-9.4690065e-03,8.6605392e-03 +1,O16,nuclide_density,absorption,7.6911962e-01,4.1945687e-01 +1,O16,nuclide_density,absorption,-6.3724795e-01,6.6341488e-01 +1,U235,nuclide_density,absorption,1.4109543e+02,1.8518890e+01 +1,U235,nuclide_density,absorption,-1.2052822e+02,3.2084002e+01 +1,,temperature,absorption,3.9995404e-05,2.1703384e-05 +1,,temperature,absorption,2.7274361e-06,8.9339257e-06 +3,,density,nu-fission,0.0000000e+00,0.0000000e+00 +3,,density,scatter,-8.7934551e-01,1.0210652e+00 +3,,density,nu-fission,0.0000000e+00,0.0000000e+00 +3,,density,scatter,-3.5236516e-03,3.5236516e-03 +3,,density,nu-fission,7.2953012e-02,2.9725863e-01 +3,,density,scatter,-2.2716410e+00,1.2546347e+00 +3,,density,nu-fission,8.9171481e-02,3.0634856e-01 +3,,density,scatter,-1.0151747e-03,9.6296775e-03 +3,,density,nu-fission,0.0000000e+00,0.0000000e+00 +3,,density,scatter,2.3163923e+00,4.8766690e+00 +3,,density,nu-fission,0.0000000e+00,0.0000000e+00 +3,,density,scatter,0.0000000e+00,0.0000000e+00 +3,,density,nu-fission,0.0000000e+00,0.0000000e+00 +3,,density,scatter,4.0762434e-01,3.8504141e+00 +3,,density,nu-fission,0.0000000e+00,0.0000000e+00 +3,,density,scatter,0.0000000e+00,0.0000000e+00 diff --git a/tests/test_diff_tally/test_diff_tally.py b/tests/test_diff_tally/test_diff_tally.py new file mode 100644 index 0000000000..36287a66ac --- /dev/null +++ b/tests/test_diff_tally/test_diff_tally.py @@ -0,0 +1,140 @@ +#!/usr/bin/env python + +import glob +import os +import sys + +import pandas as pd + +sys.path.insert(0, os.pardir) +from testing_harness import PyAPITestHarness +import openmc + +class DiffTallyTestHarness(PyAPITestHarness): + def _build_inputs(self): + # Build default materials/geometry + self._input_set.build_default_materials_and_geometry() + + # Set settings explicitly + self._input_set.settings.batches = 3 + self._input_set.settings.inactive = 0 + self._input_set.settings.particles = 100 + self._input_set.settings.source = openmc.Source(space=openmc.stats.Box( + [-160, -160, -183], [160, 160, 183])) + self._input_set.settings.temperature['multipole'] = True + + self._input_set.tallies = openmc.Tallies() + + filt_mats = openmc.MaterialFilter((1, 3)) + filt_eout = openmc.EnergyoutFilter((0.0, 0.625, 20.0e6)) + + # We want density derivatives for both water and fuel to get coverage + # for both fissile and non-fissile materials. + d1 = openmc.TallyDerivative(derivative_id=1) + d1.variable = 'density' + d1.material = 3 + d2 = openmc.TallyDerivative(derivative_id=2) + d2.variable = 'density' + d2.material = 1 + + # O-16 is a good nuclide to test against because it is present in both + # water and fuel. Some routines need to recognize that they have the + # perturbed nuclide but not the perturbed material. + d3 = openmc.TallyDerivative(derivative_id=3) + d3.variable = 'nuclide_density' + d3.material = 1 + d3.nuclide = 'O16' + + # A fissile nuclide, just for good measure. + d4 = openmc.TallyDerivative(derivative_id=4) + d4.variable = 'nuclide_density' + d4.material = 1 + d4.nuclide = 'U235' + + # Temperature derivatives. + d5 = openmc.TallyDerivative(derivative_id=5) + d5.variable = 'temperature' + d5.material = 1 + + derivs = [d1, d2, d3, d4, d5] + + # Cover the flux score. + for i in range(5): + t = openmc.Tally() + t.add_score('flux') + t.add_filter(filt_mats) + t.derivative = derivs[i] + self._input_set.tallies.append(t) + + # Cover supported scores with a collision estimator. + for i in range(5): + t = openmc.Tally() + t.add_score('total') + t.add_score('absorption') + t.add_score('scatter') + t.add_score('fission') + t.add_score('nu-fission') + t.add_filter(filt_mats) + t.add_nuclide('total') + t.add_nuclide('U235') + t.derivative = derivs[i] + self._input_set.tallies.append(t) + + # Cover an analog estimator. + for i in range(5): + t = openmc.Tally() + t.add_score('absorption') + t.add_filter(filt_mats) + t.estimator = 'analog' + t.derivative = derivs[i] + self._input_set.tallies.append(t) + + # Energyout filter and total nuclide for the density derivatives. + for i in range(2): + t = openmc.Tally() + t.add_score('nu-fission') + t.add_score('scatter') + t.add_filter(filt_mats) + t.add_filter(filt_eout) + t.add_nuclide('total') + t.add_nuclide('U235') + t.derivative = derivs[i] + self._input_set.tallies.append(t) + + # Energyout filter without total nuclide for other derivatives. + for i in range(2, 5): + t = openmc.Tally() + t.add_score('nu-fission') + t.add_score('scatter') + t.add_filter(filt_mats) + t.add_filter(filt_eout) + t.add_nuclide('U235') + t.derivative = derivs[i] + self._input_set.tallies.append(t) + + self._input_set.export() + + def _get_results(self): + # Read the statepoint and summary files. + statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name))[0] + sp = openmc.StatePoint(statepoint) + + # Extract the tally data as a Pandas DataFrame. + df = pd.DataFrame() + for t in sp.tallies.values(): + df = df.append(t.get_pandas_dataframe(), ignore_index=True) + + # Extract the relevant data as a CSV string. + cols = ('d_material', 'd_nuclide', 'd_variable', 'score', 'mean', + 'std. dev.') + return df.to_csv(None, columns=cols, index=False, float_format='%.7e') + + def _cleanup(self): + super(DiffTallyTestHarness, self)._cleanup() + f = os.path.join(os.getcwd(), 'tallies.xml') + if os.path.exists(f): os.remove(f) + + +if __name__ == '__main__': + harness = DiffTallyTestHarness('statepoint.3.h5', True) + harness.main()