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

Merge branch 'develop' into particle_owned_seeds

Browse source
This commit is contained in:
John Tramm 2019-12-05 12:05:09 -06:00 • committed by GitHub
commit 08ce3f4614
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
4 changed files with 128 additions and 74 deletions
Download patch file Download diff file Expand all files Collapse all files

10
CMakeLists.txt
View file

@ -152,6 +152,9 @@ target_compile_definitions(gsl INTERFACE GSL_THROW_ON_CONTRACT_VIOLATION)
# RPATH information
#===============================================================================
# Provide install directory variables as defined by GNU coding standards
include(GNUInstallDirs)
# This block of code ensures that dynamic libraries can be found via the RPATH
# whether the executable is the original one from the build directory or the
# installed one in CMAKE_INSTALL_PREFIX. Ref:
@ -164,16 +167,14 @@ set(CMAKE_SKIP_BUILD_RPATH FALSE)
# (but later on when installing)
set(CMAKE_BUILD_WITH_INSTALL_RPATH FALSE)
set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib")
# add the automatically determined parts of the RPATH
# which point to directories outside the build tree to the install RPATH
set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)
# the RPATH to be used when installing, but only if it's not a system directory
list(FIND CMAKE_PLATFORM_IMPLICIT_LINK_DIRECTORIES "${CMAKE_INSTALL_PREFIX}/lib" isSystemDir)
list(FIND CMAKE_PLATFORM_IMPLICIT_LINK_DIRECTORIES "${CMAKE_INSTALL_FULL_LIBDIR}" isSystemDir)
if("${isSystemDir}" STREQUAL "-1")
set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib")
set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_FULL_LIBDIR}")
endif()
#===============================================================================
@ -360,7 +361,6 @@ add_custom_command(TARGET libopenmc POST_BUILD
# Install executable, scripts, manpage, license
#===============================================================================
include(GNUInstallDirs)
set(INSTALL_CONFIGDIR ${CMAKE_INSTALL_LIBDIR}/cmake/OpenMC)
install(TARGETS openmc libopenmc pugixml faddeeva gsl
EXPORT openmc-targets

94
docs/source/methods/neutron_physics.rst
View file

@ -358,9 +358,9 @@ secondary energy and angle sampling.
For a reaction with secondary products, it is necessary to determine the
outgoing angle and energy of the products. For any reaction other than elastic
and level inelastic scattering, the outgoing energy must be determined based on
tabulated or parameterized data. The `ENDF-6 Format`_ specifies a variety of
ways that the secondary energy distribution can be represented. ENDF File 5
contains uncorrelated energy distribution whereas ENDF File 6 contains
tabulated or parameterized data. The `ENDF-6 Format <endf102>`_ specifies a
variety of ways that the secondary energy distribution can be represented. ENDF
File 5 contains uncorrelated energy distribution whereas ENDF File 6 contains
correlated energy-angle distributions. The ACE format specifies its own
representations based loosely on the formats given in ENDF-6. OpenMC's HDF5
nuclear data files use a combination of ENDF and ACE distributions; in this
@ -1357,23 +1357,29 @@ Calculating Integrated Cross Sections
The first aspect of using |sab| tables is calculating cross sections to replace
the data that would normally appear on the incident neutron data, which do not
account for thermal binding effects. For incoherent elastic and inelastic
scattering, the cross sections are stored as linearly interpolable functions on
a specified energy grid. For coherent elastic data, the cross section can be
expressed as
account for thermal binding effects. For incoherent inelastic scattering, the
cross section is stored as a linearly interpolable function on a specified
energy grid. For coherent elastic data, the cross section can be expressed as
.. math::
:label: coherent-elastic-xs
\sigma(E) = \frac{\sigma_c}{E} \sum_{E_i < E} f_i e^{-4WE_i}
\sigma(E) = \frac{1}{E} \sum_{E_i < E} s_i
where :math:`\sigma_c` is the effective bound coherent scattering cross section,
:math:`W` is the effective Debye-Waller coefficient, :math:`E_i` are the
energies of the Bragg edges, and :math:`f_i` are related to crystallographic
structure factors. Since the functional form of the cross section is just 1/E
and the proportionality constant changes only at Bragg edges, the
proportionality constants are stored and then the cross section can be
calculated analytically based on equation :eq:`coherent-elastic-xs`.
where :math:`E_i` are the energies of the Bragg edges and :math:`s_i` are
related to crystallographic structure factors. Since the functional form of the
cross section is just 1/E and the proportionality constant changes only at Bragg
edges, the proportionality constants are stored and then the cross section can
be calculated analytically based on equation :eq:`coherent-elastic-xs`. For
incoherent elastic data, the cross section can be expressed as
.. math::
:label: incoherent-elastic-xs
\sigma(E) = \frac{\sigma_b}{2} \left( \frac{1 - e^{-4EW'}}{2EW'} \right)
where :math:`\sigma_b` is the characteristic bound cross section and :math:`W'`
is the Debye-Waller integral divided by the atomic mass.
Outgoing Angle for Coherent Elastic Scattering
----------------------------------------------
@ -1388,7 +1394,7 @@ scatter then neutron is given by
.. math::
:label: coherent-elastic-probability
\frac{f_i e^{-4WE_i}}{\sum_j f_j e^{-4WE_j}}.
\frac{s_i}{\sum_j s_j}.
After a Bragg edge has been sampled, the cosine of the angle of scattering is
given analytically by
@ -1400,20 +1406,35 @@ given analytically by
where :math:`E_i` is the energy of the Bragg edge that scattered the neutron.
.. _incoherent elastic angle:
Outgoing Angle for Incoherent Elastic Scattering
------------------------------------------------
For incoherent elastic scattering, the probability distribution for the cosine
of the angle of scattering is represent as a series of equally-likely discrete
For incoherent elastic scattering, OpenMC has two methods for calculating the
cosine of the angle of scattering. The first method uses the Debye-Waller
integral, :math:`W'`, and the characteristic bound cross section as given
directly in an ENDF-6 formatted file. In this case, the cosine of the angle of
scattering can be sampled by inverting equation 7.4 from the `ENDF-6 Format
Manual <endf102>`_:
.. math::
:label: incoherent-elastic-mu-exact
\mu = \frac{1}{c} \log \left( 1 + \xi \left( e^{2c} - 1 \right) \right) - 1
where :math:`\xi` is a random number sampled on unit interval and :math:`c =
2EW'`. In the second method, the probability distribution for the cosine of the
angle of scattering is represented as a series of equally-likely discrete
cosines :math:`\mu_{i,j}` for each incoming energy :math:`E_i` on the thermal
elastic energy grid. First the outgoing angle bin :math:`j` is sampled. Then, if
the incoming energy of the neutron satisfies :math:`E_i < E < E_{i+1}` the final
cosine is
the incoming energy of the neutron satisfies :math:`E_i < E < E_{i+1}` the
cosine of the angle of scattering is
.. math::
:label: incoherent-elastic-angle
\mu = \mu_{i,j} + f (\mu_{i+1,j} - \mu_{i,j})
\mu' = \mu_{i,j} + f (\mu_{i+1,j} - \mu_{i,j})
where the interpolation factor is defined as
@ -1422,6 +1443,30 @@ where the interpolation factor is defined as
f = \frac{E - E_i}{E_{i+1} - E_i}.
To better represent the true, continuous nature of the cosine distribution, the
sampled value of :math:`mu'` is then "smeared" based on the neighboring values.
First, values of :math:`\mu` are calculated for outgoing angle bins :math:`j-1`
and :math:`j+1`:
.. math::
:label: incoherent-elastic-smear1
\mu_\text{left} = \mu_{i,j-1} + f (\mu_{i+1,j-1} - \mu_{i,j-1}) \\
\mu_\text{right} = \mu_{i,j+1} + f (\mu_{i+1,j+1} - \mu_{i,j+1}).
Then, a final cosine is calculated as:
.. math::
:label: incoherent-elastic-smear2
\mu = \mu' + \min (\mu - \mu_\text{left}, \mu + \mu_\text{right} ) \cdot
\left( \xi - \frac{1}{2} \right)
where :math:`\xi` is again a random number sampled on the unit interval. Care
must be taken to ensure that :math:`\mu` does not fall outside the interval
:math:`[-1,1]`.
Outgoing Energy and Angle for Inelastic Scattering
--------------------------------------------------
@ -1492,7 +1537,10 @@ which angular distribution data to use. Like the linear-linear interpolation
case in Law 61, the angular distribution closest to the sampled value of the
cumulative distribution function for the outgoing energy is utilized. The
actual algorithm utilized to sample the outgoing angle is shown in equation
:eq:`inelastic-angle`.
:eq:`inelastic-angle`. As in the case of incoherent elastic scattering with
discrete cosine bins, the sampled cosine is :ref:`smeared <incoherent elastic
angle>` over neighboring angle bins to better approximate a continuous
distribution.
.. _probability_tables:
@ -1660,7 +1708,7 @@ another.
.. _PREPRO: http://www-nds.iaea.org/ndspub/endf/prepro/
.. _ENDF-6 Format: https://www.oecd-nea.org/dbdata/data/manual-endf/endf102.pdf
.. _endf102: https://www.oecd-nea.org/dbdata/data/manual-endf/endf102.pdf
.. _Monte Carlo Sampler: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-9721.pdf

8
src/lattice.cpp
View file

@ -895,8 +895,12 @@ HexLattice::get_indices(Position r, Direction u) const
Position r_t = get_local_position(r, i_xyz);
// calculate distance
double d = r_t.x*r_t.x + r_t.y*r_t.y;
// check for coincidence
bool on_boundary = coincident(d, d_min);
// check for coincidence. Because the numerical error incurred
// in hex geometry is higher than other geometries, the relative
// coincidence is checked here so that coincidence is successfully
// detected on large hex lattice with particles far from the origin
// which have rounding errors larger than the FP_COINCIDENT thresdhold.
bool on_boundary = coincident(1.0, d_min/d);
if (d < d_min || on_boundary) {
// normalize r_t and find dot product
r_t /= std::sqrt(d);

90
src/secondary_thermal.cpp
View file

@ -101,14 +101,31 @@ IncoherentElasticAEDiscrete::sample(double E_in, double& E_out, double& mu,
// incoming energies.
// Sample outgoing cosine bin
int k = prn(prn_seed) * mu_out_.shape()[1];
int n_mu = mu_out_.shape()[1];
int k = prn(prn_seed) * n_mu;
// Determine outgoing cosine corresponding to E_in[i] and E_in[i+1]
double mu_ik = mu_out_(i, k);
double mu_i1k = mu_out_(i+1, k);
// Rather than use the sampled discrete mu directly, it is smeared over
// a bin of width 0.5*min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the
// discrete mu value itself.
// Cosine of angle between incoming and outgoing neutron
mu = (1 - f)*mu_ik + f*mu_i1k;
// Interpolate kth mu value between distributions at energies i and i+1
mu = mu_out_(i, k) + f*(mu_out_(i+1, k) - mu_out_(i, k));
// Inteprolate (k-1)th mu value between distributions at energies i and i+1.
// When k==0, pick a value that will smear the cosine out to a minimum of -1.
double mu_left = (k == 0) ?
-1.0 - (mu + 1.0) :
mu_out_(i, k-1) + f*(mu_out_(i+1, k-1) - mu_out_(i, k-1));
// Inteprolate (k+1)th mu value between distributions at energies i and i+1.
// When k is the last discrete value, pick a value that will smear the cosine
// out to a maximum of 1.
double mu_right = (k == n_mu - 1) ?
1.0 + (1.0 - mu) :
mu_out_(i, k+1) + f*(mu_out_(i+1, k+1) - mu_out_(i, k+1));
// Smear cosine
mu += std::min(mu - mu_left, mu_right - mu)*(prn() - 0.5);
// Energy doesn't change in elastic scattering
E_out = E_in;
@ -240,25 +257,12 @@ IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu,
double f;
get_energy_index(energy_, E_in, i, f);
// Sample between ith and [i+1]th bin
int l = f > prn(prn_seed) ? i + 1 : i;
// Determine endpoints on grid i
auto n = distribution_[i].e_out.size();
double E_i_1 = distribution_[i].e_out[0];
double E_i_J = distribution_[i].e_out[n - 1];
// Determine endpoints on grid i + 1
n = distribution_[i + 1].e_out.size();
double E_i1_1 = distribution_[i + 1].e_out[0];
double E_i1_J = distribution_[i + 1].e_out[n - 1];
double E_1 = E_i_1 + f * (E_i1_1 - E_i_1);
double E_J = E_i_J + f * (E_i1_J - E_i_J);
// Pick closer energy based on interpolation factor
int l = f > 0.5 ? i + 1 : i;
// Determine outgoing energy bin
// (First reset n_energy_out to the right value)
n = distribution_[l].n_e_out;
auto n = distribution_[l].n_e_out;
double r1 = prn(prn_seed);
double c_j = distribution_[l].e_out_cdf[0];
double c_j1;
@ -289,11 +293,12 @@ IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu,
2.0*frac*(r1 - c_j))) - p_l_j) / frac;
}
// Now interpolate between incident energy bins i and i + 1
if (l == i) {
E_out = E_1 + (E_out - E_i_1) * (E_J - E_1) / (E_i_J - E_i_1);
// Adjustment of outgoing energy
double E_l = energy_[l];
if (E_out < 0.5*E_l) {
E_out *= 2.0*E_in/E_l - 1.0;
} else {
E_out = E_1 + (E_out - E_i1_1) * (E_J - E_1) / (E_i1_J - E_i1_1);
E_out += E_in - E_l;
}
// Sample outgoing cosine bin
@ -301,31 +306,28 @@ IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu,
std::size_t k = prn(prn_seed) * n_mu;
// Rather than use the sampled discrete mu directly, it is smeared over
// a bin of width min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the
// a bin of width 0.5*min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the
// discrete mu value itself.
const auto& mu_l = distribution_[l].mu;
f = (r1 - c_j)/(c_j1 - c_j);
// Determine (k-1)th mu value
double mu_left;
if (k == 0) {
mu_left = -1.0;
} else {
mu_left = mu_l(j, k-1) + f*(mu_l(j+1, k-1) - mu_l(j, k-1));
}
// Determine kth mu value
// Interpolate kth mu value between distributions at energies j and j+1
mu = mu_l(j, k) + f*(mu_l(j+1, k) - mu_l(j, k));
// Determine (k+1)th mu value
double mu_right;
if (k == n_mu - 1) {
mu_right = 1.0;
} else {
mu_right = mu_l(j, k+1) + f*(mu_l(j+1, k+1) - mu_l(j, k+1));
}
// Inteprolate (k-1)th mu value between distributions at energies j and j+1.
// When k==0, pick a value that will smear the cosine out to a minimum of -1.
double mu_left = (k == 0) ?
mu_left = -1.0 - (mu + 1.0) :
mu_left = mu_l(j, k-1) + f*(mu_l(j+1, k-1) - mu_l(j, k-1));
// Smear angle
// Inteprolate (k+1)th mu value between distributions at energies j and j+1.
// When k is the last discrete value, pick a value that will smear the cosine
// out to a maximum of 1.
double mu_right = (k == n_mu - 1) ?
mu_right = 1.0 + (1.0 - mu) :
mu_right = mu_l(j, k+1) + f*(mu_l(j+1, k+1) - mu_l(j, k+1));
// Smear cosine
mu += std::min(mu - mu_left, mu_right - mu)*(prn(prn_seed) - 0.5);
}