diff --git a/docs/source/io_formats/nuclear_data.rst b/docs/source/io_formats/nuclear_data.rst index 87b700371a..5595f2061a 100644 --- a/docs/source/io_formats/nuclear_data.rst +++ b/docs/source/io_formats/nuclear_data.rst @@ -44,6 +44,7 @@ temperature-dependent data set. For example, the data set corresponding to - **center_of_mass** (*int*) -- Whether the reference frame for scattering is center-of-mass (1) or laboratory (0) - **n_product** (*int*) -- Number of reaction products + - **redundant** (*int*) -- Whether reaction is redundant **//reactions/reaction_/K/** diff --git a/examples/jupyter/nuclear-data.ipynb b/examples/jupyter/nuclear-data.ipynb index 726de6f3e2..91923faa1d 100644 --- a/examples/jupyter/nuclear-data.ipynb +++ b/examples/jupyter/nuclear-data.ipynb @@ -621,7 +621,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "There is also `summed_reactions` attribute for cross sections (like total) which are built from summing up other cross sections." + "There is also `redundant_reactions` attribute for cross sections (like total) which are built from summing up other cross sections." ] }, { @@ -640,7 +640,7 @@ } ], "source": [ - "pprint(list(gd157.summed_reactions.values()))" + "pprint(list(gd157.redundant_reactions.values()))" ] }, { diff --git a/include/openmc/reaction.h b/include/openmc/reaction.h index 7b4d32a94d..1015a4ad27 100644 --- a/include/openmc/reaction.h +++ b/include/openmc/reaction.h @@ -34,6 +34,7 @@ public: int mt_; //!< ENDF MT value double q_value_; //!< Reaction Q value in [eV] bool scatter_in_cm_; //!< scattering system in center-of-mass? + bool redundant_; //!< redundant reaction? std::vector xs_; //!< Cross section at each temperature std::vector products_; //!< Reaction products }; @@ -48,6 +49,7 @@ extern "C" { int reaction_mt(Reaction* rx); double reaction_q_value(Reaction* rx); bool reaction_scatter_in_cm(Reaction* rx); + bool reaction_redundant(Reaction* rx); double reaction_product_decay_rate(Reaction* rx, int product); int reaction_product_emission_mode(Reaction* rx, int product); int reaction_product_particle(Reaction* rx, int product); diff --git a/openmc/data/function.py b/openmc/data/function.py index f7636f25cb..3362040dcf 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -502,7 +502,7 @@ class Sum(EqualityMixin): """Sum of multiple functions. This class allows you to create a callable object which represents the sum - of other callable objects. This is used for summed reactions whereby the + of other callable objects. This is used for redundant reactions whereby the cross section is defined as the sum of other cross sections. Parameters diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index af34c560c6..fb00bccb9d 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -86,8 +86,8 @@ class IncidentNeutron(EqualityMixin): Resonance parameters resonance_covariance : openmc.data.ResonanceCovariance or None Covariance for resonance parameters - summed_reactions : collections.OrderedDict - Contains summed cross sections, e.g., the total cross section. The keys + redundant_reactions : collections.OrderedDict + Contains redundant cross sections, e.g., the total cross section. The keys are the MT values and the values are Reaction objects. temperatures : list of str List of string representations the temperatures of the target nuclide @@ -113,18 +113,18 @@ class IncidentNeutron(EqualityMixin): self.energy = {} self._fission_energy = None self.reactions = OrderedDict() - self.summed_reactions = OrderedDict() + self.redundant_reactions = OrderedDict() self._urr = {} self._resonances = None def __contains__(self, mt): - return mt in self.reactions or mt in self.summed_reactions + return mt in self.reactions or mt in self.redundant_reactions def __getitem__(self, mt): if mt in self.reactions: return self.reactions[mt] - elif mt in self.summed_reactions: - return self.summed_reactions[mt] + elif mt in self.redundant_reactions: + return self.redundant_reactions[mt] else: raise KeyError('No reaction with MT={}.'.format(mt)) @@ -171,8 +171,8 @@ class IncidentNeutron(EqualityMixin): return self._resonance_covariance @property - def summed_reactions(self): - return self._summed_reactions + def redundant_reactions(self): + return self._redundant_reactions @property def urr(self): @@ -237,10 +237,10 @@ class IncidentNeutron(EqualityMixin): res_cov.ResonanceCovariances) self._resonance_covariance = resonance_covariance - @summed_reactions.setter - def summed_reactions(self, summed_reactions): - cv.check_type('summed reactions', summed_reactions, Mapping) - self._summed_reactions = summed_reactions + @redundant_reactions.setter + def redundant_reactions(self, redundant_reactions): + cv.check_type('redundant reactions', redundant_reactions, Mapping) + self._redundant_reactions = redundant_reactions @urr.setter def urr(self, urr): @@ -287,8 +287,8 @@ class IncidentNeutron(EqualityMixin): # Add energy grid self.energy[strT] = data.energy[strT] - # Add normal and summed reactions - for mt in chain(data.reactions, data.summed_reactions): + # Add normal and redundant reactions + for mt in chain(data.reactions, data.redundant_reactions): if mt in self: self[mt].xs[strT] = data[mt].xs[strT] else: @@ -392,7 +392,7 @@ class IncidentNeutron(EqualityMixin): self[2].xs['0K'] = Tabulated1D(x, y) def get_reaction_components(self, mt): - """Determine what reactions make up summed reaction. + """Determine what reactions make up redundant reaction. Parameters ---------- @@ -402,7 +402,7 @@ class IncidentNeutron(EqualityMixin): Returns ------- mts : list of int - ENDF MT numbers of reactions that make up the summed reaction and + ENDF MT numbers of reactions that make up the redundant reaction and have cross sections provided. """ @@ -480,6 +480,12 @@ class IncidentNeutron(EqualityMixin): tgroup = g.create_group('total_nu') rx.derived_products[0].to_hdf5(tgroup) + # Write redundant reaction data only for reactions with photon production + for rx in self.redundant_reactions.values(): + if any(p.particle == 'photon' for p in rx.products): + rx_group = rxs_group.create_group('reaction_{:03}'.format(rx.mt)) + rx.to_hdf5(rx_group) + # Write unresolved resonance probability tables if self.urr: urr_group = g.create_group('urr') @@ -555,20 +561,23 @@ class IncidentNeutron(EqualityMixin): for name, obj in sorted(rxs_group.items()): if name.startswith('reaction_'): rx = Reaction.from_hdf5(obj, data.energy) - data.reactions[rx.mt] = rx + if rx.redundant: + data.redundant_reactions[rx.mt] = rx + else: + data.reactions[rx.mt] = rx # Read total nu data if available if rx.mt in (18, 19, 20, 21, 38) and 'total_nu' in group: tgroup = group['total_nu'] rx.derived_products.append(Product.from_hdf5(tgroup)) - # Build summed reactions. Start from the highest MT number because + # Build redundant reactions. Start from the highest MT number because # high MTs never depend on lower MTs. for mt_sum in sorted(SUM_RULES, reverse=True): if mt_sum not in data: rxs = [data[mt] for mt in SUM_RULES[mt_sum] if mt in data] if len(rxs) > 0: - data.summed_reactions[mt_sum] = rx = Reaction(mt_sum) + data.redundant_reactions[mt_sum] = rx = Reaction(mt_sum) if rx.mt == 18 and 'total_nu' in group: tgroup = group['total_nu'] rx.derived_products.append(Product.from_hdf5(tgroup)) @@ -641,15 +650,17 @@ class IncidentNeutron(EqualityMixin): absorption_xs = ace.xss[ace.jxs[1] + 2 * n_energy:ace.jxs[1] + 3 * n_energy] - # Create summed reactions (total and absorption) + # Create redundant reactions (total and absorption) total = Reaction(1) total.xs[strT] = Tabulated1D(energy, total_xs) - data.summed_reactions[1] = total + total.redundant = True + data.redundant_reactions[1] = total if np.count_nonzero(absorption_xs) > 0: absorption = Reaction(27) absorption.xs[strT] = Tabulated1D(energy, absorption_xs) - data.summed_reactions[27] = absorption + absorption.redundant = True + data.redundant_reactions[27] = absorption # Read each reaction n_reaction = ace.nxs[4] + 1 @@ -671,19 +682,24 @@ class IncidentNeutron(EqualityMixin): 'cross section is given.'.format(mt)) continue - # Create summed reaction with appropriate cross section + # Create redundant reaction with appropriate cross section rx = Reaction(mt) mts = data.get_reaction_components(mt) if len(mts) == 0: warn('Photon production is present for MT={} but no ' 'reaction components exist.'.format(mt)) continue - rx.xs[strT] = Sum([data.reactions[mt_i].xs[strT] - for mt_i in mts]) - # Determine summed cross section + xss = [data.reactions[mt_i].xs[strT] for mt_i in mts] + idx = min([xs._threshold_idx if hasattr(xs, '_threshold_idx') + else 0 for xs in xss]) + rx.xs[strT] = Tabulated1D(energy[idx:], Sum(xss)(energy[idx:])) + rx.xs[strT]._threshold_idx = idx + rx.redundant = True + + # Determine redundant cross section rx.products += _get_photon_products_ace(ace, rx) - data.summed_reactions[mt] = rx + data.redundant_reactions[mt] = rx # Read unresolved resonance probability tables urr = ProbabilityTables.from_ace(ace) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 6e76d10bac..61c903158b 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -372,8 +372,8 @@ class IncidentPhoton(EqualityMixin): excitation energy), 's_collision' (collision stopping power in [eV cm\ :sup:`2`/g]), and 's_radiative' (radiative stopping power in [eV cm\ :sup:`2`/g]) - summed_reactions : collections.OrderedDict - Contains summed cross sections. The keys are MT values and the values + redundant_reactions : collections.OrderedDict + Contains redundant cross sections. The keys are MT values and the values are instances of :class:`PhotonReaction`. """ @@ -382,19 +382,19 @@ class IncidentPhoton(EqualityMixin): self.atomic_number = atomic_number self._atomic_relaxation = None self.reactions = OrderedDict() - self.summed_reactions = OrderedDict() + self.redundant_reactions = OrderedDict() self.compton_profiles = {} self.stopping_powers = {} self.bremsstrahlung = {} def __contains__(self, mt): - return mt in self.reactions or mt in self.summed_reactions + return mt in self.reactions or mt in self.redundant_reactions def __getitem__(self, mt): if mt in self.reactions: return self.reactions[mt] - elif mt in self.summed_reactions: - return self.summed_reactions[mt] + elif mt in self.redundant_reactions: + return self.redundant_reactions[mt] else: raise KeyError('No reaction with MT={}.'.format(mt)) @@ -654,7 +654,7 @@ class IncidentPhoton(EqualityMixin): return data - def export_to_hdf5(self, path, mode='a'): + def export_to_hdf5(self, path, mode='a', libver='earliest'): """Export incident photon data to an HDF5 file. Parameters @@ -667,7 +667,7 @@ class IncidentPhoton(EqualityMixin): """ # Open file and write version - f = h5py.File(path, mode, libver='latest') + f = h5py.File(path, mode, libver=libver) f.attrs['filetype'] = np.string_('data_photon') if 'version' not in f.attrs: f.attrs['version'] = np.array(HDF5_VERSION) diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index 54f0de61e7..547208dae6 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -576,14 +576,7 @@ def _get_photon_products_ace(ace, rx): # Determine corresponding reaction neutron_mt = photon_mts[i] // 1000 - # Restrict to photons that match the requested MT. Note that if the - # photon is assigned to MT=18 but the file splits fission into - # MT=19,20,21,38, we assign the photon product to each of the individual - # reactions - if neutron_mt == 18: - if rx.mt not in (18, 19, 20, 21, 38): - continue - elif neutron_mt != rx.mt: + if neutron_mt != rx.mt: continue # Create photon product and assign to reactions @@ -668,7 +661,7 @@ def _get_photon_products_endf(ev, rx): Returns ------- - photons : list of openmc.Products + products : list of openmc.Products Photons produced from reaction with given MT """ @@ -785,6 +778,8 @@ class Reaction(EqualityMixin): Indicates whether scattering kinematics should be performed in the center-of-mass or laboratory reference frame. grid above the threshold value in barns. + redundant : bool + Indicates whether or not this is a redundant reaction mt : int The ENDF MT number for this reaction. q_value : float @@ -803,6 +798,7 @@ class Reaction(EqualityMixin): def __init__(self, mt): self._center_of_mass = True + self._redundant = False self._q_value = 0. self._xs = {} self._products = [] @@ -820,6 +816,10 @@ class Reaction(EqualityMixin): def center_of_mass(self): return self._center_of_mass + @property + def redundant(self): + return self._redundant + @property def q_value(self): return self._q_value @@ -841,6 +841,11 @@ class Reaction(EqualityMixin): cv.check_type('center of mass', center_of_mass, (bool, np.bool_)) self._center_of_mass = center_of_mass + @redundant.setter + def redundant(self, redundant): + cv.check_type('redundant', redundant, (bool, np.bool_)) + self._redundant = redundant + @q_value.setter def q_value(self, q_value): cv.check_type('Q value', q_value, Real) @@ -882,6 +887,7 @@ class Reaction(EqualityMixin): group.attrs['label'] = np.string_(self.mt) group.attrs['Q_value'] = self.q_value group.attrs['center_of_mass'] = 1 if self.center_of_mass else 0 + group.attrs['redundant'] = 1 if self.redundant else 0 for T in self.xs: Tgroup = group.create_group(T) if self.xs[T] is not None: @@ -918,6 +924,7 @@ class Reaction(EqualityMixin): rx = cls(mt) rx.q_value = group.attrs['Q_value'] rx.center_of_mass = bool(group.attrs['center_of_mass']) + rx.redundant = bool(group.attrs['redundant']) # Read cross section at each temperature for T, Tgroup in group.items(): diff --git a/scripts/openmc-get-photon-data b/scripts/openmc-get-photon-data index ae268a4cc0..39e44b177f 100755 --- a/scripts/openmc-get-photon-data +++ b/scripts/openmc-get-photon-data @@ -1,4 +1,4 @@ -#!/usr/bin/env python +#!/usr/bin/env python3 """ Download ENDF/B-VII.1 ENDF data from NNDC for photo-atomic and atomic diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp index d5a47957bb..863ffd8b08 100644 --- a/src/distribution_energy.cpp +++ b/src/distribution_energy.cpp @@ -182,11 +182,13 @@ double ContinuousTabular::sample(double E) const // Interpolation for energy E1 and EK int n_energy_out = distribution_[i].e_out.size(); - double E_i_1 = distribution_[i].e_out[0]; + int n_discrete = distribution_[i].n_discrete; + double E_i_1 = distribution_[i].e_out[n_discrete]; double E_i_K = distribution_[i].e_out[n_energy_out - 1]; n_energy_out = distribution_[i+1].e_out.size(); - double E_i1_1 = distribution_[i+1].e_out[0]; + n_discrete = distribution_[i+1].n_discrete; + double E_i1_1 = distribution_[i+1].e_out[n_discrete]; double E_i1_K = distribution_[i+1].e_out[n_energy_out - 1]; double E_1 = E_i_1 + r*(E_i1_1 - E_i_1); @@ -194,25 +196,38 @@ double ContinuousTabular::sample(double E) const // Determine outgoing energy bin n_energy_out = distribution_[l].e_out.size(); + n_discrete = distribution_[l].n_discrete; double r1 = prn(); double c_k = distribution_[l].c[0]; - double c_k1; - int k; - for (k = 0; k < n_energy_out - 2; ++k) { - c_k1 = distribution_[l].c[k+1]; - if (r1 < c_k1) break; - c_k = c_k1; + int k = 0; + int end = n_energy_out - 2; + + // Discrete portion + for (int j = 0; j < n_discrete; ++j) { + k = j; + c_k = distribution_[l].c[k]; + if (r1 < c_k) { + end = j; + break; + } } - // Check to make sure 1 <= k <= NP - 1 - k = std::max(0, std::min(k, n_energy_out - 2)); + // Continuous portion + double c_k1; + for (int j = n_discrete; j < end; ++j) { + k = j; + c_k1 = distribution_[l].c[k+1]; + if (r1 < c_k1) break; + k = j + 1; + c_k = c_k1; + } double E_l_k = distribution_[l].e_out[k]; double p_l_k = distribution_[l].p[k]; double E_out; if (distribution_[l].interpolation == Interpolation::histogram) { // Histogram interpolation - if (p_l_k > 0.0) { + if (p_l_k > 0.0 && k >= n_discrete) { E_out = E_l_k + (r1 - c_k)/p_l_k; } else { E_out = E_l_k; @@ -233,7 +248,7 @@ double ContinuousTabular::sample(double E) const } // Now interpolate between incident energy bins i and i + 1 - if (!histogram_interp && n_energy_out > 1) { + if (!histogram_interp && n_energy_out > 1 && k >= n_discrete) { if (l == i) { return E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1); } else { diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 47a2d700a1..9eb50e6632 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -476,7 +476,8 @@ contains if (MTs % data(i) /= N_FISSION .and. MTs % data(i) /= N_F .and. & MTs % data(i) /= N_NF .and. MTs % data(i) /= N_2NF .and. & MTs % data(i) /= N_3NF .and. MTs % data(i) < 200 .and. & - MTs % data(i) /= N_LEVEL .and. MTs % data(i) /= ELASTIC) then + MTs % data(i) /= N_LEVEL .and. MTs % data(i) /= ELASTIC .and. & + .not. this % reactions(i) % redundant) then call index_inelastic_scatter % push_back(i) end if @@ -623,29 +624,10 @@ contains this % reaction_index(this % reactions(i) % MT) = i associate (rx => this % reactions(i)) - ! Skip total inelastic level scattering, gas production cross sections - ! (MT=200+), etc. - if (rx % MT == N_LEVEL .or. rx % MT == N_NONELASTIC) cycle - if (rx % MT > N_5N2P .and. rx % MT < N_P0) cycle - - ! Skip level cross sections if total is available - if (rx % MT >= N_P0 .and. rx % MT <= N_PC .and. find(MTs, N_P) /= -1) cycle - if (rx % MT >= N_D0 .and. rx % MT <= N_DC .and. find(MTs, N_D) /= -1) cycle - if (rx % MT >= N_T0 .and. rx % MT <= N_TC .and. find(MTs, N_T) /= -1) cycle - if (rx % MT >= N_3HE0 .and. rx % MT <= N_3HEC .and. find(MTs, N_3HE) /= -1) cycle - if (rx % MT >= N_A0 .and. rx % MT <= N_AC .and. find(MTs, N_A) /= -1) cycle - if (rx % MT >= N_2N0 .and. rx % MT <= N_2NC .and. find(MTs, N_2N) /= -1) cycle - do t = 1, n_temperature j = rx % xs_threshold(t) n = rx % xs_size(t) - ! Add contribution to total cross section - do k = j, j + n - 1 - this % xs(t) % value(XS_TOTAL,k) = this % xs(t) % & - value(XS_TOTAL,k) + rx % xs(t, k - j + 1) - end do - ! Calculate photon production cross section do k = 1, rx % products_size() if (rx % product_particle(k) == PHOTON) then @@ -658,6 +640,28 @@ contains end if end do + ! Skip total inelastic level scattering, gas production cross sections + ! (MT=200+), etc. + if (rx % MT == N_LEVEL .or. rx % MT == N_NONELASTIC) cycle + if (rx % MT > N_5N2P .and. rx % MT < N_P0) cycle + + ! Skip level cross sections if total is available + if (rx % MT >= N_P0 .and. rx % MT <= N_PC .and. find(MTs, N_P) /= -1) cycle + if (rx % MT >= N_D0 .and. rx % MT <= N_DC .and. find(MTs, N_D) /= -1) cycle + if (rx % MT >= N_T0 .and. rx % MT <= N_TC .and. find(MTs, N_T) /= -1) cycle + if (rx % MT >= N_3HE0 .and. rx % MT <= N_3HEC .and. find(MTs, N_3HE) /= -1) cycle + if (rx % MT >= N_A0 .and. rx % MT <= N_AC .and. find(MTs, N_A) /= -1) cycle + if (rx % MT >= N_2N0 .and. rx % MT <= N_2NC .and. find(MTs, N_2N) /= -1) cycle + + ! Skip redundant reactions, which are used for photon production + if (rx % redundant) cycle + + ! Add contribution to total cross section + do k = j, j + n - 1 + this % xs(t) % value(XS_TOTAL,k) = this % xs(t) % & + value(XS_TOTAL,k) + rx % xs(t, k - j + 1) + end do + ! Add contribution to absorption cross section if (is_disappearance(rx % MT)) then do k = j, j + n - 1 diff --git a/src/photon_header.F90 b/src/photon_header.F90 index bec559b995..a284c7cc02 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -163,14 +163,18 @@ contains call this % coherent_int_form_factor % from_hdf5(dset_id) call close_dataset(dset_id) - dset_id = open_dataset(rgroup, 'anomalous_real') - call this % coherent_anomalous_real % from_hdf5(dset_id) - call close_dataset(dset_id) + if (object_exists(group_id, 'anomalous_real')) then + dset_id = open_dataset(rgroup, 'anomalous_real') + call this % coherent_anomalous_real % from_hdf5(dset_id) + call close_dataset(dset_id) + end if - dset_id = open_dataset(rgroup, 'anomalous_imag') - call this % coherent_anomalous_imag % from_hdf5(dset_id) - call close_dataset(dset_id) - call close_group(rgroup) + if (object_exists(group_id, 'anomalous_imag')) then + dset_id = open_dataset(rgroup, 'anomalous_imag') + call this % coherent_anomalous_imag % from_hdf5(dset_id) + call close_dataset(dset_id) + call close_group(rgroup) + end if ! Read incoherent scattering rgroup = open_group(group_id, 'incoherent') diff --git a/src/photon_physics.F90 b/src/photon_physics.F90 index e17e0f47d0..ff631b5cf1 100644 --- a/src/photon_physics.F90 +++ b/src/photon_physics.F90 @@ -557,6 +557,8 @@ contains real(8) :: c, c_l, c_max type(BremsstrahlungData), pointer :: mat + if (p % material == MATERIAL_VOID) return + if (p % E < energy_cutoff(PHOTON)) return ! Get bremsstrahlung data for this material and particle type diff --git a/src/reaction.cpp b/src/reaction.cpp index 48eca15e08..3539fed1ae 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -14,9 +14,11 @@ Reaction::Reaction(hid_t group, const std::vector& temperatures) { read_attribute(group, "Q_value", q_value_); read_attribute(group, "mt", mt_); - int cm; - read_attribute(group, "center_of_mass", cm); - scatter_in_cm_ = (cm == 1); + int tmp; + read_attribute(group, "center_of_mass", tmp); + scatter_in_cm_ = (tmp == 1); + read_attribute(group, "redundant", tmp); + redundant_ = (tmp == 1); // Read cross section and threshold_idx data for (auto t : temperatures) { @@ -85,6 +87,8 @@ double reaction_q_value(Reaction* rx) { return rx->q_value_; } bool reaction_scatter_in_cm(Reaction* rx) { return rx->scatter_in_cm_; } +bool reaction_redundant(Reaction* rx) { return rx->redundant_; } + double reaction_product_decay_rate(Reaction* rx, int product) { return rx->products_[product - 1].decay_rate_; diff --git a/src/reaction_header.F90 b/src/reaction_header.F90 index 915bf48337..9d6add0f0a 100644 --- a/src/reaction_header.F90 +++ b/src/reaction_header.F90 @@ -20,11 +20,13 @@ module reaction_header integer(C_INT) :: MT ! ENDF MT value real(C_DOUBLE) :: Q_value ! Reaction Q value logical(C_BOOL) :: scatter_in_cm ! scattering system in center-of-mass? + logical(C_BOOL) :: redundant ! redundant reaction? contains procedure :: from_hdf5 procedure :: mt_ procedure :: q_value_ procedure :: scatter_in_cm_ + procedure :: redundant_ procedure :: product_decay_rate procedure :: product_emission_mode procedure :: product_particle @@ -64,6 +66,12 @@ module reaction_header logical(C_BOOL) :: b end function + function reaction_redundant(ptr) result(b) bind(C) + import C_PTR, C_BOOL + type(C_PTR), value :: ptr + logical(C_BOOL) :: b + end function + pure function reaction_product_decay_rate(ptr, product) result(rate) bind(C) import C_PTR, C_INT, C_DOUBLE type(C_PTR), value :: ptr @@ -159,6 +167,7 @@ contains this % MT = reaction_mt(this % ptr) this % Q_value = reaction_q_value(this % ptr) this % scatter_in_cm = reaction_scatter_in_cm(this % ptr) + this % redundant = reaction_redundant(this % ptr) end subroutine from_hdf5 function mt_(this) result(mt) @@ -182,6 +191,13 @@ contains cm = reaction_scatter_in_cm(this % ptr) end function + function redundant_(this) result(redundant) + class (Reaction), intent(in) :: this + logical(C_BOOL) :: redundant + + redundant = reaction_redundant(this % ptr) + end function + pure function product_decay_rate(this, product) result(rate) class(Reaction), intent(in) :: this integer(C_INT), intent(in) :: product diff --git a/src/secondary_correlated.cpp b/src/secondary_correlated.cpp index 1acaf6551c..e238634c3d 100644 --- a/src/secondary_correlated.cpp +++ b/src/secondary_correlated.cpp @@ -177,11 +177,13 @@ void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu) const // Interpolation for energy E1 and EK int n_energy_out = distribution_[i].e_out.size(); - double E_i_1 = distribution_[i].e_out[0]; + int n_discrete = distribution_[i].n_discrete; + double E_i_1 = distribution_[i].e_out[n_discrete]; double E_i_K = distribution_[i].e_out[n_energy_out - 1]; n_energy_out = distribution_[i+1].e_out.size(); - double E_i1_1 = distribution_[i+1].e_out[0]; + n_discrete = distribution_[i+1].n_discrete; + double E_i1_1 = distribution_[i+1].e_out[n_discrete]; double E_i1_K = distribution_[i+1].e_out[n_energy_out - 1]; double E_1 = E_i_1 + r*(E_i1_1 - E_i_1); @@ -189,24 +191,37 @@ void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu) const // Determine outgoing energy bin n_energy_out = distribution_[l].e_out.size(); + n_discrete = distribution_[l].n_discrete; double r1 = prn(); double c_k = distribution_[l].c[0]; - double c_k1; - int k; - for (k = 0; k < n_energy_out - 2; ++k) { - c_k1 = distribution_[l].c[k+1]; - if (r1 < c_k1) break; - c_k = c_k1; + int k = 0; + int end = n_energy_out - 2; + + // Discrete portion + for (int j = 0; j < n_discrete; ++j) { + k = j; + c_k = distribution_[l].c[k]; + if (r1 < c_k) { + end = j; + break; + } } - // Check to make sure 1 <= k <= NP - 1 - k = std::max(0, std::min(k, n_energy_out - 2)); + // Continuous portion + double c_k1; + for (int j = n_discrete; j < end; ++j) { + k = j; + c_k1 = distribution_[l].c[k+1]; + if (r1 < c_k1) break; + k = j + 1; + c_k = c_k1; + } double E_l_k = distribution_[l].e_out[k]; double p_l_k = distribution_[l].p[k]; if (distribution_[l].interpolation == Interpolation::histogram) { // Histogram interpolation - if (p_l_k > 0.0) { + if (p_l_k > 0.0 && k >= n_discrete) { E_out = E_l_k + (r1 - c_k)/p_l_k; } else { E_out = E_l_k; @@ -228,10 +243,12 @@ void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu) const } // Now interpolate between incident energy bins i and i + 1 - if (l == i) { - E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1); - } else { - E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1); + if (k >= n_discrete){ + if (l == i) { + E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1); + } else { + E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1); + } } // Find correlated angular distribution for closest outgoing energy bin diff --git a/src/secondary_kalbach.cpp b/src/secondary_kalbach.cpp index 3e79072100..cec5344def 100644 --- a/src/secondary_kalbach.cpp +++ b/src/secondary_kalbach.cpp @@ -144,11 +144,13 @@ void KalbachMann::sample(double E_in, double& E_out, double& mu) const // Interpolation for energy E1 and EK int n_energy_out = distribution_[i].e_out.size(); - double E_i_1 = distribution_[i].e_out[0]; + int n_discrete = distribution_[i].n_discrete; + double E_i_1 = distribution_[i].e_out[n_discrete]; double E_i_K = distribution_[i].e_out[n_energy_out - 1]; n_energy_out = distribution_[i+1].e_out.size(); - double E_i1_1 = distribution_[i+1].e_out[0]; + n_discrete = distribution_[i+1].n_discrete; + double E_i1_1 = distribution_[i+1].e_out[n_discrete]; double E_i1_K = distribution_[i+1].e_out[n_energy_out - 1]; double E_1 = E_i_1 + r*(E_i1_1 - E_i_1); @@ -156,25 +158,38 @@ void KalbachMann::sample(double E_in, double& E_out, double& mu) const // Determine outgoing energy bin n_energy_out = distribution_[l].e_out.size(); + n_discrete = distribution_[l].n_discrete; double r1 = prn(); double c_k = distribution_[l].c[0]; - double c_k1; - int k; - for (k = 0; k < n_energy_out - 2; ++k) { - c_k1 = distribution_[l].c[k+1]; - if (r1 < c_k1) break; - c_k = c_k1; + int k = 0; + int end = n_energy_out - 2; + + // Discrete portion + for (int j = 0; j < n_discrete; ++j) { + k = j; + c_k = distribution_[l].c[k]; + if (r1 < c_k) { + end = j; + break; + } } - // Check to make sure 1 <= k <= NP - 1 - k = std::max(0, std::min(k, n_energy_out - 2)); + // Continuous portion + double c_k1; + for (int j = n_discrete; j < end; ++j) { + k = j; + c_k1 = distribution_[l].c[k+1]; + if (r1 < c_k1) break; + k = j + 1; + c_k = c_k1; + } double E_l_k = distribution_[l].e_out[k]; double p_l_k = distribution_[l].p[k]; double km_r, km_a; if (distribution_[l].interpolation == Interpolation::histogram) { // Histogram interpolation - if (p_l_k > 0.0) { + if (p_l_k > 0.0 && k >= n_discrete) { E_out = E_l_k + (r1 - c_k)/p_l_k; } else { E_out = E_l_k; @@ -205,10 +220,12 @@ void KalbachMann::sample(double E_in, double& E_out, double& mu) const } // Now interpolate between incident energy bins i and i + 1 - if (l == i) { - E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1); - } else { - E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1); + if (k >= n_discrete) { + if (l == i) { + E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1); + } else { + E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1); + } } // Sampled correlated angle from Kalbach-Mann parameters diff --git a/tests/regression_tests/photon_production/inputs_true.dat b/tests/regression_tests/photon_production/inputs_true.dat index 58617b1dc9..32e3bdcf38 100644 --- a/tests/regression_tests/photon_production/inputs_true.dat +++ b/tests/regression_tests/photon_production/inputs_true.dat @@ -6,7 +6,7 @@ - + @@ -37,14 +37,14 @@ - - 14 + + 9 2 1 2 - flux + current diff --git a/tests/regression_tests/photon_production/results_true.dat b/tests/regression_tests/photon_production/results_true.dat index 84204f9b80..6f53e2d334 100644 --- a/tests/regression_tests/photon_production/results_true.dat +++ b/tests/regression_tests/photon_production/results_true.dat @@ -1,3 +1,3 @@ tally 1: -sum = 1.371553E-08 -sum_sq = 1.881158E-16 +sum = 1.550000E-02 +sum_sq = 2.402500E-04 diff --git a/tests/regression_tests/photon_production/test.py b/tests/regression_tests/photon_production/test.py index 14c321c43c..68992c867b 100644 --- a/tests/regression_tests/photon_production/test.py +++ b/tests/regression_tests/photon_production/test.py @@ -17,7 +17,7 @@ class SourceTestHarness(PyAPITestHarness): cyl = openmc.XCylinder(boundary_type='vacuum', R=1.0e-6) x_plane_left = openmc.XPlane(boundary_type='vacuum', x0=-1.0) x_plane_center = openmc.XPlane(boundary_type='transmission', x0=1.0) - x_plane_right = openmc.XPlane(boundary_type='vacuum', x0=11.0) + x_plane_right = openmc.XPlane(boundary_type='vacuum', x0=1.0e9) inner_cyl_left = openmc.Cell() inner_cyl_right = openmc.Cell() @@ -46,11 +46,11 @@ class SourceTestHarness(PyAPITestHarness): settings.source = source settings.export_to_xml() - cell_filter = openmc.CellFilter(inner_cyl_right) + surface_filter = openmc.SurfaceFilter(cyl) particle_filter = openmc.ParticleFilter('photon') tally = openmc.Tally() - tally.filters = [cell_filter, particle_filter] - tally.scores = ['flux'] + tally.filters = [surface_filter, particle_filter] + tally.scores = ['current'] tallies = openmc.Tallies([tally]) tallies.export_to_xml() @@ -65,6 +65,6 @@ class SourceTestHarness(PyAPITestHarness): return outstr -def test_source(): +def test_photon_production(): harness = SourceTestHarness('statepoint.1.h5') harness.main() diff --git a/tests/regression_tests/photon_source/test.py b/tests/regression_tests/photon_source/test.py index 30961ab5aa..61367d3756 100644 --- a/tests/regression_tests/photon_source/test.py +++ b/tests/regression_tests/photon_source/test.py @@ -56,6 +56,6 @@ class SourceTestHarness(PyAPITestHarness): return outstr -def test_source(): +def test_photon_source(): harness = SourceTestHarness('statepoint.1.h5') harness.main() diff --git a/tools/ci/travis-before-script.sh b/tools/ci/travis-before-script.sh index 6d44521a6e..f267cb0200 100755 --- a/tools/ci/travis-before-script.sh +++ b/tools/ci/travis-before-script.sh @@ -7,7 +7,7 @@ sh -e /etc/init.d/xvfb start # Download NNDC HDF5 data if [[ ! -e $HOME/nndc_hdf5/cross_sections.xml ]]; then - wget https://anl.box.com/shared/static/na85do11dfh0lb9utye2il5o6yaxx8hi.xz -O - | tar -C $HOME -xvJ + wget https://anl.box.com/shared/static/a0eflty17atnpd0pp7460exagr3nuhm7.xz -O - | tar -C $HOME -xvJ fi # Download ENDF/B-VII.1 distribution