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

Merge pull request #1068 from amandalund/photon-production-fix

Browse source 
Photon production fix
This commit is contained in:
Paul Romano 2018-09-13 20:14:15 -05:00 • committed by GitHub
commit 7a907f60f9
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
21 changed files with 238 additions and 133 deletions
Download patch file Download diff file Expand all files Collapse all files

1
docs/source/io_formats/nuclear_data.rst
View file

@ -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
**/<nuclide name>/reactions/reaction_<mt>/<TTT>K/**

4
examples/jupyter/nuclear-data.ipynb
View file

@ -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()))"
]
},
{

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 redundant_; //!< redundant 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_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);

2
openmc/data/function.py
View file

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

70
openmc/data/neutron.py
View file

@ -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,6 +480,12 @@ class IncidentNeutron(EqualityMixin):
tgroup = g.create_group('total_nu')
rx.derived_products[0].to_hdf5(tgroup)
# 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)
# Write unresolved resonance probability tables
if self.urr:
urr_group = g.create_group('urr')
@ -555,20 +561,23 @@ 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.redundant:
data.redundant_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))
# 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))
@ -641,15 +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)
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)
data.summed_reactions[27] = absorption
absorption.redundant = True
data.redundant_reactions[27] = absorption
# Read each reaction
n_reaction = ace.nxs[4] + 1
@ -671,19 +682,24 @@ class IncidentNeutron(EqualityMixin):
'cross section is given.'.format(mt))
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:
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])
# Determine summed cross section
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
rx.redundant = True
# 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)

16
openmc/data/photon.py
View file

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

25
openmc/data/reaction.py
View file

@ -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
@ -668,7 +661,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
"""
@ -785,6 +778,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.
redundant : bool
Indicates whether or not this is a redundant reaction
mt : int
The ENDF MT number for this reaction.
q_value : float
@ -803,6 +798,7 @@ class Reaction(EqualityMixin):
def __init__(self, mt):
self._center_of_mass = True
self._redundant = False
self._q_value = 0.
self._xs = {}
self._products = []
@ -820,6 +816,10 @@ class Reaction(EqualityMixin):
def center_of_mass(self):
return self._center_of_mass
@property
def redundant(self):
return self._redundant
@property
def q_value(self):
return self._q_value
@ -841,6 +841,11 @@ class Reaction(EqualityMixin):
cv.check_type('center of mass', center_of_mass, (bool, np.bool_))
self._center_of_mass = center_of_mass
@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):
cv.check_type('Q value', q_value, Real)
@ -882,6 +887,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['redundant'] = 1 if self.redundant else 0
for T in self.xs:
Tgroup = group.create_group(T)
if self.xs[T] is not None:
@ -918,6 +924,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.redundant = bool(group.attrs['redundant'])
# Read cross section at each temperature
for T, Tgroup in group.items():

2
scripts/openmc-get-photon-data
View file

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

39
src/distribution_energy.cpp
View file

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

44
src/nuclide_header.F90
View file

@ -476,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) % redundant) then
call index_inelastic_scatter % push_back(i)
end if
@ -623,29 +624,10 @@ contains
this % reaction_index(this % reactions(i) % MT) = i
associate (rx => this % reactions(i))
! 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
do t = 1, n_temperature
j = rx % xs_threshold(t)
n = rx % xs_size(t)
! Add contribution to total cross section
do k = j, j + n - 1
this % xs(t) % value(XS_TOTAL,k) = this % xs(t) % &
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
@ -658,6 +640,28 @@ contains
end if
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 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 redundant reactions, which are used for photon production
if (rx % redundant) cycle
! Add contribution to total cross section
do k = j, j + n - 1
this % xs(t) % value(XS_TOTAL,k) = this % xs(t) % &
value(XS_TOTAL,k) + rx % xs(t, k - j + 1)
end do
! Add contribution to absorption cross section
if (is_disappearance(rx % MT)) then
do k = j, j + n - 1

18
src/photon_header.F90
View file

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

2
src/photon_physics.F90
View file

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

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, "redundant", tmp);
redundant_ = (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_redundant(Reaction* rx) { return rx->redundant_; }
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) :: redundant ! redundant reaction?
contains
procedure :: from_hdf5
procedure :: mt_
procedure :: q_value_
procedure :: scatter_in_cm_
procedure :: redundant_
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_redundant(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 % redundant = reaction_redundant(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 redundant_(this) result(redundant)
class (Reaction), intent(in) :: this
logical(C_BOOL) :: redundant
redundant = reaction_redundant(this % ptr)
end function
pure function product_decay_rate(this, product) result(rate)
class(Reaction), intent(in) :: this
integer(C_INT), intent(in) :: product

47
src/secondary_correlated.cpp
View file

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

47
src/secondary_kalbach.cpp
View file

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

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

10
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()
@ -65,6 +65,6 @@ class SourceTestHarness(PyAPITestHarness):
return outstr
def test_source():
def test_photon_production():
harness = SourceTestHarness('statepoint.1.h5')
harness.main()

2
tests/regression_tests/photon_source/test.py
View file

@ -56,6 +56,6 @@ class SourceTestHarness(PyAPITestHarness):
return outstr
def test_source():
def test_photon_source():
harness = SourceTestHarness('statepoint.1.h5')
harness.main()

2
tools/ci/travis-before-script.sh
View file

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