crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 13:45:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Clean up

Browse source
This commit is contained in:
amandalund 2018-08-30 11:39:12 -05:00
parent d81f1b9c7b
commit 6db62984d7
7 changed files with 72 additions and 125 deletions
Download patch file Download diff file Expand all files Collapse all files

2
include/openmc/reaction.h
View file

@ -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<TemperatureXS> xs_; //!< Cross section at each temperature
std::vector<ReactionProduct> 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);

32
openmc/data/neutron.py
View file

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

14
openmc/data/reaction.py
View file

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

80
src/nuclide_header.F90
View file

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

43
src/physics.F90
View file

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

10
src/reaction.cpp
View file

@ -14,9 +14,11 @@ Reaction::Reaction(hid_t group, const std::vector<int>& 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_;

16
src/reaction_header.F90
View file

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