From bad7851822be15540ad18cc0f49c7ad7edafda1f Mon Sep 17 00:00:00 2001 From: amandalund Date: Thu, 23 Aug 2018 16:29:40 -0500 Subject: [PATCH 01/10] Account for photon production of summed reactions --- openmc/data/neutron.py | 24 ++++++++- openmc/data/reaction.py | 2 +- src/distribution_energy.cpp | 2 +- src/nuclide_header.F90 | 49 +++++++++++++++++++ src/physics.F90 | 43 ++++++++++++++-- src/secondary_correlated.cpp | 2 +- src/secondary_kalbach.cpp | 2 +- .../photon_production/inputs_true.dat | 8 +-- .../photon_production/results_true.dat | 4 +- .../photon_production/test.py | 8 +-- 10 files changed, 124 insertions(+), 20 deletions(-) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index af34c560c6..f0484221d6 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -492,6 +492,14 @@ 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 @@ -562,6 +570,14 @@ class IncidentNeutron(EqualityMixin): 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 # high MTs never depend on lower MTs. for mt_sum in sorted(SUM_RULES, reverse=True): @@ -678,8 +694,12 @@ class IncidentNeutron(EqualityMixin): 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]) + + 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 # 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 54f0de61e7..1f0028fa5c 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -668,7 +668,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 """ diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp index cadeeb3df1..5406c8136a 100644 --- a/src/distribution_energy.cpp +++ b/src/distribution_energy.cpp @@ -211,7 +211,7 @@ double ContinuousTabular::sample(double E) const double E_out; if (distribution_[l].interpolation == Interpolation::histogram) { // Histogram interpolation - if (p_l_k > 0.0) { + if (p_l_k > 0.0 && k >= distribution_[l].n_discrete) { E_out = E_l_k + (r1 - c_k)/p_l_k; } else { E_out = E_l_k; diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 8cd225bc98..8be9a9e4e9 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -96,6 +96,7 @@ 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 @@ -492,6 +493,30 @@ 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. @@ -702,6 +727,30 @@ 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 d5f6302a82..5f70709058 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -565,11 +565,12 @@ contains ! SAMPLE_PHOTON_PRODUCT !=============================================================================== - subroutine sample_photon_product(i_nuclide, E, i_reaction, i_product) + subroutine sample_photon_product(i_nuclide, E, i_reaction, i_product, summed) 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 @@ -590,6 +591,7 @@ 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) @@ -613,6 +615,33 @@ 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 @@ -1718,6 +1747,7 @@ 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 / & @@ -1732,11 +1762,16 @@ contains do i = 1, nu ! Sample the reaction and product - call sample_photon_product(i_nuclide, p % E, i_reaction, i_product) + call sample_photon_product(i_nuclide, p % E, i_reaction, i_product, summed) ! Sample the outgoing energy and angle - call nuclides(i_nuclide) % reactions(i_reaction) % & - product_sample(i_product, p % E, E, mu) + 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 ! Sample the new direction uvw = rotate_angle(p % coord(1) % uvw, mu) diff --git a/src/secondary_correlated.cpp b/src/secondary_correlated.cpp index b1653f3170..9156830bb0 100644 --- a/src/secondary_correlated.cpp +++ b/src/secondary_correlated.cpp @@ -205,7 +205,7 @@ void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu) const 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 >= distribution_[l].n_discrete) { E_out = E_l_k + (r1 - c_k)/p_l_k; } else { E_out = E_l_k; diff --git a/src/secondary_kalbach.cpp b/src/secondary_kalbach.cpp index e8941ff93d..e41c7e9aa4 100644 --- a/src/secondary_kalbach.cpp +++ b/src/secondary_kalbach.cpp @@ -173,7 +173,7 @@ void KalbachMann::sample(double E_in, double& E_out, double& mu) const 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 >= distribution_[l].n_discrete) { E_out = E_l_k + (r1 - c_k)/p_l_k; } else { E_out = E_l_k; 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..150cd388a9 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() From c4311bb5d22f21130aafcb5d2bdc249d199b203b Mon Sep 17 00:00:00 2001 From: amandalund Date: Wed, 29 Aug 2018 13:19:58 -0500 Subject: [PATCH 02/10] Fix calculation of photon production cross section to not skip over some reactions with photon production --- src/nuclide_header.F90 | 29 +++++++++++++++++------------ 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index f7b6589f9b..5ed627abc3 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -649,6 +649,23 @@ contains this % reaction_index(this % reactions(i) % MT) = i associate (rx => this % 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 + ! Skip total inelastic level scattering, gas production cross sections ! (MT=200+), etc. if (rx % MT == N_LEVEL .or. rx % MT == N_NONELASTIC) cycle @@ -672,18 +689,6 @@ contains 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 - 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 - ! Add contribution to absorption cross section if (is_disappearance(rx % MT)) then do k = j, j + n - 1 From 6db62984d727095244befada0974571e9bdb5d72 Mon Sep 17 00:00:00 2001 From: amandalund Date: Thu, 30 Aug 2018 11:39:12 -0500 Subject: [PATCH 03/10] Clean up --- include/openmc/reaction.h | 2 + openmc/data/neutron.py | 32 +++++++--------- openmc/data/reaction.py | 14 +++++++ src/nuclide_header.F90 | 80 ++++++++------------------------------- src/physics.F90 | 43 ++------------------- src/reaction.cpp | 10 +++-- src/reaction_header.F90 | 16 ++++++++ 7 files changed, 72 insertions(+), 125 deletions(-) diff --git a/include/openmc/reaction.h b/include/openmc/reaction.h index 7b4d32a94d..9e84752275 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 f0484221d6..b209d28aec 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 1f0028fa5c..2bffdc43ac 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 5ed627abc3..6330cce43e 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 3f24d0a50c..992e487675 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 48eca15e08..a625b9686c 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 915bf48337..f58fa1035d 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 From 3f0af69aef521b2e92323d6bf87c1335b56c2ffd Mon Sep 17 00:00:00 2001 From: amandalund Date: Tue, 11 Sep 2018 08:46:12 -0500 Subject: [PATCH 04/10] Fix energy distributions to correctly sample discrete lines --- src/distribution_energy.cpp | 39 +++++++++++++++++++++--------- src/photon_header.F90 | 18 ++++++++------ src/secondary_correlated.cpp | 47 ++++++++++++++++++++++++------------ src/secondary_kalbach.cpp | 47 ++++++++++++++++++++++++------------ 4 files changed, 102 insertions(+), 49 deletions(-) diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp index d0d0feb26b..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 && k >= distribution_[l].n_discrete) { + 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/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/secondary_correlated.cpp b/src/secondary_correlated.cpp index ec741a1e60..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 && k >= distribution_[l].n_discrete) { + 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 0779230688..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 && k >= distribution_[l].n_discrete) { + 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 From 9d335ec91877515f1703114238383811e43f6f7f Mon Sep 17 00:00:00 2001 From: amandalund Date: Wed, 12 Sep 2018 00:04:24 -0500 Subject: [PATCH 05/10] Don't duplicate fission cross sections in photon production summed reactions; update nndc data link --- openmc/data/neutron.py | 4 ++++ src/photon_physics.F90 | 2 ++ tools/ci/travis-before-script.sh | 2 +- 3 files changed, 7 insertions(+), 1 deletion(-) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index b209d28aec..2b41292501 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -682,6 +682,10 @@ class IncidentNeutron(EqualityMixin): 'cross section is given.'.format(mt)) continue + # Don't duplicate fission cross sections + if mt == 18: + continue + # Create summed reaction with appropriate cross section rx = Reaction(mt) mts = data.get_reaction_components(mt) 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/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 From dfd7bb9437836ae9bbf783f948ea0545642b6ecb Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 13 Sep 2018 07:18:47 -0500 Subject: [PATCH 06/10] Rename summed -> redundant to make ENDF language --- docs/source/io_formats/nuclear_data.rst | 1 + examples/jupyter/nuclear-data.ipynb | 4 +- include/openmc/reaction.h | 4 +- openmc/data/function.py | 2 +- openmc/data/neutron.py | 62 ++++++++++++------------- openmc/data/photon.py | 12 ++--- openmc/data/reaction.py | 22 ++++----- src/nuclide_header.F90 | 6 +-- src/reaction.cpp | 6 +-- src/reaction_header.F90 | 14 +++--- 10 files changed, 67 insertions(+), 66 deletions(-) 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 9e84752275..1015a4ad27 100644 --- a/include/openmc/reaction.h +++ b/include/openmc/reaction.h @@ -34,7 +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? + bool redundant_; //!< redundant reaction? std::vector xs_; //!< Cross section at each temperature std::vector products_; //!< Reaction products }; @@ -49,7 +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); + 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 2b41292501..6d7218d40e 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,8 +480,8 @@ 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(): + # 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) @@ -561,8 +561,8 @@ class IncidentNeutron(EqualityMixin): for name, obj in sorted(rxs_group.items()): if name.startswith('reaction_'): rx = Reaction.from_hdf5(obj, data.energy) - if rx.summed: - data.summed_reactions[rx.mt] = rx + if rx.redundant: + data.redundant_reactions[rx.mt] = rx else: data.reactions[rx.mt] = rx @@ -571,13 +571,13 @@ class IncidentNeutron(EqualityMixin): 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)) @@ -650,17 +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) - total.summed = True - 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) - absorption.summed = True - data.summed_reactions[27] = absorption + absorption.redundant = True + data.redundant_reactions[27] = absorption # Read each reaction n_reaction = ace.nxs[4] + 1 @@ -686,7 +686,7 @@ class IncidentNeutron(EqualityMixin): if mt == 18: 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: @@ -699,11 +699,11 @@ 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 + rx.redundant = True - # Determine summed cross section + # 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..3de325d329 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)) diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index 2bffdc43ac..26b51363f6 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -785,8 +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 + redundant : bool + Indicates whether or not this is a redundant reaction mt : int The ENDF MT number for this reaction. q_value : float @@ -805,7 +805,7 @@ class Reaction(EqualityMixin): def __init__(self, mt): self._center_of_mass = True - self._summed = False + self._redundant = False self._q_value = 0. self._xs = {} self._products = [] @@ -824,8 +824,8 @@ class Reaction(EqualityMixin): return self._center_of_mass @property - def summed(self): - return self._summed + def redundant(self): + return self._redundant @property def q_value(self): @@ -848,10 +848,10 @@ 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 + @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): @@ -894,7 +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 + 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: @@ -931,7 +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']) + rx.redundant = bool(group.attrs['redundant']) # 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 8d1db8f2e5..9eb50e6632 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -477,7 +477,7 @@ contains 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 .and. & - .not. this % reactions(i) % summed) then + .not. this % reactions(i) % redundant) then call index_inelastic_scatter % push_back(i) end if @@ -653,8 +653,8 @@ contains 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 summed reactions, which are used for photon production - if (rx % summed) 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 diff --git a/src/reaction.cpp b/src/reaction.cpp index a625b9686c..3539fed1ae 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -17,8 +17,8 @@ Reaction::Reaction(hid_t group, const std::vector& temperatures) int tmp; read_attribute(group, "center_of_mass", tmp); scatter_in_cm_ = (tmp == 1); - read_attribute(group, "summed", tmp); - summed_ = (tmp == 1); + read_attribute(group, "redundant", tmp); + redundant_ = (tmp == 1); // Read cross section and threshold_idx data for (auto t : temperatures) { @@ -87,7 +87,7 @@ 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_; } +bool reaction_redundant(Reaction* rx) { return rx->redundant_; } double reaction_product_decay_rate(Reaction* rx, int product) { diff --git a/src/reaction_header.F90 b/src/reaction_header.F90 index f58fa1035d..9d6add0f0a 100644 --- a/src/reaction_header.F90 +++ b/src/reaction_header.F90 @@ -20,13 +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? + logical(C_BOOL) :: redundant ! redundant reaction? contains procedure :: from_hdf5 procedure :: mt_ procedure :: q_value_ procedure :: scatter_in_cm_ - procedure :: summed_ + procedure :: redundant_ procedure :: product_decay_rate procedure :: product_emission_mode procedure :: product_particle @@ -66,7 +66,7 @@ module reaction_header logical(C_BOOL) :: b end function - function reaction_summed(ptr) result(b) bind(C) + function reaction_redundant(ptr) result(b) bind(C) import C_PTR, C_BOOL type(C_PTR), value :: ptr logical(C_BOOL) :: b @@ -167,7 +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) + this % redundant = reaction_redundant(this % ptr) end subroutine from_hdf5 function mt_(this) result(mt) @@ -191,11 +191,11 @@ contains cm = reaction_scatter_in_cm(this % ptr) end function - function summed_(this) result(summed) + function redundant_(this) result(redundant) class (Reaction), intent(in) :: this - logical(C_BOOL) :: summed + logical(C_BOOL) :: redundant - summed = reaction_summed(this % ptr) + redundant = reaction_redundant(this % ptr) end function pure function product_decay_rate(this, product) result(rate) From 973016c2d8904fd1b0144463e091572119a5dd03 Mon Sep 17 00:00:00 2001 From: amandalund Date: Thu, 13 Sep 2018 11:27:40 -0500 Subject: [PATCH 07/10] Address @paulromano comments in #1068 --- openmc/data/neutron.py | 2 +- tests/regression_tests/photon_production/test.py | 2 +- tests/regression_tests/photon_source/test.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index 6d7218d40e..842dfe7351 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -482,7 +482,7 @@ class IncidentNeutron(EqualityMixin): # 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]): + 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) diff --git a/tests/regression_tests/photon_production/test.py b/tests/regression_tests/photon_production/test.py index 150cd388a9..68992c867b 100644 --- a/tests/regression_tests/photon_production/test.py +++ b/tests/regression_tests/photon_production/test.py @@ -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() From 6280cb45c303eee57b969f8982da133312514d02 Mon Sep 17 00:00:00 2001 From: amandalund Date: Thu, 13 Sep 2018 11:40:05 -0500 Subject: [PATCH 08/10] Fix photon HDF5 library to be compatible with earlier version --- openmc/data/photon.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 3de325d329..5ba5f105a2 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -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 From d02d565aafe66419fc3d220cf71d0ce87c8706c4 Mon Sep 17 00:00:00 2001 From: amandalund Date: Thu, 13 Sep 2018 11:44:47 -0500 Subject: [PATCH 09/10] One more fix for HDF5 version compatibility --- openmc/data/photon.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 5ba5f105a2..61c903158b 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -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) From 3ea5d09e12fc7606790be89333f1cfaef4e67b75 Mon Sep 17 00:00:00 2001 From: amandalund Date: Thu, 13 Sep 2018 16:14:13 -0500 Subject: [PATCH 10/10] Don't assign MT=18 photon production to partial fission MTs --- openmc/data/neutron.py | 4 ---- openmc/data/reaction.py | 9 +-------- scripts/openmc-get-photon-data | 2 +- 3 files changed, 2 insertions(+), 13 deletions(-) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index 842dfe7351..fb00bccb9d 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -682,10 +682,6 @@ class IncidentNeutron(EqualityMixin): 'cross section is given.'.format(mt)) continue - # Don't duplicate fission cross sections - if mt == 18: - continue - # Create redundant reaction with appropriate cross section rx = Reaction(mt) mts = data.get_reaction_components(mt) diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index 26b51363f6..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 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