From d29923f6bff3dd14fc882268fc8e94bb2f0d7f8a Mon Sep 17 00:00:00 2001 From: amandalund Date: Mon, 9 Jul 2018 21:33:48 -0500 Subject: [PATCH] Address #1023 comments --- docs/source/io_formats/nuclear_data.rst | 4 +- docs/source/methods/photon_physics.rst | 62 ++++++++++++++----------- openmc/data/photon.py | 11 +++-- src/constants.F90 | 2 +- src/input_xml.F90 | 2 +- src/material_header.F90 | 4 +- src/photon_physics.F90 | 28 +++++------ src/physics.F90 | 16 +++---- src/tallies/tally.F90 | 1 - 9 files changed, 69 insertions(+), 61 deletions(-) diff --git a/docs/source/io_formats/nuclear_data.rst b/docs/source/io_formats/nuclear_data.rst index 9bc9d295b..c5d8ac347 100644 --- a/docs/source/io_formats/nuclear_data.rst +++ b/docs/source/io_formats/nuclear_data.rst @@ -181,7 +181,7 @@ Incident Photon Data :Datasets: - **I** (*double*) -- Mean excitation energy in [eV] - **energy** (*double[]*) -- Energies in [eV] - - **s_collision** (*double[]*) -- Collisiong stopping power in [eV-cm\ :sup:`2`\ /g] + - **s_collision** (*double[]*) -- Collision stopping power in [eV-cm\ :sup:`2`\ /g] - **s_radiative** (*double[]*) -- Radiative stopping power in [eV-cm\ :sup:`2`\ /g] **//subshells/** @@ -199,7 +199,7 @@ Incident Photon Data :Attributes: - **threshold_idx** (*int*) -- Index on the energy - grid that the reaction threshold + grid of the reaction threshold ------------------------------- Thermal Neutron Scattering Data diff --git a/docs/source/methods/photon_physics.rst b/docs/source/methods/photon_physics.rst index 836c854b0..bf4892a78 100644 --- a/docs/source/methods/photon_physics.rst +++ b/docs/source/methods/photon_physics.rst @@ -5,7 +5,7 @@ Photon Physics ============== Photons, being neutral particles, behave much in the same manner as neutrons, -traveling in straight lines and experiencing occassional collisions which change +traveling in straight lines and experiencing occasional collisions which change their energy and direction. Photons undergo four basic interactions as they pass through matter: coherent (Rayleigh) scattering, incoherent (Compton) scattering, photoelectric effect, and pair/triplet production. Photons with energy in the @@ -31,11 +31,11 @@ cross section is .. math:: :label: thomson - \frac{d\sigma}{d\mu} = \pi r_0^2 ( 1 + \mu^2 ) + \frac{d\sigma}{d\mu} = \pi r_e^2 ( 1 + \mu^2 ) -where :math:`\mu` is the cosine of the scattering angle and :math:`r_0` is the -classical radius of the electron. Thomson scattering can generally occur when -the photon energy is much less than rest mass energy of the particle. +where :math:`\mu` is the cosine of the scattering angle and :math:`r_e` is the +classical electron radius. Thomson scattering can generally occur when the +photon energy is much less than the rest mass energy of the particle. In practice, most elastic scattering of photons off electrons happens not with free electrons but those bound in atoms. This process is known as Rayleigh @@ -50,30 +50,33 @@ The differential cross section for Rayleigh scattering is given by .. math:: :label: coherent-xs - \frac{d\sigma(E,E',\mu)}{d\mu} = \pi r_0^2 ( 1 + \mu^2 ) \left [ ( F(x, Z) + - F'(E) )^2 + F''(E)^2 \right ] + \frac{d\sigma(E,E',\mu)}{d\mu} &= \pi r_e^2 ( 1 + \mu^2 )~\left| F(x,Z) + + F' + iF'' \right|^2 \\ + &= \pi r_e^2 ( 1 + \mu^2 ) \left [ ( F(x,Z) + + F'(E) )^2 + F''(E)^2 \right ] where :math:`F(x,Z)` is a form factor as a function of the momentum transfer -:math:`x` and the atomic number :math:`Z` and :math:`F' + iF''` is a factor that -accounts for `anomalous scattering`_ which can occur near absorption edges. In a -Monte Carlo simulation, when coherent scattering occurs, we only need to sample -the scattering angle using the differential cross section in :eq:`coherent-xs` -since the energy of the photon does not change. In OpenMC, anomalous scattering -is ignored such that differential cross section becomes +:math:`x` and the atomic number :math:`Z` and the term :math:`F' + iF''` +accounts for `anomalous scattering`_ which can occur near absorption edges. In +a Monte Carlo simulation, when coherent scattering occurs, we only need to +sample the scattering angle using the differential cross section in +:eq:`coherent-xs` since the energy of the photon does not change. In OpenMC, +anomalous scattering is ignored such that the differential cross section +becomes .. math:: :label: coherent-xs-openmc - \frac{d\sigma(E,E',\mu)}{d\mu} = \pi r_0^2 ( 1 + \mu^2 ) F(x, Z)^2 + \frac{d\sigma(E,E',\mu)}{d\mu} = \pi r_e^2 ( 1 + \mu^2 ) F(x, Z)^2 -To construct a proper probability density, we need to normalize the differential -cross section in :eq:`coherent-xs-openmc` by the integrated coherent scattering -cross section: +To construct a proper probability density, we need to normalize the +differential cross section in :eq:`coherent-xs-openmc` by the integrated +coherent scattering cross section: .. math:: :label: coherent-pdf-1 - p(\mu) d\mu = \frac{\pi r_0^2}{\sigma(E)} ( 1 + \mu^2 ) F(x, Z)^2 d\mu. + p(\mu) d\mu = \frac{\pi r_e^2}{\sigma(E)} ( 1 + \mu^2 ) F(x, Z)^2 d\mu. Since the form factor is given in terms of the momentum transfer, it is more convenient to change variables of the probability density to :math:`x^2`. The @@ -84,23 +87,24 @@ momentum transfer is traditionally expressed as x = \kappa \alpha \sqrt{1 - \mu} -where the coefficient :math:`\kappa` can be shown to be +where :math:`\alpha` is the ratio of the photon energy to the electron rest +mass, and the coefficient :math:`\kappa` can be shown to be .. math:: :label: kappa \kappa = \frac{m_e c^2}{\sqrt{2}hc} \approx 29.14329, -:math:`m_e` is the mass of the electron, :math:`c` is the speed of light +where :math:`m_e` is the mass of the electron, :math:`c` is the speed of light in a vacuum, and :math:`h` is Planck's constant. Using :eq:`momentum-transfer`, -we have that :math:`\mu = 1 - [x/(\kappa\alpha)]^2` and :math:`d\mu/dx^2 = +we have :math:`\mu = 1 - [x/(\kappa\alpha)]^2` and :math:`d\mu/dx^2 = -1/(\kappa\alpha)^2`. The probability density in :math:`x^2` is .. math:: :label: coherent-pdf-x2 p(x^2) dx^2 = p(\mu) \left | \frac{d\mu}{dx^2} \right | dx^2 = \frac{2\pi - r_0^2 A(\bar{x}^2,Z)}{(\kappa\alpha)^2 \sigma(E)} \left ( + r_e^2 A(\bar{x}^2,Z)}{(\kappa\alpha)^2 \sigma(E)} \left ( \frac{1 + \mu^2}{2} \right ) \left ( \frac{F(x, Z)^2}{A(\bar{x}^2, Z)} \right ) dx^2 where :math:`\bar{x}` is the maximum value of :math:`x` that occurs for @@ -116,7 +120,7 @@ and :math:`A(x^2, Z)` is the integral of the square of the form factor: .. math:: :label: coherent-int-ff - A(x^2, Z) = \int_0^{x^2} F(\chi, Z)^2 d\chi^2. + A(x^2, Z) = \int_0^{x^2} F(x,Z)^2 dx^2. As you see, we have multiplied and divided the probability density by the integral of the squared form factor so that the density in :eq:`coherent-pdf-x2` @@ -127,7 +131,7 @@ run-time to do a table search on the cumulative distribution function: .. math:: :label: coherent-form-factor-cdf - \frac{\int_0^{x^2} F(\chi,Z)^2 d\chi^2}{\int_0^{\bar{x}^2} F(x,Z)^2 dx^2} + \frac{\int_0^{x^2} F(x,Z)^2 dx^2}{\int_0^{\bar{x}^2} F(x,Z)^2 dx^2} Once a trial :math:`x^2` value has been selected, we can calculate :math:`\mu` and perform rejection sampling using the Thomson scattering differential cross @@ -162,7 +166,7 @@ the two authors who discovered it: .. math:: :label: klein-nishina - \frac{d\sigma_{KN}}{d\mu} = \pi r_0^2 \left ( \frac{\alpha'}{\alpha} \right + \frac{d\sigma_{KN}}{d\mu} = \pi r_e^2 \left ( \frac{\alpha'}{\alpha} \right ) \left [ \frac{\alpha'}{\alpha} + \frac{\alpha}{\alpha'} + \mu^2 - 1 \right ] @@ -188,7 +192,7 @@ differential cross section for incoherent scattering is given by .. math:: :label: incoherent-xs - \frac{d\sigma}{d\mu} = \frac{d\sigma_{KN}}{d\mu} S(x,Z) = \pi r_0^2 \left ( + \frac{d\sigma}{d\mu} = \frac{d\sigma_{KN}}{d\mu} S(x,Z) = \pi r_e^2 \left ( \frac{\alpha'}{\alpha} \right )^2 \left [ \frac{\alpha'}{\alpha} + \frac{\alpha}{\alpha'} + \mu^2 - 1 \right ] S(x,Z) @@ -213,6 +217,10 @@ Doppler Energy Broadening LA-UR-04-0487_ and LA-UR-04-0488_ +Compton Electrons ++++++++++++++++++ + + Photoelectric Effect -------------------- @@ -253,7 +261,7 @@ Thick-Target Bremsstrahlung Approximation +++++++++++++++++++++++++++++++++++++++++ -.. _Koblinger: http://www.tandfonline.com/doi/abs/10.13182/NSE75-A26646 +.. _Koblinger: https://doi.org/10.13182/NSE75-A26663 .. _anomalous scattering: http://pd.chem.ucl.ac.uk/pdnn/diff1/anomscat.htm diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 4ed6ef8c8..f0d8a2f02 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -110,10 +110,11 @@ class AtomicRelaxation(EqualityMixin): """Atomic relaxation data. This class stores the binding energy, number of electrons, and electron - transitions possible from ioniziation for each subshell with an atom. All of - the data originates from an ENDF-6 atomic relaxation sub-library - (NSUB=6). Instances of this class are not normally instantiated directly but - rather created using the factory method :math:`AtomicRelaxation.from_endf`. + transitions possible from ioniziation for each electron subshell of an + atom. All of the data originates from an ENDF-6 atomic relaxation + sub-library (NSUB=6). Instances of this class are not normally instantiated + directly but rather created using the factory method + :math:`AtomicRelaxation.from_endf`. Parameters ---------- @@ -353,7 +354,7 @@ class IncidentPhoton(EqualityMixin): Dictionary of bremsstrahlung DCS data with keys 'electron_energy' (incident electron kinetic energy values in [eV]), 'photon_energy' (ratio of the energy of the emitted photon to the incident electron - kinetic energy), and 'dcs' (cross sectin values in [b]). The cross + kinetic energy), and 'dcs' (cross section values in [b]). The cross sections are in scaled form: :math:`(\beta^2/Z^2) E_k (d\sigma/dE_k)`, where :math:`E_k` is the energy of the emitted photon. compton_profiles : dict diff --git a/src/constants.F90 b/src/constants.F90 index e1b298d5a..2ee5fb2cd 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -78,7 +78,7 @@ module constants MASS_NEUTRON = 1.00866491588_8, & ! mass of a neutron in amu MASS_NEUTRON_EV = 939.5654133e6_8, & ! mass of a neutron in eV/c^2 MASS_PROTON = 1.007276466879_8, & ! mass of a proton in amu - MASS_ELECTRON = 0.5109989461e6_8, & ! electron mass energy equivalent in eV + MASS_ELECTRON_EV = 0.5109989461e6_8, & ! electron mass energy equivalent in eV/c^2 FINE_STRUCTURE = 137.035999139_8, & ! inverse fine structure constant PLANCK_C = 1.2398419739062977e4_8,& ! Planck's constant times c in eV-Angstroms AMU = 1.660539040e-27_8, & ! 1 amu in kg diff --git a/src/input_xml.F90 b/src/input_xml.F90 index c66cc2660..b0831bd90 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2610,7 +2610,7 @@ contains do l = 1, filt % n_bins if (filt % particles(l) /= NEUTRON) then call warning("Particle filter other than NEUTRON used with & - &photon transport turn off. All tallies for particle & + &photon transport turned off. All tallies for particle & &type " // trim(to_str(filt % particles(l))) // " will have no scores") end if end do diff --git a/src/material_header.F90 b/src/material_header.F90 index 9efe932d4..03ecf69d9 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -809,7 +809,7 @@ contains ! Issy-les-Moulineaux, France (2011). if (positron_) then do i = 1, n_e - t = log(ONE + 1.0e6_8*ttb_e_grid(i)/(Z_eq_sq*MASS_ELECTRON)) + t = log(ONE + 1.0e6_8*ttb_e_grid(i)/(Z_eq_sq*MASS_ELECTRON_EV)) r = ONE - exp(-1.2359e-1_8*t + 6.1274e-2_8*t**2 - 3.1516e-2_8*t**3 + & 7.7446e-3_8*t**4 - 1.0595e-3_8*t**5 + 7.0568e-5_8*t**6 - & 1.808e-6_8*t**7) @@ -846,7 +846,7 @@ contains x = x_l + (k - k_l) * (x_r - x_l) / (k_r - k_l) ! Ratio of the velocity of the charged particle to the speed of light - beta = sqrt(e*(e + TWO*MASS_ELECTRON)) / (e + MASS_ELECTRON) + beta = sqrt(e*(e + TWO*MASS_ELECTRON_EV)) / (e + MASS_ELECTRON_EV) ! Compute the integrand of the PDF f(j) = x / (beta**2 * stopping_power(j) * w) diff --git a/src/photon_physics.F90 b/src/photon_physics.F90 index 79388a8ce..3a84a9a01 100644 --- a/src/photon_physics.F90 +++ b/src/photon_physics.F90 @@ -107,20 +107,20 @@ contains ! Note that the parameter used here does not correspond exactly to the ! momentum transfer q in ENDF-102 Eq. (27.2). Rather, this is the ! parameter as defined by Hubbell, where the actual data comes from - x = MASS_ELECTRON/PLANCK_C*alpha*sqrt(HALF*(ONE - mu)) + x = MASS_ELECTRON_EV/PLANCK_C*alpha*sqrt(HALF*(ONE - mu)) ! Calculate S(x, Z) and S(x_max, Z) form_factor_x = el % incoherent_form_factor % evaluate(x) if (form_factor_xmax == ZERO) then form_factor_xmax = el % incoherent_form_factor % evaluate(& - MASS_ELECTRON/PLANCK_C*alpha) + MASS_ELECTRON_EV/PLANCK_C*alpha) end if ! Perform rejection on form factor if (prn() < form_factor_x / form_factor_xmax) then if (use_doppler_) then call compton_doppler(el, alpha, mu, e_out, i_shell) - alpha_out = e_out/MASS_ELECTRON + alpha_out = e_out/MASS_ELECTRON_EV else i_shell = 0 end if @@ -168,16 +168,16 @@ contains e_b = el % binding_energy(i_shell) ! Determine p_z,max - e = alpha*MASS_ELECTRON + e = alpha*MASS_ELECTRON_EV if (e < e_b) then - e_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON + e_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON_EV exit end if pz_max = -FINE_STRUCTURE*(e_b - (e - e_b)*alpha*(ONE - mu)) / & sqrt(TWO*e*(e - e_b)*(ONE - mu) + e_b**2) if (pz_max < ZERO) then - e_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON + e_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON_EV exit end if @@ -230,7 +230,7 @@ contains quad = b**2 - FOUR*a*c if (quad < 0) then - e_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON + e_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON_EV exit end if quad = sqrt(quad) @@ -280,7 +280,7 @@ contains do ! Determine maximum value of x^2 - x2_max = (MASS_ELECTRON/PLANCK_C*alpha)**2 + x2_max = (MASS_ELECTRON_EV/PLANCK_C*alpha)**2 ! Determine F(x^2_max, Z) F_max = el % coherent_int_form_factor % evaluate(x2_max) @@ -517,20 +517,20 @@ contains end do ! Compute the kinetic energy of the electron and the positron - E_electron = (alpha*e - ONE)*MASS_ELECTRON - E_positron = (alpha*(ONE - e) - ONE)*MASS_ELECTRON + E_electron = (alpha*e - ONE)*MASS_ELECTRON_EV + E_positron = (alpha*(ONE - e) - ONE)*MASS_ELECTRON_EV ! Sample the scattering angle of the electron. The cosine of the polar ! angle of the direction relative to the incident photon is sampled from ! p(mu) = C/(1 - beta*mu)^2 using the inverse transform method. - beta = sqrt(E_electron*(E_electron + TWO*MASS_ELECTRON)) & - / (E_electron + MASS_ELECTRON) + beta = sqrt(E_electron*(E_electron + TWO*MASS_ELECTRON_EV)) & + / (E_electron + MASS_ELECTRON_EV) rn = TWO*prn() - ONE mu_electron = (rn + beta)/(rn*beta + ONE) ! Sample the scattering angle of the positron - beta = sqrt(E_positron*(E_positron + TWO*MASS_ELECTRON)) & - / (E_positron + MASS_ELECTRON) + beta = sqrt(E_positron*(E_positron + TWO*MASS_ELECTRON_EV)) & + / (E_positron + MASS_ELECTRON_EV) rn = TWO*prn() - ONE mu_positron = (rn + beta)/(rn*beta + ONE) diff --git a/src/physics.F90 b/src/physics.F90 index 02872ef03..4f2196d7d 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -195,7 +195,7 @@ contains p % event_nuclide = i_element ! Calculate photon energy over electron rest mass equivalent - alpha = p % E/MASS_ELECTRON + alpha = p % E/MASS_ELECTRON_EV ! For tallying purposes, this routine might be called directly. In that ! case, we need to sample a reaction via the cutoff variable @@ -226,7 +226,7 @@ contains end if ! Create Compton electron - E_electron = (alpha - alpha_out)*MASS_ELECTRON - e_b + E_electron = (alpha - alpha_out)*MASS_ELECTRON_EV - e_b mu_electron = (alpha - alpha_out*mu) & / sqrt(alpha**2 + alpha_out**2 - TWO*alpha*alpha_out*mu) phi = TWO*PI*prn() @@ -241,7 +241,7 @@ contains end if phi = phi + PI - p % E = alpha_out*MASS_ELECTRON + p % E = alpha_out*MASS_ELECTRON_EV p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu, phi) p % event_MT = INCOHERENT return @@ -274,8 +274,8 @@ contains SAMPLE_MU: do r = prn() if (FOUR * (ONE - r) * r >= prn()) then - rel_vel = sqrt(E_electron * (E_electron + TWO * MASS_ELECTRON))& - / (E_electron + MASS_ELECTRON) + rel_vel = sqrt(E_electron * (E_electron + TWO * MASS_ELECTRON_EV))& + / (E_electron + MASS_ELECTRON_EV) mu = (TWO * r + rel_vel - ONE) / & (TWO * rel_vel * r - rel_vel + ONE) exit SAMPLE_MU @@ -354,7 +354,7 @@ contains ! energy locally (electron_treatment = ELECTRON_LED) or creates secondary ! bremsstrahlung photons from electron deflections with charged particles ! (electron_treatment = ELECTRON_TTB). Two annihilation photons of energy -! MASS_ELECTRON (0.511 MeV) are created and travel in opposite directions. +! MASS_ELECTRON_EV (0.511 MeV) are created and travel in opposite directions. !=============================================================================== subroutine sample_positron_reaction(p) @@ -380,8 +380,8 @@ contains uvw(3) = sqrt(ONE - mu*mu)*sin(phi) ! Create annihilation photon pair traveling in opposite directions - call p % create_secondary( uvw, MASS_ELECTRON, PHOTON, .true.) - call p % create_secondary(-uvw, MASS_ELECTRON, PHOTON, .true.) + call p % create_secondary( uvw, MASS_ELECTRON_EV, PHOTON, .true.) + call p % create_secondary(-uvw, MASS_ELECTRON_EV, PHOTON, .true.) p % E = ZERO p % alive = .false. diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 0c9448b4e..8ddd42bce 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -298,7 +298,6 @@ contains if (t % find_filter(FILTER_ENERGYOUT) > 0) then ! Normally, we only need to make contributions to one scoring ! bin. However, in the case of fission, since multiple fission - ! neutrons were emitted with different energies, multiple ! outgoing energy bins may have been scored to. The following ! logic treats this special case and results to multiple bins