From ff3d03427536a78e5967e856c033a513bebe4c9c Mon Sep 17 00:00:00 2001 From: samuelshaner Date: Thu, 15 Dec 2016 20:23:07 +0000 Subject: [PATCH] added tracklength estimator capability for decay-rate tally and fixed other issues --- openmc/material.py | 30 +++-- openmc/mgxs/library.py | 9 ++ openmc/mgxs/mdgxs.py | 73 ++++++++++- openmc/mgxs/mgxs.py | 118 ++++++++++++++++++ openmc/mgxs_library.py | 72 +++++++++-- openmc/tallies.py | 11 +- src/input_xml.F90 | 7 -- src/nuclide_header.F90 | 3 +- src/tally.F90 | 273 ++++++++++++++++++++++++++++++----------- 9 files changed, 486 insertions(+), 110 deletions(-) diff --git a/openmc/material.py b/openmc/material.py index a7a05989b..f5c5e80cc 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -129,7 +129,7 @@ class Material(object): string = 'Material\n' string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id) string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name) - string += '{0: <16}{1}{2}\n'.format('\Temperature', '=\t', + string += '{0: <16}{1}{2}\n'.format('\tTemperature', '=\t', self._temperature) string += '{0: <16}{1}{2}'.format('\tDensity', '=\t', self._density) @@ -251,7 +251,7 @@ class Material(object): Parameters ---------- - units : {'g/cm3', 'g/cc', 'km/cm3', 'atom/b-cm', 'atom/cm3', 'sum', 'macro'} + units : {'g/cm3', 'g/cc', 'kg/cm3', 'atom/b-cm', 'atom/cm3', 'sum', 'macro'} Physical units of density. density : float, optional Value of the density. Must be specified unless units is given as @@ -641,37 +641,45 @@ class Material(object): nucs = [] nuc_densities = [] nuc_density_types = [] + for nuclide in nuclides.items(): nuc, nuc_density, nuc_density_type = nuclide[1] nucs.append(nuc) nuc_densities.append(nuc_density) nuc_density_types.append(nuc_density_type) + nucs = np.array(nucs) + nuc_densities = np.array(nuc_densities) + nuc_density_types = np.array(nuc_density_types) + if sum_density: density = np.sum(nuc_densities) + percent_in_atom = np.all(nuc_density_types == 'ao') density_in_atom = density > 0. sum_percent = 0. awrs = [] - for n, nuclide in enumerate(nuclides.items()): + for nuclide in nuclides.items(): awr = openmc.data.atomic_mass(nuclide[0]) if awr is not None: - awrs.append(awr / openmc.data.NEUTRON_MASS) + awrs.append(awr) else: raise ValueError(nuclide[0] + " is invalid") + # Convert the weight amounts to atomic amounts + if not percent_in_atom: + for n, nuc in enumerate(nucs): + nuc_densities[n] *= self.average_molar_mass / awrs[n] + # Now that we have the awr, lets finish calculating densities sum_percent = np.sum(nuc_densities) nuc_densities = nuc_densities / sum_percent + if not density_in_atom: - sum_percent = 0. - for n, nuc in enumerate(nucs): - x = nuc_densities[n] - sum_percent += x * awrs[n] - sum_percent = 1. / sum_percent - density = -density * sum_percent * \ - openmc.data.AVOGADRO / openmc.data.NEUTRON_MASS * 1.E-24 + density = -density / self.average_molar_mass * 1.E-24 \ + * openmc.data.AVOGADRO + nuc_densities = density * nuc_densities nuclides = OrderedDict() diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index 7e0289064..a7440a5c5 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -524,6 +524,15 @@ class Library(object): for domain in self.domains: for mgxs_type in self.mgxs_types: mgxs = self.get_mgxs(domain, mgxs_type) + + if mgxs_type in openmc.mgxs.MDGXS_TYPES: + if self.num_delayed_groups == 0: + mgxs.delayed_groups = None + else: + delayed_groups \ + = list(range(1,self.num_delayed_groups+1)) + mgxs.delayed_groups = delayed_groups + for tally in mgxs.tallies.values(): tallies_file.append(tally, merge=merge) diff --git a/openmc/mgxs/mdgxs.py b/openmc/mgxs/mdgxs.py index 5b2f451c2..1da7f99e2 100644 --- a/openmc/mgxs/mdgxs.py +++ b/openmc/mgxs/mdgxs.py @@ -966,6 +966,7 @@ class ChiDelayed(MDGXS): super(ChiDelayed, self).__init__(domain, domain_type, energy_groups, delayed_groups, by_nuclide, name) self._rxn_type = 'chi-delayed' + self._estimator = 'analog' @property def scores(self): @@ -987,10 +988,6 @@ class ChiDelayed(MDGXS): def tally_keys(self): return ['delayed-nu-fission-in', 'delayed-nu-fission-out'] - @property - def estimator(self): - return 'analog' - @property def rxn_rate_tally(self): if self._rxn_rate_tally is None: @@ -1018,6 +1015,29 @@ class ChiDelayed(MDGXS): return self._xs_tally + def get_homogenized_mgxs(self, other_mgxs): + """Construct a homogenized mgxs with other mgxs objects. + + Parameters + ---------- + other_mgxs : openmc.mgxs.MGXS or Iterable of openmc.mgxs.MGXS + The MGXS to homogenize with this one. + + Returns + ------- + openmc.mgxs.MGXS + A new homogenized MGXS + + Raises + ------ + ValueError + If the other_mgxs are of a different type. + + """ + + return self._get_homogenized_mgxs(other_mgxs, 'delayed-nu-fission-in') + + def get_slice(self, nuclides=[], groups=[], delayed_groups=[]): """Build a sliced ChiDelayed for the specified nuclides and energy groups. @@ -1586,6 +1606,28 @@ class Beta(MDGXS): return self._xs_tally + def get_homogenized_mgxs(self, other_mgxs): + """Construct a homogenized mgxs with other mgxs objects. + + Parameters + ---------- + other_mgxs : openmc.mgxs.MGXS or Iterable of openmc.mgxs.MGXS + The MGXS to homogenize with this one. + + Returns + ------- + openmc.mgxs.MGXS + A new homogenized MGXS + + Raises + ------ + ValueError + If the other_mgxs are of a different type. + + """ + + return self._get_homogenized_mgxs(other_mgxs, 'nu-fission') + class DecayRate(MDGXS): r"""The decay rate for delayed neutron precursors. @@ -1703,7 +1745,6 @@ class DecayRate(MDGXS): super(DecayRate, self).__init__(domain, domain_type, energy_groups, delayed_groups, by_nuclide, name) self._rxn_type = 'decay-rate' - self._estimator = 'analog' @property def scores(self): @@ -1738,6 +1779,28 @@ class DecayRate(MDGXS): return self._xs_tally + def get_homogenized_mgxs(self, other_mgxs): + """Construct a homogenized mgxs with other mgxs objects. + + Parameters + ---------- + other_mgxs : openmc.mgxs.MGXS or Iterable of openmc.mgxs.MGXS + The MGXS to homogenize with this one. + + Returns + ------- + openmc.mgxs.MGXS + A new homogenized MGXS + + Raises + ------ + ValueError + If the other_mgxs are of a different type. + + """ + + return self._get_homogenized_mgxs(other_mgxs, 'delayed-nu-fission') + @add_metaclass(ABCMeta) class MatrixMDGXS(MDGXS): diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index 3ad67484a..f8b546a2c 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -996,6 +996,102 @@ class MGXS(object): avg_xs.sparse = self.sparse return avg_xs + def _get_homogenized_mgxs(self, other_mgxs, denom_score='flux'): + """Construct a homogenized mgxs with other mgxs objects. + + Parameters + ---------- + other_mgxs : openmc.mgxs.MGXS or Iterable of openmc.mgxs.MGXS + The MGXS to homogenize with this one. + denom_score : str + The denominator score in the denominator of computing the MGXS. + + Returns + ------- + openmc.mgxs.MGXS + A new homogenized MGXS + + Raises + ------ + ValueError + If the other_mgxs are of a different type. + + """ + + # Check type of denom score + cv.check_type('denom_score', denom_score, str) + + # Construct a collection of the subdomain filter bins to average across + if isinstance(other_mgxs, openmc.mgxs.MGXS): + other_mgxs = [other_mgxs] + + cv.check_iterable_type('other_mgxs', other_mgxs, openmc.mgxs.MGXS) + for mgxs in other_mgxs: + if mgxs.rxn_type != self.rxn_type: + msg = 'Not able to homogenize two MGXS with different rxn types' + raise ValueError(msg) + + # Clone this MGXS to initialize the homogenized version + homogenized_mgxs = copy.deepcopy(self) + homogenized_mgxs._derived = True + name = '{} + '.format(self.domain.name) + + # Get the domain filter + filter_type = _DOMAIN_TO_FILTER[self.domain_type] + self_filter = self.rxn_rate_tally.find_filter(filter_type) + + # Get the rxn rate and denom tallies + rxn_rate_tally = self.rxn_rate_tally + denom_tally = self.tallies[denom_score] + + for mgxs in other_mgxs: + + # Swap the domain filter bins for the other mgxs rxn rate tally + other_rxn_rate_tally = copy.deepcopy(mgxs.rxn_rate_tally) + other_filter = other_rxn_rate_tally.find_filter(filter_type) + other_filter._bins = self_filter._bins + + # Swap the domain filter bins for the denom tally + other_denom_tally = copy.deepcopy(mgxs.tallies[denom_score]) + other_filter = other_denom_tally.find_filter(filter_type) + other_filter._bins = self_filter._bins + + # Add the rxn rate and denom tallies + rxn_rate_tally += other_rxn_rate_tally + denom_tally += other_denom_tally + + # Update the name for the homogenzied MGXS + name += '{} + '.format(mgxs.domain.name) + + # Set the properties of the homogenized MGXS + homogenized_mgxs._rxn_rate_tally = rxn_rate_tally + homogenized_mgxs.tallies[denom_score] = denom_tally + homogenized_mgxs._domain.name = name[:-3] + + return homogenized_mgxs + + def get_homogenized_mgxs(self, other_mgxs): + """Construct a homogenized mgxs with other mgxs objects. + + Parameters + ---------- + other_mgxs : openmc.mgxs.MGXS or Iterable of openmc.mgxs.MGXS + The MGXS to homogenize with this one. + + Returns + ------- + openmc.mgxs.MGXS + A new homogenized MGXS + + Raises + ------ + ValueError + If the other_mgxs are of a different type. + + """ + + return self._get_homogenized_mgxs(other_mgxs, 'flux') + def get_slice(self, nuclides=[], groups=[]): """Build a sliced MGXS for the specified nuclides and energy groups. @@ -4561,6 +4657,28 @@ class Chi(MGXS): return self._xs_tally + def get_homogenized_mgxs(self, other_mgxs): + """Construct a homogenized mgxs with other mgxs objects. + + Parameters + ---------- + other_mgxs : openmc.mgxs.MGXS or Iterable of openmc.mgxs.MGXS + The MGXS to homogenize with this one. + + Returns + ------- + openmc.mgxs.MGXS + A new homogenized MGXS + + Raises + ------ + ValueError + If the other_mgxs are of a different type. + + """ + + return self._get_homogenized_mgxs(other_mgxs, 'nu-fission-in') + def get_slice(self, nuclides=[], groups=[]): """Build a sliced Chi for the specified nuclides and energy groups. diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py index 331d9ed12..e92c7eee3 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -15,8 +15,8 @@ from openmc.checkvalue import check_type, check_value, check_greater_than, \ # Supported incoming particle MGXS angular treatment representations _REPRESENTATIONS = ['isotropic', 'angle'] _SCATTER_TYPES = ['tabular', 'legendre', 'histogram'] -_XS_SHAPES = ["[G][G'][Order]", "[G]", "[G']", "[G][G']", "[DG]", "[G][DG]", - "[G'][DG]", "[G][G'][DG]"] +_XS_SHAPES = ["[G][G'][Order]", "[G]", "[G']", "[G][G']", "[DG]", "[DG][G]", + "[DG][G']", "[DG][G][G']"] class XSdata(object): @@ -145,11 +145,11 @@ class XSdata(object): [DG]: beta, decay_rate - [G][DG]: delayed_nu_fission, beta, decay_rate + [DG][G]: delayed_nu_fission, beta, decay_rate - [G'][DG]: chi_delayed + [DG][G']: chi_delayed - [G][G'][DG]: delayed_nu_fission + [DG][G][G']: delayed_nu_fission """ @@ -296,11 +296,11 @@ class XSdata(object): self._xs_shapes["[G][G']"] = (self.energy_groups.num_groups, self.energy_groups.num_groups) self._xs_shapes["[DG]"] = (self.num_delayed_groups,) - self._xs_shapes["[G][DG]"] = (self.energy_groups.num_groups, + self._xs_shapes["[DG][G]"] = (self.energy_groups.num_groups, self.num_delayed_groups) - self._xs_shapes["[G'][DG]"] = (self.energy_groups.num_groups, + self._xs_shapes["[DG'][G']"] = (self.energy_groups.num_groups, self.num_delayed_groups) - self._xs_shapes["[G][G'][DG]"] = (self.energy_groups.num_groups, + self._xs_shapes["[DG][G][G']"] = (self.energy_groups.num_groups, self.energy_groups.num_groups, self.num_delayed_groups) @@ -634,7 +634,7 @@ class XSdata(object): """ # Get the accepted shapes for this xs - shapes = [self.xs_shapes["[G']"], self.xs_shapes["[G'][DG]"]] + shapes = [self.xs_shapes["[G']"], self.xs_shapes["[DG][G']"]] # Convert to a numpy array so we can easily get the shape for checking chi_delayed = np.asarray(chi_delayed) @@ -664,7 +664,7 @@ class XSdata(object): """ # Get the accepted shapes for this xs - shapes = [self.xs_shapes["[DG]"], self.xs_shapes["[G][DG]"]] + shapes = [self.xs_shapes["[DG]"], self.xs_shapes["[DG][G]"]] # Convert to a numpy array so we can easily get the shape for checking beta = np.asarray(beta) @@ -694,7 +694,7 @@ class XSdata(object): """ # Get the accepted shapes for this xs - shapes = [self.xs_shapes["[DG]"], self.xs_shapes["[G][DG]"]] + shapes = [self.xs_shapes["[DG]"], self.xs_shapes["[DG][G]"]] # Convert to a numpy array so we can easily get the shape for checking decay_rate = np.asarray(decay_rate) @@ -856,7 +856,7 @@ class XSdata(object): """ # Get the accepted shapes for this xs - shapes = [self.xs_shapes["[G][DG]"], self.xs_shapes["[G][G'][DG]"]] + shapes = [self.xs_shapes["[G][DG]"], self.xs_shapes["[DG][G][G']"]] # Convert to a numpy array so we can easily get the shape for checking delayed_nu_fission = np.asarray(delayed_nu_fission) @@ -1644,6 +1644,54 @@ class XSdata(object): self._multiplicity_matrix[i] = \ np.nan_to_num(self._multiplicity_matrix[i]) + def set_inverse_velocity_mgxs(self, inverse_velocity, temperature=294., + nuclide='total', xs_type='macro', + subdomain=None): + """This method allows for an openmc.mgxs.InverseVelocity + to be used to set the absorption cross section for this XSdata object. + + Parameters + ---------- + inverse_velocity : openmc.mgxs.InverseVelocity + MGXS Object containing the inverse velocity for the domain of + interest. + temperature : float + Temperature (in Kelvin) of the data. Defaults to room temperature + (294K). + nuclide : str + Individual nuclide (or 'total' if obtaining material-wise data) + to gather data for. Defaults to 'total'. + xs_type: {'macro', 'micro'} + Provide the macro or micro cross section in units of cm^-1 or + barns. Defaults to 'macro'. + subdomain : iterable of int + If the MGXS contains a mesh domain type, the subdomain parameter + specifies which mesh cell (i.e., [i, j, k] index) to use. + + See also + -------- + openmc.mgxs.Library.create_mg_library() + openmc.mgxs.Library.get_xsdata + + """ + + check_type('inverse_velocity', inverse_velocity, openmc.mgxs.InverseVelocity) + check_value('energy_groups', inverse_velocity.energy_groups, + [self.energy_groups]) + check_value('domain_type', inverse_velocity.domain_type, + openmc.mgxs.DOMAIN_TYPES) + check_type('temperature', temperature, Real) + check_value('temperature', temperature, self.temperatures) + + i = np.where(self.temperatures == temperature)[0][0] + if self.representation == 'isotropic': + self._inverse_velocity[i] = inverse_velocity.get_xs\ + (nuclides=nuclide, xs_type=xs_type, + subdomains=subdomain) + elif self.representation == 'angle': + msg = 'Angular-Dependent MGXS have not yet been implemented' + raise ValueError(msg) + def to_hdf5(self, file): """Write XSdata to an HDF5 file diff --git a/openmc/tallies.py b/openmc/tallies.py index af57a1ffd..20e03129b 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -705,7 +705,7 @@ class Tally(object): """ - # Two tallys must have the same number of filters + # Two tallies must have the same number of filters if len(self.filters) != len(other.filters): return False @@ -3486,9 +3486,14 @@ class Tallies(cv.CheckedList): for d in derivs: root_element.append(d.to_xml_element()) - def export_to_xml(self): + def export_to_xml(self, path='tallies.xml'): """Create a tallies.xml file that can be used for a simulation. + Parameters + ---------- + path : str + Path to file to write. Defaults to 'tallies.xml'. + """ root_element = ET.Element("tallies") @@ -3501,5 +3506,5 @@ class Tallies(cv.CheckedList): # Write the XML Tree to the tallies.xml file tree = ET.ElementTree(root_element) - tree.write("tallies.xml", xml_declaration=True, + tree.write(path, xml_declaration=True, encoding='utf-8', method="xml") diff --git a/src/input_xml.F90 b/src/input_xml.F90 index cdb9c253f..c39ef0412 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3668,13 +3668,6 @@ contains end if case ('decay-rate') t % score_bins(j) = SCORE_DECAY_RATE - - ! Set tally estimator to analog for CE mode - ! (MG mode has all data available without a collision being - ! necessary) - if (run_CE) then - t % estimator = ESTIMATOR_ANALOG - end if case ('delayed-nu-fission') t % score_bins(j) = SCORE_DELAYED_NU_FISSION if (t % find_filter(FILTER_ENERGYOUT) > 0) then diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index e57d7da62..a329b5bcc 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -621,7 +621,8 @@ module nuclide_header case (EMISSION_DELAYED) if (this % n_precursor > 0) then - if (present(group)) then + if (present(group) .and. group < & + size(this % reactions(this % index_fission(1)) % products)) then ! If delayed group specified, determine yield immediately associate(p => this % reactions(this % index_fission(1)) % products(1 + group)) nu = p % yield % evaluate(E) diff --git a/src/tally.F90 b/src/tally.F90 index e3fdd1ab5..8d06dc9ac 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -679,8 +679,8 @@ contains end if end if - case (SCORE_DECAY_RATE) + ! make sure the correct energy is used if (t % estimator == ESTIMATOR_TRACKLENGTH) then E = p % E @@ -691,12 +691,134 @@ contains ! Set the delayedgroup filter index dg_filter = t % find_filter(FILTER_DELAYEDGROUP) - if (survival_biasing) then - ! No fission events occur if survival biasing is on -- need to - ! calculate fraction of absorptions that would have resulted in - ! delayed-nu-fission - if (micro_xs(p % event_nuclide) % absorption > ZERO .and. & - nuclides(p % event_nuclide) % fissionable) then + if (t % estimator == ESTIMATOR_ANALOG) then + if (survival_biasing) then + ! No fission events occur if survival biasing is on -- need to + ! calculate fraction of absorptions that would have resulted in + ! delayed-nu-fission + if (micro_xs(p % event_nuclide) % absorption > ZERO .and. & + nuclides(p % event_nuclide) % fissionable) then + + ! Check if the delayed group filter is present + if (dg_filter > 0) then + select type(filt => t % filters(dg_filter) % obj) + type is (DelayedGroupFilter) + + ! Loop over all delayed group bins and tally to them + ! individually + do d_bin = 1, filt % n_bins + + ! Get the delayed group for this bin + d = filt % groups(d_bin) + + ! Compute the yield for this delayed group + yield = nuclides(p % event_nuclide) & + % nu(E, EMISSION_DELAYED, d) + + associate (rxn => nuclides(p % event_nuclide) % & + reactions(nuclides(p % event_nuclide) % index_fission(1))) + + ! Compute the score + score = p % absorb_wgt * yield * & + micro_xs(p % event_nuclide) % fission & + / micro_xs(p % event_nuclide) % absorption & + * rxn % products(1 + d) % decay_rate + end associate + + ! Tally to bin + call score_fission_delayed_dg(t, d_bin, score, score_index) + end do + cycle SCORE_LOOP + end select + else + + ! If the delayed group filter is not present, compute the score + ! by accumulating the absorbed weight times the decay rate times + ! the fraction of the delayed-nu-fission xs to the absorption xs + ! for all delayed groups. + score = ZERO + + associate (rxn => nuclides(p % event_nuclide) % & + reactions(nuclides(p % event_nuclide) % index_fission(1))) + + ! We need to be careful not to overshoot the number of delayed + ! groups since this could cause the range of the rxn % products + ! array to be exceeded. Hence, we use the size of this array + ! and not the MAX_DELAYED_GROUPS constant for this loop. + do d = 1, size(rxn % products) - 2 + + score = score + rxn % products(1 + d) % decay_rate * & + p % absorb_wgt * micro_xs(p % event_nuclide) % fission *& + nuclides(p % event_nuclide) % nu(E, EMISSION_DELAYED, d)& + / micro_xs(p % event_nuclide) % absorption + end do + end associate + end if + end if + else + + ! Skip any non-fission events + if (.not. p % fission) cycle SCORE_LOOP + ! If there is no outgoing energy filter, than we only need to + ! score to one bin. For the score to be 'analog', we need to + ! score the number of particles that were banked in the fission + ! bank. Since this was weighted by 1/keff, we multiply by keff + ! to get the proper score. Loop over the neutrons produced from + ! fission and check which ones are delayed. If a delayed neutron is + ! encountered, add its contribution to the fission bank to the + ! score. + + score = ZERO + + ! loop over number of particles banked + do k = 1, p % n_bank + + ! get the delayed group + g = fission_bank(n_bank - p % n_bank + k) % delayed_group + + ! Case for tallying delayed emissions + if (g /= 0) then + + ! Accumulate the decay rate times delayed nu fission score + associate (rxn => nuclides(p % event_nuclide) % & + reactions(nuclides(p % event_nuclide) % index_fission(1))) + + ! determine score based on bank site weight and keff. + score = score + keff * fission_bank(n_bank - p % n_bank + k) & + % wgt * rxn % products(1 + g) % decay_rate + end associate + + ! if the delayed group filter is present, tally to corresponding + ! delayed group bin if it exists + if (dg_filter > 0) then + + ! declare the delayed group filter type + select type(filt => t % filters(dg_filter) % obj) + type is (DelayedGroupFilter) + + ! loop over delayed group bins until the corresponding bin is + ! found + do d_bin = 1, filt % n_bins + d = filt % groups(d_bin) + + ! check whether the delayed group of the particle is equal to + ! the delayed group of this bin + if (d == g) then + call score_fission_delayed_dg(t, d_bin, score, score_index) + end if + end do + end select + + ! Reset the score to zero + score = ZERO + end if + end if + end do + end if + else + + ! Check if tally is on a single nuclide + if (i_nuclide > 0) then ! Check if the delayed group filter is present if (dg_filter > 0) then @@ -711,17 +833,14 @@ contains d = filt % groups(d_bin) ! Compute the yield for this delayed group - yield = nuclides(p % event_nuclide) & - % nu(E, EMISSION_DELAYED, d) + yield = nuclides(i_nuclide) % nu(E, EMISSION_DELAYED, d) - associate (rxn => nuclides(p % event_nuclide) % & - reactions(nuclides(p % event_nuclide) % index_fission(1))) + associate (rxn => nuclides(i_nuclide) % & + reactions(nuclides(i_nuclide) % index_fission(1))) - ! Compute the score - score = p % absorb_wgt * yield * & - micro_xs(p % event_nuclide) % fission & - / micro_xs(p % event_nuclide) % absorption & - * rxn % products(1 + d) % decay_rate + ! Compute the score and tally to bin + score = micro_xs(i_nuclide) % fission * yield * flux * & + atom_density * rxn % products(1 + d) % decay_rate end associate ! Tally to bin @@ -746,78 +865,90 @@ contains ! and not the MAX_DELAYED_GROUPS constant for this loop. do d = 1, size(rxn % products) - 2 - score = score + rxn % products(1 + d) % decay_rate * & - p % absorb_wgt * micro_xs(p % event_nuclide) % fission *& - nuclides(p % event_nuclide) % nu(E, EMISSION_DELAYED, d)& - / micro_xs(p % event_nuclide) % absorption + score = score + micro_xs(i_nuclide) % fission * flux * & + nuclides(i_nuclide) % nu(E, EMISSION_DELAYED) * & + atom_density * rxn % products(1 + d) % decay_rate end do end associate end if - end if - else - ! Skip any non-fission events - if (.not. p % fission) cycle SCORE_LOOP - ! If there is no outgoing energy filter, than we only need to - ! score to one bin. For the score to be 'analog', we need to - ! score the number of particles that were banked in the fission - ! bank. Since this was weighted by 1/keff, we multiply by keff - ! to get the proper score. Loop over the neutrons produced from - ! fission and check which ones are delayed. If a delayed neutron is - ! encountered, add its contribution to the fission bank to the - ! score. + ! Tally is on total nuclides + else - score = ZERO + ! Check if the delayed group filter is present + if (dg_filter > 0) then + select type(filt => t % filters(dg_filter) % obj) + type is (DelayedGroupFilter) - ! loop over number of particles banked - do k = 1, p % n_bank + ! Loop over all nuclides in the current material + do l = 1, materials(p % material) % n_nuclides - ! get the delayed group - g = fission_bank(n_bank - p % n_bank + k) % delayed_group + ! Get atom density + atom_density_ = materials(p % material) % atom_density(l) - ! Case for tallying delayed emissions - if (g /= 0) then + ! Get index in nuclides array + i_nuc = materials(p % material) % nuclide(l) - ! Accumulate the decay rate times delayed nu fission score - associate (rxn => nuclides(p % event_nuclide) % & - reactions(nuclides(p % event_nuclide) % index_fission(1))) + if (nuclides(i_nuc) % fissionable) then - ! determine score based on bank site weight and keff. - score = score + keff * fission_bank(n_bank - p % n_bank + k) & - % wgt * rxn % products(1 + g) % decay_rate - end associate + ! Loop over all delayed group bins and tally to them + ! individually + do d_bin = 1, filt % n_bins - ! if the delayed group filter is present, tally to corresponding - ! delayed group bin if it exists - if (dg_filter > 0) then + ! Get the delayed group for this bin + d = filt % groups(d_bin) - ! declare the delayed group filter type - select type(filt => t % filters(dg_filter) % obj) - type is (DelayedGroupFilter) + ! Get the yield for the desired nuclide and delayed group + yield = nuclides(i_nuc) % nu(E, EMISSION_DELAYED, d) - ! loop over delayed group bins until the corresponding bin is - ! found - do d_bin = 1, filt % n_bins - d = filt % groups(d_bin) + associate (rxn => nuclides(i_nuc) % & + reactions(nuclides(i_nuc) % index_fission(1))) - ! check whether the delayed group of the particle is equal to - ! the delayed group of this bin - if (d == g) then + ! Compute the score + score = micro_xs(i_nuc) % fission * yield * flux * & + atom_density_ * rxn % products(1 + d) % decay_rate + end associate + + ! Tally to bin call score_fission_delayed_dg(t, d_bin, score, score_index) - end if - end do - end select + end do + end if + end do + cycle SCORE_LOOP + end select + else - ! Reset the score to zero - score = ZERO - end if + score = ZERO + + ! Loop over all nuclides in the current material + do l = 1, materials(p % material) % n_nuclides + + ! Get atom density + atom_density_ = materials(p % material) % atom_density(l) + + ! Get index in nuclides array + i_nuc = materials(p % material) % nuclide(l) + + if (nuclides(i_nuc) % fissionable) then + + associate (rxn => nuclides(i_nuc) % & + reactions(nuclides(i_nuc) % index_fission(1))) + + ! We need to be careful not to overshoot the number of delayed + ! groups since this could cause the range of the rxn % products + ! array to be exceeded. Hence, we use the size of this array + ! and not the MAX_DELAYED_GROUPS constant for this loop. + do d = 1, size(rxn % products) - 2 + + ! Accumulate the contribution from each nuclide + score = score + micro_xs(i_nuc) % fission * nuclides(i_nuc) %& + nu(E, EMISSION_DELAYED) * atom_density_ * flux * & + rxn % products(1 + d) % decay_rate + end do + end associate + end if + end do end if - end do - - ! If the delayed group filter is present, cycle because the - ! score_fission_delayed_dg(...) has already tallied the score - if (dg_filter > 0) then - cycle SCORE_LOOP end if end if