diff --git a/CMakeLists.txt b/CMakeLists.txt index 68e5bdaef1..2ccc457f79 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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 diff --git a/docs/source/methods/neutron_physics.rst b/docs/source/methods/neutron_physics.rst index 384e79bdd6..e898e63b0e 100644 --- a/docs/source/methods/neutron_physics.rst +++ b/docs/source/methods/neutron_physics.rst @@ -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 `_ 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 `_: + +.. 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 ` 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 diff --git a/src/lattice.cpp b/src/lattice.cpp index 6f10ed3584..092ff3d116 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -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); diff --git a/src/secondary_thermal.cpp b/src/secondary_thermal.cpp index 13b5746fc6..f063e5685b 100644 --- a/src/secondary_thermal.cpp +++ b/src/secondary_thermal.cpp @@ -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); }