From f1a0d7821f30ead22df69b7a8013132435f0534d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 15 Nov 2019 16:26:55 -0600 Subject: [PATCH 1/7] Fix RPATH when CMAKE_INSTALL_LIBDIR is not just lib --- CMakeLists.txt | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 68e5bdaef1..8c02068a3b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -152,6 +152,8 @@ target_compile_definitions(gsl INTERFACE GSL_THROW_ON_CONTRACT_VIOLATION) # RPATH information #=============================================================================== +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 +166,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 +360,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 From 3379e67bb7fb1be4dd367396918e3ad023247653 Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Tue, 19 Nov 2019 17:00:42 -0500 Subject: [PATCH 2/7] use relative distances for coincidence test in hex lattice --- src/lattice.cpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) 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); From aed61a4573cf7bf6dba53619cc1e2ffb05942c36 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 20 Nov 2019 15:45:26 -0600 Subject: [PATCH 3/7] Add comment regarding GNUInstallDirs --- CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 8c02068a3b..2ccc457f79 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -152,6 +152,7 @@ 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 From d3251e62f86c59037c519e681864978b5bcf6aea Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 27 Nov 2019 07:51:04 -0600 Subject: [PATCH 4/7] Fix cosine smearing for S(a,b) --- src/secondary_thermal.cpp | 62 ++++++++++++++++++++++++--------------- 1 file changed, 38 insertions(+), 24 deletions(-) diff --git a/src/secondary_thermal.cpp b/src/secondary_thermal.cpp index 8c9cc69676..3089c7a50f 100644 --- a/src/secondary_thermal.cpp +++ b/src/secondary_thermal.cpp @@ -98,14 +98,31 @@ IncoherentElasticAEDiscrete::sample(double E_in, double& E_out, double& mu) cons // incoming energies. // Sample outgoing cosine bin - int k = prn() * mu_out_.shape()[1]; + int n_mu = mu_out_.shape()[1]; + int k = prn() * 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; @@ -296,31 +313,28 @@ IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu) const std::size_t k = prn() * 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() - 0.5); } From eb12d612e943487e157f3a956002bdd765b33479 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 27 Nov 2019 11:11:01 -0600 Subject: [PATCH 5/7] Change outgoing energy adjustment for incoherent inelastic to match MCNP6 --- src/secondary_thermal.cpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/secondary_thermal.cpp b/src/secondary_thermal.cpp index 3089c7a50f..207d4d0d83 100644 --- a/src/secondary_thermal.cpp +++ b/src/secondary_thermal.cpp @@ -252,8 +252,8 @@ IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu) const double f; get_energy_index(energy_, E_in, i, f); - // Sample between ith and [i+1]th bin - int l = f > prn() ? i + 1 : i; + // Pick closer energy based on interpolation factor + int l = f > 0.5 ? i + 1 : i; // Determine endpoints on grid i auto n = distribution_[i].e_out.size(); @@ -301,11 +301,12 @@ IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu) const 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 From 428e6248ac1ac6fbe162ef553bbb00908f390c6a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 27 Nov 2019 14:23:16 -0600 Subject: [PATCH 6/7] Remove unused variables in incoherent inelastic sampling --- src/secondary_thermal.cpp | 15 +-------------- 1 file changed, 1 insertion(+), 14 deletions(-) diff --git a/src/secondary_thermal.cpp b/src/secondary_thermal.cpp index 207d4d0d83..4642e8ef84 100644 --- a/src/secondary_thermal.cpp +++ b/src/secondary_thermal.cpp @@ -255,22 +255,9 @@ IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu) const // Pick closer energy based on interpolation factor int l = f > 0.5 ? 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); - // 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(); double c_j = distribution_[l].e_out_cdf[0]; double c_j1; From 20bd6b5d35f623fb08933db7652daa854a4dda76 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 27 Nov 2019 14:23:34 -0600 Subject: [PATCH 7/7] Update documentation for S(a,b) --- docs/source/methods/neutron_physics.rst | 94 +++++++++++++++++++------ 1 file changed, 71 insertions(+), 23 deletions(-) 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