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()