From bad7851822be15540ad18cc0f49c7ad7edafda1f Mon Sep 17 00:00:00 2001 From: amandalund Date: Thu, 23 Aug 2018 16:29:40 -0500 Subject: [PATCH] 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()