crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Account for photon production of summed reactions

Browse source
This commit is contained in:
amandalund 2018-08-23 16:29:40 -05:00
parent fa9a6610af
commit bad7851822
10 changed files with 124 additions and 20 deletions
Download patch file Download diff file Expand all files Collapse all files

24
openmc/data/neutron.py
View file

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

2
openmc/data/reaction.py
View file

@ -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
"""

2
src/distribution_energy.cpp
View file

@ -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;

49
src/nuclide_header.F90
View file

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

43
src/physics.F90
View file

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

2
src/secondary_correlated.cpp
View file

@ -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;

2
src/secondary_kalbach.cpp
View file

@ -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;

8
tests/regression_tests/photon_production/inputs_true.dat
View file

@ -6,7 +6,7 @@
<surface boundary="vacuum" coeffs="0.0 0.0 1e-06" id="9" type="x-cylinder" />
<surface boundary="vacuum" coeffs="-1.0" id="10" type="x-plane" />
<surface coeffs="1.0" id="11" type="x-plane" />
<surface boundary="vacuum" coeffs="11.0" id="12" type="x-plane" />
<surface boundary="vacuum" coeffs="1000000000.0" id="12" type="x-plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
@ -37,14 +37,14 @@
</settings>
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<filter id="1" type="cell">
<bins>14</bins>
<filter id="1" type="surface">
<bins>9</bins>
</filter>
<filter id="2" type="particle">
<bins>2</bins>
</filter>
<tally id="1">
<filters>1 2</filters>
<scores>flux</scores>
<scores>current</scores>
</tally>
</tallies>

4
tests/regression_tests/photon_production/results_true.dat
View file

@ -1,3 +1,3 @@
tally 1:
sum = 1.371553E-08
sum_sq = 1.881158E-16
sum = 1.550000E-02
sum_sq = 2.402500E-04

8
tests/regression_tests/photon_production/test.py
View file

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