diff --git a/include/openmc/reaction.h b/include/openmc/reaction.h index 7b4d32a94..9e8475227 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 summed_; //!< summed 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_summed(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/neutron.py b/openmc/data/neutron.py index f0484221d..b209d28ae 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -480,6 +480,12 @@ class IncidentNeutron(EqualityMixin): tgroup = g.create_group('total_nu') rx.derived_products[0].to_hdf5(tgroup) + # Write summed reaction data only for reactions with photon production + for rx in self.summed_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') @@ -492,14 +498,6 @@ class IncidentNeutron(EqualityMixin): fer_group = g.create_group('fission_energy_release') self.fission_energy.to_hdf5(fer_group) - # Write summed reaction data only for reactions with photon production - for rx in self.summed_reactions.values(): - if any([p.particle == 'photon' for p in rx.products]): - if not 'summed_reactions' in g: - srxs_group = g.create_group('summed_reactions') - rx_group = srxs_group.create_group('reaction_{:03}'.format(rx.mt)) - rx.to_hdf5(rx_group) - f.close() @classmethod @@ -563,22 +561,17 @@ 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.summed: + data.summed_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)) - # Read summed reaction data - if 'summed_reactions' in group: - srxs_group = group['summed_reactions'] - for name, obj in sorted(srxs_group.items()): - if name.startswith('reaction_'): - rx = Reaction.from_hdf5(obj, data.energy) - data.summed_reactions[rx.mt] = rx - - # Build summed reactions. Start from the highest MT number because + # Build summed 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: @@ -660,11 +653,13 @@ class IncidentNeutron(EqualityMixin): # Create summed reactions (total and absorption) total = Reaction(1) total.xs[strT] = Tabulated1D(energy, total_xs) + total.summed = True data.summed_reactions[1] = total if np.count_nonzero(absorption_xs) > 0: absorption = Reaction(27) absorption.xs[strT] = Tabulated1D(energy, absorption_xs) + absorption.summed = True data.summed_reactions[27] = absorption # Read each reaction @@ -700,6 +695,7 @@ class IncidentNeutron(EqualityMixin): else 0 for xs in xss]) rx.xs[strT] = Tabulated1D(energy[idx:], Sum(xss)(energy[idx:])) rx.xs[strT]._threshold_idx = idx + rx.summed = True # Determine summed cross section rx.products += _get_photon_products_ace(ace, rx) diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index 1f0028fa5..2bffdc43a 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -785,6 +785,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. + summed : bool + Indicates whether or not this is a summed reactions mt : int The ENDF MT number for this reaction. q_value : float @@ -803,6 +805,7 @@ class Reaction(EqualityMixin): def __init__(self, mt): self._center_of_mass = True + self._summed = False self._q_value = 0. self._xs = {} self._products = [] @@ -820,6 +823,10 @@ class Reaction(EqualityMixin): def center_of_mass(self): return self._center_of_mass + @property + def summed(self): + return self._summed + @property def q_value(self): return self._q_value @@ -841,6 +848,11 @@ class Reaction(EqualityMixin): cv.check_type('center of mass', center_of_mass, (bool, np.bool_)) self._center_of_mass = center_of_mass + @summed.setter + def summed(self, summed): + cv.check_type('summed', summed, (bool, np.bool_)) + self._summed = summed + @q_value.setter def q_value(self, q_value): cv.check_type('Q value', q_value, Real) @@ -882,6 +894,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['summed'] = 1 if self.summed else 0 for T in self.xs: Tgroup = group.create_group(T) if self.xs[T] is not None: @@ -918,6 +931,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.summed = bool(group.attrs['summed']) # Read cross section at each temperature for T, Tgroup in group.items(): diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 5ed627abc..6330cce43 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -94,7 +94,6 @@ module nuclide_header ! Reactions type(Reaction), allocatable :: reactions(:) - type(Reaction), allocatable :: summed_reactions(:) integer, allocatable :: index_inelastic_scatter(:) ! Array that maps MT values to index in reactions; used at tally-time. Note @@ -477,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) % summed) then call index_inelastic_scatter % push_back(i) end if @@ -491,30 +491,6 @@ contains this % index_inelastic_scatter = & index_inelastic_scatter % data(1: index_inelastic_scatter % size()) - ! Read summed reactions if present. These are needed for photon production. - if (object_exists(group_id, 'summed_reactions')) then - ! Get MT values for summed reactions based on group names - call MTs % clear() - rxs_group = open_group(group_id, 'summed_reactions') - call get_groups(rxs_group, grp_names) - do j = 1, size(grp_names) - if (starts_with(grp_names(j), "reaction_")) then - call MTs % push_back(int(str_to_int(grp_names(j)(10:12)))) - end if - end do - - ! Read summed reactions - allocate(this % summed_reactions(MTs % size())) - do i = 1, size(this % summed_reactions) - rx_group = open_group(rxs_group, 'reaction_' // trim(& - zero_padded(MTs % data(i), 3))) - - call this % summed_reactions(i) % from_hdf5(rx_group, temps_to_read) - call close_group(rx_group) - end do - call close_group(rxs_group) - end if - ! Read unresolved resonance probability tables if present if (object_exists(group_id, 'urr')) then this % urr_present = .true. @@ -664,24 +640,22 @@ contains end do end if end do - 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 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 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) + ! Skip summed reactions, which are used for photon production + if (rx % summed) cycle ! Add contribution to total cross section do k = j, j + n - 1 @@ -730,30 +704,6 @@ contains end associate ! rx end do ! reactions - ! Add contribution to photon production cross section from summed reactions - if (allocated(this % summed_reactions)) then - do i = 1, size(this % summed_reactions) - associate (rx => this % summed_reactions(i)) - do t = 1, n_temperature - j = rx % xs_threshold(t) - n = rx % xs_size(t) - - ! Calculate photon production cross section - do k = 1, rx % products_size() - if (rx % product_particle(k) == PHOTON) then - do l = 1, n - this % xs(t) % value(XS_PHOTON_PROD,l+j-1) = & - this % xs(t) % value(XS_PHOTON_PROD,l+j-1) + & - rx % xs(t, l) * rx % product_yield(k, & - this % grid(t) % energy(l+j-1)) - end do - end if - end do - end do - end associate - end do - end if - ! Determine number of delayed neutron precursors if (this % fissionable) then associate (rx => this % reactions(this % index_fission(1))) diff --git a/src/physics.F90 b/src/physics.F90 index 3f24d0a50..992e48767 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -565,12 +565,11 @@ contains ! SAMPLE_PHOTON_PRODUCT !=============================================================================== - subroutine sample_photon_product(i_nuclide, E, i_reaction, i_product, summed) + subroutine sample_photon_product(i_nuclide, E, i_reaction, i_product) integer, intent(in) :: i_nuclide ! index in nuclides array real(8), intent(in) :: E ! energy of neutron integer, intent(out) :: i_reaction ! index in nuc % reactions array integer, intent(out) :: i_product ! index in reaction % products array - logical, intent(out) :: summed ! whether a summed reaction was sampled integer :: i_grid integer :: i_temp @@ -591,7 +590,6 @@ contains f = micro_xs(i_nuclide) % interp_factor cutoff = prn() * micro_xs(i_nuclide) % photon_prod prob = ZERO - summed = .false. ! Loop through each reaction type REACTION_LOOP: do i_reaction = 1, size(nuc % reactions) @@ -615,33 +613,6 @@ contains end do end associate end do REACTION_LOOP - - ! Loop through each summed reaction type - if (allocated(nuc % summed_reactions)) then - SUMMED_REACTION_LOOP: do i_reaction = 1, size(nuc % summed_reactions) - associate (rx => nuc % summed_reactions(i_reaction)) - threshold = rx % xs_threshold(i_temp) - - ! if energy is below threshold for this reaction, skip it - if (i_grid < threshold) cycle - - do i_product = 1, rx % products_size() - if (rx % product_particle(i_product) == PHOTON) then - summed = .true. - - ! add to cumulative probability - yield = rx % product_yield(i_product, E) - prob = prob + ((ONE - f) * rx % xs(i_temp, i_grid - threshold + 1) & - + f*(rx % xs(i_temp, i_grid - threshold + 2))) * yield - - if (prob > cutoff) return - last_valid_reaction = i_reaction - last_valid_product = i_product - end if - end do - end associate - end do SUMMED_REACTION_LOOP - end if end associate i_reaction = last_valid_reaction @@ -1500,7 +1471,6 @@ contains real(8) :: uvw(3) integer :: nu integer :: i - logical :: summed ! Sample the number of photons produced nu_t = p % wgt * micro_xs(i_nuclide) % photon_prod / & @@ -1515,16 +1485,11 @@ contains do i = 1, nu ! Sample the reaction and product - call sample_photon_product(i_nuclide, p % E, i_reaction, i_product, summed) + call sample_photon_product(i_nuclide, p % E, i_reaction, i_product) ! Sample the outgoing energy and angle - if (summed) then - call nuclides(i_nuclide) % summed_reactions(i_reaction) % & - product_sample(i_product, p % E, E, mu) - else - call nuclides(i_nuclide) % reactions(i_reaction) % & - product_sample(i_product, p % E, E, mu) - end if + call nuclides(i_nuclide) % reactions(i_reaction) % & + product_sample(i_product, p % E, E, mu) ! Sample the new direction uvw = rotate_angle(p % coord(1) % uvw, mu) diff --git a/src/reaction.cpp b/src/reaction.cpp index 48eca15e0..a625b9686 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, "summed", tmp); + summed_ = (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_summed(Reaction* rx) { return rx->summed_; } + 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 915bf4833..f58fa1035 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) :: summed ! summed reaction? contains procedure :: from_hdf5 procedure :: mt_ procedure :: q_value_ procedure :: scatter_in_cm_ + procedure :: summed_ 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_summed(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 % summed = reaction_summed(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 summed_(this) result(summed) + class (Reaction), intent(in) :: this + logical(C_BOOL) :: summed + + summed = reaction_summed(this % ptr) + end function + pure function product_decay_rate(this, product) result(rate) class(Reaction), intent(in) :: this integer(C_INT), intent(in) :: product