diff --git a/.gitignore b/.gitignore index 7eb8f286ac..3832322966 100644 --- a/.gitignore +++ b/.gitignore @@ -69,6 +69,7 @@ scripts/G4EMLOW*/ # Images *.ppm *.voxel +*.vti # PyCharm project configuration files .idea diff --git a/docs/source/methods/photon_physics.rst b/docs/source/methods/photon_physics.rst index 4d58ee0db2..028a71de72 100644 --- a/docs/source/methods/photon_physics.rst +++ b/docs/source/methods/photon_physics.rst @@ -85,26 +85,26 @@ momentum transfer is traditionally expressed as .. math:: :label: momentum-transfer - x = \kappa \alpha \sqrt{1 - \mu} + x = a k \sqrt{1 - \mu} -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 +where :math:`k` is the ratio of the photon energy to the electron rest +mass, and the coefficient :math:`a` can be shown to be .. math:: - :label: kappa + :label: omega - \kappa = \frac{m_e c^2}{\sqrt{2}hc} \approx 29.14329, + a = \frac{m_e c^2}{\sqrt{2}hc} \approx 29.14329~\unicode{x212B}, 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 :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 +we have :math:`\mu = 1 - [x/(ak)]^2` and :math:`d\mu/dx^2 = +-1/(ak)^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_e^2 A(\bar{x}^2,Z)}{(\kappa\alpha)^2 \sigma(E)} \left ( + r_e^2 A(\bar{x}^2,Z)}{(ak)^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 @@ -113,7 +113,7 @@ where :math:`\bar{x}` is the maximum value of :math:`x` that occurs for .. math:: :label: xmax - \bar{x} = \kappa \alpha \sqrt{2} = \frac{m_e c^2}{hc} \alpha, + \bar{x} = a k \sqrt{2} = \frac{m_e c^2}{hc} k, and :math:`A(x^2, Z)` is the integral of the square of the form factor: @@ -154,6 +154,8 @@ section. The complete algorithm is as follows: 6. If :math:`\xi_2 < (1 + \mu^2)/2`, accept :math:`\mu`. Otherwise, repeat the sampling at step 3. +.. _incoherent-sampling: + Incoherent (Compton) Scattering ------------------------------- @@ -166,12 +168,11 @@ the two authors who discovered it: .. math:: :label: klein-nishina - \frac{d\sigma_{KN}}{d\mu} = \pi r_e^2 \left ( \frac{\alpha'}{\alpha} \right - )^2 \left [ \frac{\alpha'}{\alpha} + \frac{\alpha}{\alpha'} + \mu^2 - 1 - \right ] + \frac{d\sigma_{KN}}{d\mu} = \pi r_e^2 \left ( \frac{k'}{k} \right)^2 \left + [ \frac{k'}{k} + \frac{k}{k'} + \mu^2 - 1 \right ] -where :math:`\alpha` and :math:`\alpha'` are the ratios of the incoming and -exiting photon energies to the electron rest mass energy equivalent (0.511 MeV), +where :math:`k` and :math:`k'` are the ratios of the incoming and exiting +photon energies to the electron rest mass energy equivalent (0.511 MeV), respectively. Although it appears that the outgoing energy and angle are separate, there is actually a one-to-one relationship between them such that only one needs to be sampled: @@ -179,32 +180,31 @@ only one needs to be sampled: .. math:: :label: compton-energy-angle - \alpha' = \frac{\alpha}{1 + \alpha(1 - \mu)}. + k' = \frac{k}{1 + k(1 - \mu)}. -Note that when :math:`\alpha'/\alpha` goes to one, i.e., scattering is elastic, -the Klein-Nishina cross section becomes identical to the Thomson cross -section. In general though, the scattering is inelastic and is known as Compton -scattering. When a photon interacts with a bound electron in an atom, the -Klein-Nishina formula must be modified to account for the binding effects. As in -the case of coherent scattering, this is done by means of a form factor. The -differential cross section for incoherent scattering is given by +Note that when :math:`k'/k` goes to one, i.e., scattering is elastic, the +Klein-Nishina cross section becomes identical to the Thomson cross section. In +general though, the scattering is inelastic and is known as Compton scattering. +When a photon interacts with a bound electron in an atom, the Klein-Nishina +formula must be modified to account for the binding effects. As in the case of +coherent scattering, this is done by means of a form factor. The 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_e^2 \left ( - \frac{\alpha'}{\alpha} \right )^2 \left [ \frac{\alpha'}{\alpha} + - \frac{\alpha}{\alpha'} + \mu^2 - 1 \right ] S(x,Z) + \frac{k'}{k} \right )^2 \left [ \frac{k'}{k} + \frac{k}{k'} + \mu^2 - 1 + \right ] S(x,Z) where :math:`S(x,Z)` is the form factor. The approach in OpenMC is to first sample the Klein-Nishina cross section and then perform rejection sampling on the form factor. As in other codes, `Kahn's rejection method`_ is used for -:math:`\alpha < 3` and a direct method by Koblinger_ is used for :math:`\alpha -\ge 3`. The complete algorithm is as follows: +:math:`k < 3` and a direct method by Koblinger_ is used for :math:`k \ge 3`. +The complete algorithm is as follows: -1. If :math:`\alpha < 3`, sample :math:`\mu` from the Klein-Nishina cross - section using Kahn's rejection method. Otherwise, use Koblinger's direct - method. +1. If :math:`k < 3`, sample :math:`\mu` from the Klein-Nishina cross section + using Kahn's rejection method. Otherwise, use Koblinger's direct method. 2. Calculate :math:`x` and :math:`\bar{x}` using :eq:`momentum-transfer` and :eq:`xmax`, respectively. @@ -215,19 +215,370 @@ the form factor. As in other codes, `Kahn's rejection method`_ is used for Doppler Energy Broadening +++++++++++++++++++++++++ -LA-UR-04-0487_ and LA-UR-04-0488_ +Bound electrons are not at rest but have a momentum distribution that will +cause the energy of the scattered photon to be Doppler broadened. More tightly +bound electrons have a wider momentum distribution, so the energy spectrum of +photons scattering off inner shell electrons will be broadened the most. +In addition, scattering from bound electrons places a limit on the maximum +scattered photon energy: + +.. math:: + :label: max-energy-out + + E'_{\text{max}} = E - E_{b,i}, + +where :math:`E_{b,i}` is the binding energy of the :math:`i`-th subshell. + +Compton profiles :math:`J_i(p_z)` are used to account for the binding effects. +The quantity :math:`p_z = {\bf p} \cdot {\bf q}/q` is the projection of the +initial electron momentum on :math:`{\bf q}`, where the scattering vector +:math:`{\bf q} = {\bf p} - {\bf p'}` is the momentum gained by the photon, +:math:`{\bf p}` is the initial momentum of the electron, and :math:`{\bf p'}` +is the momentum of the scattered electron. Applying the conservation of energy +and momentum, :math:`p_z` can be written in terms of the photon energy and +scattering angle: + +.. math:: + :label: pz + + p_z = \frac{E - E' - EE'(1 - \mu)/(m_e c^2)}{-\alpha \sqrt{E^2 + E'^2 - + 2EE'\mu}}, + +where :math:`\alpha` is the fine structure constant. The maximum momentum +transferred, :math:`p_{z,\text{max}}`, can be calculated from :eq:`pz` using +:math:`E' = E'_{\text{max}}`. The Compton profile of the :math:`i`-th electron +subshell is defined as + +.. math:: + :label: compton-profile + + J_i(p_z) = \int \int \rho_i({\bf p}) dp_x dp_y, + +where :math:`\rho_i({\bf p})` is the initial electron momentum distribution. +:math:`J_i(p_z)` can be interpreted as the probability density function of +:math:`p_z`. + +The Doppler broadened energy of the Compton-scattered photon can be sampled by +selecting an electron shell, sampling a value of :math:`p_z` using the Compton +profile, and calculating the scattered photon energy. The theory and methods +used to do this are described in detail in LA-UR-04-0487_ and LA-UR-04-0488_. +The sampling algorithm is summarized below: + +1. Sample :math:`\mu` from :eq:`incoherent-xs` using the algorithm described in + :ref:`incoherent-sampling`. + +2. Sample the electron subshell :math:`i` using the number of electrons per + shell as the probability mass function. + +3. Sample :math:`p_z` using :math:`J_i(p_z)` as the PDF. + +4. Calculate :math:`E'` by solving :eq:`pz` for :math:`E'` using the sampled + value of :math:`p_z`. + +5. If :math:`p_z < p_{z,\text{max}}` for shell :math:`i`, accept :math:`E'`. + Otherwise repeat from step 2. Compton Electrons +++++++++++++++++ +Because the Compton-scattered photons can transfer a large fraction of their +energy to the kinetic energy of the recoil electron, which may in turn go on to +lose its energy as bremsstrahlung radiation, it is necessary to accurately +model the angular and energy distributions of Compton electrons. The energy of +the recoil electron ejected from the :math:`i`-th subshell is given by + +.. math:: + :label: compton-electron-energy + + E_{-} = E - E' - E_{b,i}. + +The direction of the electron is assumed to be in the direction of the momentum +transfer, with the cosine of the polar angle given by + +.. math:: + :label: compton-electron-mu + + \mu_{-} = \frac{E - E'\mu}{\sqrt{E^2 +E'^2 - 2EE'\mu}} + +and the azimuthal angle :math:`\phi_{-} = \phi + \pi`, where :math:`\phi` is +the azimuthal angle of the photon. The vacancy left by the ejected electron is +filled through atomic relaxation. Photoelectric Effect -------------------- +In the photoelectric effect, the incident photon is absorbed by an atomic +electron, which is then emitted from the :math:`i`-th shell with kinetic energy + +.. math:: + :label: photoelectron-kinetic-energy + + E_{-} = E - E_{b,i}. + +Photoelectric emission is only possible when the photon energy exceeds the +binding energy of the shell. These binding energies are often referred to as +edge energies because the otherwise continuously decreasing cross section has +discontinuities at these points, creating the characteristic sawtooth shape. +The photoelectric effect dominates at low energies and is more important for +heavier elements. + +When simulating the photoelectric effect, the first step is to sample the +electron shell. The shell :math:`i` where the ionization occurs can be +considered a discrete random variable with probability mass function + +.. math:: + :label: photoelectron-shell-pdf + + p_i = \frac{\sigma_{\text{pe},i}}{\sigma_{\text{pe}}}, + +where :math:`\sigma_{\text{pe},i}` is the cross section of the :math:`i`-th +shell, and the total photoelectric cross section of the atom, +:math:`\sigma_{\text{pe}}`, is the sum over the shell cross sections. Once the +shell has been sampled, the energy of the photoelectron is calculated using +:eq:`photoelectron-kinetic-energy`. + +To determine the direction of the photoelectron, we implement the method +described in Kaltiaisenaho_, which models the angular distribution of the +photoelectrons using the K-shell cross section derived by Sauter (K-shell +electrons are the most tightly bound, and they contribute the most to +:math:`\sigma_{\text{pe}}`). The non-relativistic Sauter distribution for +unpolarized photons can be approximated as + +.. math:: + :label: sauter + + \frac{d\sigma_{\text{pe}}}{d\mu_{-}} \propto + \frac{1 - \mu_{-}^2}{(1 - \beta_{-} \mu_{-})^4}, + +where :math:`\beta_{-}` is the ratio of the velocity of the electron to the +speed of light, + +.. math:: + :label: beta-2 + + \beta_{-} = \frac{\sqrt{(E_{-}(E_{-} + 2m_e c^2)}}{E_{-} + m_e c^2}. + +To sample :math:`\mu_{-}` from the Sauter distribution, we first express +:eq:`sauter` in the form: + +.. math:: + :label: photoelectron-mu-pdf + + f(\mu_{-}) = \frac{3}{2} \psi(\mu_{-}) g(\mu_{-}), + +where + +.. math:: + :label: mu-pdf-factors + + \psi(\mu_{-}) &= \frac{(1 - \beta_{-}^2)(1 - \mu_{-}^2)}{(1 - + \beta_{-}\mu_{-})^2}, \\ + g(\mu_{-}) &= \frac{1 - \beta_{-}^2}{2 (1 - \beta_{-}\mu_{-})^2}. + +In the interval :math:`[-1, 1]`, :math:`g(\mu_{-})` is a normalized PDF and +:math:`\psi(\mu_{-})` satisfies the condition :math:`0 < \psi(\mu_{-}) < 1`. +The following algorithm can now be used to sample :math:`\mu_{-}`: + +1. Using the inverse transform method, sample :math:`\mu_{-}` from + :math:`g(\mu_{-})` using the sampling formula + + .. math:: + + \mu_{-} = \frac{2\xi_1 + \beta_{-} - 1}{2\beta_{-}\xi_1 - \beta_{-} + 1}. + +2. If :math:`\xi_2 \le \psi(\mu_{-})`, accept :math:`\mu_{-}`. Otherwise, + repeat the sampling from step 1. + +The azimuthal angle is sampled uniformly on :math:`[0, 2\pi)`. + +The atom is left in an excited state with a vacancy in the :math:`i`-th shell +and decays to its ground state through a cascade of transitions that produce +fluorescent photons and Auger electrons. Pair Production --------------- +In electron-positron pair production, a photon is absorbed in the vicinity of +an atomic nucleus or an electron and an electron and positron are created. Pair +production is the dominant interaction with matter at high photon energies and +is more important for high-Z elements. When it takes place in the field of a +nucleus, energy is essentially conserved among the incident photon and the +resulting charged particles. Therefore, in order for pair production to occur, +the photon energy must be greater than the sum of the rest mass energies of the +electron and positron, i.e., :math:`E_{\text{threshold,pp}} = 2 m_e c^2 = +1.022` MeV. + +The photon can also interact in the field of an atomic electron. This process +is referred to as "triplet production" because the target electron is ejected +from the atom and three charged particles emerge from the interaction. In this +case, the recoiling electron also absorbs some energy, so the energy threshold +for triplet production is greater than that of pair production from atomic +nuclei, with :math:`E_{\text{threshold,tp}} = 4 m_e c^2 = 2.044` MeV. The ratio +of the triplet production cross section to the pair production cross section is +approximately 1/Z, so triplet production becomes increasingly unimportant for +high-Z elements. Though it can be significant in lighter elements, the momentum +of the recoil electron becomes negligible in the energy regime where pair +production dominates. For our purposes, it is a good approximation to treat +triplet production as pair production and only simulate the electron-positron +pair. + +Accurately modeling the creation of electron-positron pair is important because +the charged particles can go on to lose much of their energy as bremsstrahlung +radiation, and the subsequent annihilation of the positron with an electron +produces two additional photons. We sample the energy and direction of the +charged particles using a semiempirical model described in Salvat_. The +Bethe-Heitler differential cross section, given by + +.. math:: + :label: bethe-heitler + + \frac{d\sigma_{\text{pp}}}{d\epsilon} = \alpha r_e^2 Z^2 + \left[ (\epsilon^2 + (1-\epsilon)^2) (\Phi_1 - 4f_C) + + \frac{2}{3}\epsilon(1-\epsilon)(\Phi_2 - 4f_C) \right], + +is used as a starting point, where :math:`\alpha` is the fine structure +constant, :math:`f_C` is the Coulomb correction function, :math:`\Phi_1` and +:math:`\Phi_2` are screening functions, and :math:`\epsilon = (E_{-} + m_e +c^2)/E` is the electron reduced energy (i.e., the fraction of the photon energy +given to the electron). :math:`\epsilon` can take values between +:math:`\epsilon_{\text{min}} = k^{-1}` (when the kinetic energy of the electron +is zero) and :math:`\epsilon_{\text{max}} = 1 - k^{-1}` (when the kinetic +energy of the positron is zero). + +The Coulomb correction, given by + +.. math:: + :label: coulomb-correction + + f_C = \alpha^{2}Z^{2} \big[&(1 + \alpha^{2}Z^{2})^{-1} + 0.202059 + - 0.03693\alpha^{2}Z^{2} + 0.00835\alpha^{4}Z^{4} \\ + &- 0.00201\alpha^{6}Z^{6} + 0.00049\alpha^{8}Z^{8} + - 0.00012\alpha^{10}Z^{10} + 0.00003\alpha^{12}Z^{12}\big] + +is introduced to correct for the fact that the Bethe-Heitler differential cross +section was derived using the Born approximation, which treats the Coulomb +interaction as a small perturbation. + +The screening functions :math:`\Phi_1` and :math:`\Phi_2` account for the +screening of the Coulomb field of the atomic nucleus by outer electrons. Since +they are given by integrals which include the atomic form factor, they must be +computed numerically for a realistic form factor. However, by assuming +exponential screening and using a simplified form factor, analytical +approximations of the screening functions can be derived: + +.. math:: + :label: screening-functions + + \Phi_1 &= 2 - 2\ln(1 + b^2) - 4b\arctan(b^{-1}) + 4\ln(Rm_{e}c/\hbar) \\ + \Phi_2 &= \frac{4}{3} - 2\ln(1 + b^2) + 2b^2 \left[ 4 - 4b\arctan(b^{-1}) + - 3\ln(1 + b^{-2}) \right] + 4\ln(Rm_{e}c/\hbar) + +where + +.. math:: + :label: b + + b = \frac{Rm_{e}c}{2k\epsilon(1 - \epsilon)\hbar}. + +and :math:`R` is the screening radius. + +The differential cross section in :eq:`bethe-heitler` with the approximations +described above will not be accurate at low energies: the lower boundary of +:math:`\epsilon` will be shifted above :math:`\epsilon_{\text{min}}` and the +upper boundary of :math:`\epsilon` will be shifted below +:math:`\epsilon_{\text{max}}`. To offset this behavior, a correcting factor +:math:`F_0(k, Z)` is used: + +.. math:: + :label: correcting-factor + + F_0(k, Z) =~& (0.1774 + 12.10\alpha Z - 11.18\alpha^{2}Z^{2})(2/k)^{1/2} \\ + &+ (8.523 + 73.26\alpha Z - 44.41\alpha^{2}Z^{2})(2/k) \\ + &- (13.52 + 121.1\alpha Z - 96.41\alpha^{2}Z^{2})(2/k)^{3/2} \\ + &+ (8.946 + 62.05\alpha Z - 63.41\alpha^{2}Z^{2})(2/k)^{2}. + +To aid sampling, the differential cross section used to sample :math:`\epsilon` +(minus the normalization constant) can now be expressed in the form + +.. math:: + :label: pp-pdf + + \frac{d\sigma_{\text{pp}}}{d\epsilon} = + u_1 \frac{\phi_1(\epsilon)}{\phi_1(1/2)} \pi_1(\epsilon) + + u_2 \frac{\phi_2(\epsilon)}{\phi_2(1/2)} \pi_2(\epsilon) + +where + +.. math:: + :label: u + + u_1 &= \frac{2}{3} \left(\frac{1}{2} - \frac{1}{k}\right)^2 \phi_1(1/2), \\ + u_2 &= \phi_2(1/2), + +.. math:: + :label: phi + + \phi_1(\epsilon) &= \frac{1}{2}(3\Phi_1 - \Phi_2) - 4f_{C}(Z) + F_0(k, Z), \\ + \phi_2(\epsilon) &= \frac{1}{4}(3\Phi_1 + \Phi_2) - 4f_{C}(Z) + F_0(k, Z), + +and + +.. math:: + :label: pi + + \pi_1(\epsilon) &= \frac{3}{2} \left(\frac{1}{2} - \frac{1}{k}\right)^{-3} + \left(\frac{1}{2} - \epsilon\right)^2, \\ + \pi_2(\epsilon) &= \frac{1}{2} \left(\frac{1}{2} - \frac{1}{k}\right)^{-1}. + +The functions in :eq:`phi` are non-negative and maximum at :math:`\epsilon = +1/2`. In the interval :math:`(\epsilon_{\text{min}}, \epsilon_{\text{max}})`, +the functions in :eq:`pi` are normalized PDFs and +:math:`\phi_i(\epsilon)/\phi_i(1/2)` satisfies the condition :math:`0 < +\phi_i(\epsilon)/\phi_i(1/2) < 1`. The following algorithm can now be used to +sample the reduced electron energy :math:`\epsilon`: + +1. Sample :math:`i` according to the point probabilities + :math:`p(i=1) = u_1/(u_1 + u_2)` and :math:`p(i=2) = u_2/(u_1 + u_2)`. + +2. Using the inverse transform method, sample :math:`\epsilon` from + :math:`\pi_i(\epsilon)` using the sampling formula + + .. math:: + + \epsilon &= \frac{1}{2} + \left(\frac{1}{2} - \frac{1}{k}\right) + (2\xi_1 - 1)^{1/3} ~~~~&\text{if}~~ i = 1 \\ + \epsilon &= \frac{1}{k} + \left(\frac{1}{2} - + \frac{1}{k}\right) 2\xi_1 ~~~~&\text{if}~~ i = 2. + +3. If :math:`\xi_2 \le \phi_i(\epsilon)/\phi_i(1/2)`, accept + :math:`\epsilon`. Otherwise, repeat the sampling from step 1. + +Because charged particles have a much smaller range than the mean free path of +photons and because they immediately undergo multiple scattering events which +randomize their direction, it is sufficient to use a simplified model to sample +the direction of the electron and positron. The cosines of the polar angles are +sampled using the leading order term of the Sauter–Gluckstern–Hull +distribution, + +.. math:: + :label: sauter–gluckstern–hull + + p(\mu_{\pm}) = C(1 - \beta_{\pm}\mu_{\pm})^{-2}, + +where :math:`C` is a normalization constant and :math:`\beta_{\pm}` is the +ratio of the velocity of the charged particle to the speed of light given in +:eq:`beta-2`. + +The inverse transform method is used to sample :math:`\mu_{-}` and +:math:`\mu_{+}` from :eq:`sauter–gluckstern–hull`, using the sampling formula + +.. math:: + :label: sample-mu + + \mu_{\pm} = \frac{2\xi - 1 + \beta_{\pm}}{(2\xi - 1)\beta_{\pm} + 1}. + +The azimuthal angles for the electron and positron are sampled independently +and uniformly on :math:`[0, 2\pi)`. ------------------- Secondary Processes @@ -248,19 +599,297 @@ additional photons. Atomic Relaxation ----------------- +When an electron is ejected from an atom and a vacancy is left in an inner +shell, an electron from a higher energy level will fill the vacancy. This +results in either a radiative transition, in which a photon with a +characteristic energy (fluorescence photon) is emitted, or non-radiative +transition, in which an electron from a shell that is farther out (Auger +electron) is emitted. If a non-radiative transition occurs, the new vacancy is +filled in the same manner, and as the process repeats a shower of photons and +electrons can be produced. + +The energy of a fluorescence photon is the equal to the energy difference +between the transition states, i.e., + +.. math:: + :label: fluorescence-photon-energy + + E = E_{b,v} - E_{b,i}, + +where :math:`E_{b,v}` is the binding energy of the vacancy shell and +:math:`E_{b,i}` is the binding energy of the shell from which the electron +transitioned. The energy of an Auger electron is given by + +.. math:: + :label: auger-electron-energy + + E_{-} = E_{b,v} - E_{b,i} - E_{b,a}, + +where :math:`E_{b,a}` is the binding energy of the shell from which the Auger +electron is emitted. While Auger electrons are low-energy so their range and +bremsstrahlung yield is small, fluorescence photons can travel far before +depositing their energy, so the relaxation process should be modeled in detail. + +Transition energies and probabilities are needed for each subshell to simulate +atomic relaxation. Starting with the initial shell vacancy, the following +recursive algorithm is used to fill vacancies and create fluorescence photons +and Auger electrons: + +1. If there are no transitions for the vacancy shell, create a fluorescence + photon assuming it is from a captured free electron and terminate. + +2. Sample a transition using the transition probabilities for the vacancy + shell as the probability mass function. + +3. Create either a fluorescence photon or Auger electron, sampling the + direction of the particle isotropically. + +4. If a non-radiative transition occurred, repeat from step 1 for the vacancy + left by the emitted Auger electron. + +5. Repeat from step 1 for vacancy left by the transition electron. Electron-Positron Annihilation ------------------------------ +When a positron collides with an electron, both particles are annihilated and +generally two photons with equal energy are created. If the kinetic energy of +the positron is high enough, the two photons can have different energies, and +the higher-energy photon is emitted preferentially in the direction of flight +of the positron. It is also possible to produce a single photon if the +interaction occurs with a bound electron, and in some cases three (or, rarely, +even more) photons can be emitted. However, the annihilation cross section is +largest for low-energy positrons, and as the positron energy decreases, the +angular distribution of the emitted photons becomes isotropic. + +In OpenMC, we assume the most likely case in which a low-energy positron (which +has already lost most of its energy to bremsstrahlung radiation) interacts with +an electron which is free and at rest. Two photons with energy equal to the +electron rest mass energy :math:`m_e c^2 = 0.511` MeV are emitted isotropically +in opposite directions. Bremsstrahlung -------------- +When a charged particle is decelerated in the field of an atom, some of its +kinetic energy is converted into electromagnetic radiation known as +bremsstrahlung, or 'braking radiation'. In each event, an electron or positron +with kinetic energy :math:`T` generates a photon with an energy :math:`E` +between :math:`0` and :math:`T`. Bremsstrahlung is described by a cross section +that is differential in photon energy, in the direction of the emitted photon, +and in the final direction of the charged particle. However, in Monte Carlo +simulations it is typical to integrate over the angular variables to obtain a +single differential cross section with respect to photon energy, which is often +expressed in the form + +.. math:: + :label: bremsstrahlung-dcs + + \frac{d\sigma_{\text{br}}}{dE} = \frac{Z^2}{\beta^2} \frac{1}{E} + \chi(Z, T, \kappa), + +where :math:`\kappa = E/T` is the reduced photon energy and :math:`\chi(Z, T, +\kappa)` is the scaled bremsstrahlung cross section, which is experimentally +measured. + +Because electrons are attracted to atomic nuclei whereas positrons are +repulsed, the cross section for positrons is smaller, though it approaches that +of electrons in the high energy limit. To obtain the positron cross section, we +multiply :eq:`bremsstrahlung-dcs` by the :math:`\kappa`-independent factor used +in Salvat_, + +.. math:: + :label: positron-factor + + F_{\text{p}}(Z,T) = + & 1 - \text{exp}(-1.2359\times 10^{-1}t + 6.1274\times 10^{-2}t^2 - 3.1516\times 10^{-2}t^3 \\ + & + 7.7446\times 10^{-3}t^4 - 1.0595\times 10^{-3}t^5 + 7.0568\times 10^{-5}t^6 \\ + & - 1.8080\times 10^{-6}t^7), + +where + +.. math:: + :label: positron-factor-t + + t = \ln\left(1 + \frac{10^6}{Z^2}\frac{T}{\text{m}_\text{e}c^2} \right). + +:math:`F_{\text{p}}(Z,T)` is the ratio of the radiative stopping powers for +positrons and electrons. Stopping power describes the average energy loss per +unit path length of a charged particle as it passes through matter: + +.. math:: + :label: stopping-power + + -\frac{dT}{ds} = n \int E \frac{d\sigma}{dE} dE \equiv S(T), + +where :math:`n` is the number density of the material and :math:`d\sigma/dE` is +the cross section differential in energy loss. The total stopping power +:math:`S(T)` can be separated into two components: the radiative stopping +power :math:`S_{\text{rad}}(T)`, which refers to energy loss due to +bremsstrahlung, and the collision stopping power :math:`S_{\text{col}}(T)`, +which refers to the energy loss due to inelastic collisions with bound +electrons in the material that result in ionization and excitation. To obtain +the radiative stopping power for positrons, the radiative stopping power for +electrons is multiplied by :eq:`positron-factor`. Currently, the collision +stopping power for electrons is also used for positrons. + +While the models for photon interactions with matter described above can safely +assume interactions occur with free atoms, sampling the target atom based on +the macroscopic cross sections, molecular effects cannot necessarily be +disregarded for charged particle treatment. For compounds and mixtures, the +bremsstrahlung cross section is calculated using Bragg's additivity rule as + +.. math:: + :label: material-bremsstrahlung-dcs + + \frac{d\sigma_{\text{br}}}{dE} = \frac{1}{\beta^2 E} \sum_i \gamma_i Z^2_i + \chi(Z_i, T, \kappa), + +where the sum is over the constituent elements and :math:`\gamma_i` is the +atomic fraction of the :math:`i`-th element. Similarly, the radiative stopping +power is calculated using Bragg's additivity rule as + +.. math:: + :label: material-radiative-stopping-power + + S_{\text{rad}}(T) = \sum_i w_i S_{\text{rad},i}(T), + +where :math:`w_i` is the mass fraction of the :math:`i`-th element. The +collision stopping power, however, is a function of certain quantities such as +the mean excitation energy :math:`I` and the density effect correction +:math:`\delta_F` that depend on molecular properties. These quantities cannot +simply be summed over constituent elements in a compound, but should instead be +calculated for the material. Currently, we use Bragg's additivity rule to +calculate the collision stopping power as well, but this is not a good +approximation and should be fixed in the future. + .. _ttb: Thick-Target Bremsstrahlung Approximation +++++++++++++++++++++++++++++++++++++++++ +Since charged particles lose their energy on a much shorter distance scale than +neutral particles, not much error should be introduced by neglecting to +transport electrons. However, the bremsstrahlung emitted from high energy +electrons and positrons can travel far from the interaction site. Thus, even +without a full electron transport mode it is necessary to model bremsstrahlung. +We use a thick-target bremsstrahlung (TTB) approximation based on the models in +Salvat_ and Kaltiaisenaho_ for generating bremsstrahlung photons, which assumes +the charged particle loses all its energy in a single homogeneous material +region. + +To model bremsstrahlung using the TTB approximation, we need to know the number +of photons emitted by the charged particle and the energy distribution of the +photons. These quantities can be calculated using the continuous slowing down +approximation (CSDA). The CSDA assumes charged particles lose energy +continuously along their trajectory with a rate of energy loss equal to the +total stopping power, ignoring fluctuations in the energy loss. The +approximation is useful for expressing average quantities that describe how +charged particles slow down in matter. For example, the CSDA range approximates +the average path length a charged particle travels as it slows to rest: + +.. math:: + :label: csda-range + + R(T) = \int^T_0 \frac{dT'}{S(T')}. + +Actual path lengths will fluctuate around :math:`R(T)`. The average number of +photons emitted per unit path length is given by the inverse bremsstrahlung +mean free path: + +.. math:: + :label: inverse-bremsstrahlung-mfp + + \lambda_{\text{br}}^{-1}(T,E_{\text{cut}}) + = n\int_{E_{\text{cut}}}^T\frac{d\sigma_{\text{br}}}{dE}dE + = n\frac{Z^2}{\beta^2}\int_{\kappa_{\text{cut}}}^1\frac{1}{\kappa} + \chi(Z,T,\kappa)d\kappa. + +The lower limit of the integral in :eq:`inverse-bremsstrahlung-mfp` is non-zero +because the bremsstrahlung differential cross section diverges for small photon +energies but is finite for photon energies above some cutoff energy +:math:`E_{\text{cut}}`. The mean free path +:math:`\lambda_{\text{br}}^{-1}(T,E_{\text{cut}})` is used to calculate the +photon number yield, defined as the average number of photons emitted with +energy greater than :math:`E_{\text{cut}}` as the charged particle slows down +from energy :math:`T` to :math:`E_{\text{cut}}`. The photon number yield is +given by + +.. math:: + :label: photon-number-yield + + Y(T,E_{\text{cut}}) = \int^{R(T)}_{R(E_{\text{cut}})} + \lambda_{\text{br}}^{-1}(T',E_{\text{cut}})ds = \int_{E_{\text{cut}}}^T + \frac{\lambda_{\text{br}}^{-1}(T',E_{\text{cut}})}{S(T')}dT'. + +:math:`Y(T,E_{\text{cut}})` can be used to construct the energy spectrum of +bremsstrahlung photons: the number of photons created with energy between +:math:`E_1` and :math:`E_2` by a charged particle with initial kinetic energy +:math:`T` as it comes to rest is given by :math:`Y(T,E_1) - Y(T,E_2)`. + +To simulate the emission of bremsstrahlung photons, the total stopping power +and bremsstrahlung differential cross section for positrons and electrons must +be calculated for a given material using :eq:`material-bremsstrahlung-dcs` and +:eq:`material-radiative-stopping-power`. These quantities are used to build the +tabulated bremsstrahlung energy PDF and CDF for that material for each incident +energy :math:`T_k` on the energy grid. The following algorithm is then applied +to sample the photon energies: + +1. For an incident charged particle with energy :math:`T`, sample the number of + emitted photons as + + .. math:: + + N = \lfloor Y(T,E_{\text{cut}}) + \xi_1 \rfloor. + +2. Rather than interpolate the PDF between indices :math:`k` and :math:`k+1` + for which :math:`T_k < T < T_{k+1}`, which is computationally expensive, use + the composition method and sample from the PDF at either :math:`k` or + :math:`k+1`. Using linear interpolation on a logarithmic scale, the PDF can + be expressed as + + .. math:: + + p_{\text{br}}(T,E) = \pi_k p_{\text{br}}(T_k,E) + \pi_{k+1} + p_{\text{br}}(T_{k+1},E), + + where the interpolation weights are + + .. math:: + + \pi_k = \frac{\ln T_{k+1} - \ln T}{\ln T_{k+1} - \ln T_k},~~~ + \pi_{k+1} = \frac{\ln T - \ln T_k}{\ln T_{k+1} - \ln T_k}. + + Sample either the index :math:`i = k` or :math:`i = k+1` according to the + point probabilities :math:`\pi_{k}` and :math:`\pi_{k+1}`. + +3. Determine the maximum value of the CDF :math:`P_{\text{br,max}}`. + +3. Sample the photon energies using the inverse transform method with the + tabulated CDF :math:`P_{\text{br}}(T_i, E)` i.e., + + .. math:: + + E = E_j \left[ (1 + a_j) \frac{\xi_2 P_{\text{br,max}} - + P_{\text{br}}(T_i, E_j)} {E_j p_{\text{br}}(T_i, E_j)} + 1 + \right]^{\frac{1}{1 + a_j}} + + where the interpolation factor :math:`a_j` is given by + + .. math:: + + a_j = \frac{\ln p_{\text{br}}(T_i,E_{j+1}) - \ln p_{\text{br}}(T_i,E_j)} + {\ln E_{j+1} - \ln E_j} + + and :math:`P_{\text{br}}(T_i, E_j) \le \xi_2 P_{\text{br,max}} \le + P_{\text{br}}(T_i, E_{j+1})`. + +We ignore the range of the electron or positron, i.e., the bremsstrahlung +photons are produced in the same location that the charged particle was +created. The direction of the photons is assumed to be the same as the +direction of the incident charged particle, which is a reasonable approximation +at higher energies when the bremsstrahlung radiation is emitted at small +angles. .. _Koblinger: https://doi.org/10.13182/NSE75-A26663 @@ -273,3 +902,7 @@ Thick-Target Bremsstrahlung Approximation .. _LA-UR-04-0487: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-04-0487.pdf .. _LA-UR-04-0488: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-04-0488.pdf + +.. _Kaltiaisenaho: https://aaltodoc.aalto.fi/bitstream/handle/123456789/21004/master_Kaltiaisenaho_Toni_2016.pdf + +.. _Salvat: http://www.oecd-nea.org/globalsearch/download.php?doc=77434 diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index 468c25980c..fa64605513 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -1,9 +1,11 @@ #ifndef OPENMC_HDF5_INTERFACE_H #define OPENMC_HDF5_INTERFACE_H +#include // for min #include #include #include +#include // for strlen #include #include #include @@ -162,13 +164,13 @@ void read_attribute(hid_t obj_id, const char* name, xt::xarray& arr) std::size_t size = 1; for (const auto x : shape) size *= x; - T* buffer = new T[size]; + std::vector buffer(size); // Read data from attribute - read_attr(obj_id, name, H5TypeMap::type_id, buffer); + read_attr(obj_id, name, H5TypeMap::type_id, buffer.data()); // Adapt array into xarray - arr = xt::adapt(buffer, size, xt::acquire_ownership(), shape); + arr = xt::adapt(buffer, shape); } // overload for std::string @@ -193,13 +195,13 @@ read_attribute(hid_t obj_id, const char* name, std::vector& vec) // Allocate a C char array to get strings auto n = attribute_typesize(obj_id, name); - char buffer[m][n+1]; + char buffer[m][n]; // Read char data in attribute read_attr_string(obj_id, name, n, buffer[0]); for (int i = 0; i < m; ++i) { - vec.emplace_back(&buffer[i][0]); + vec.emplace_back(&buffer[i][0], std::min(strlen(buffer[i]), n)); } } @@ -244,13 +246,13 @@ void read_dataset(hid_t dset, xt::xarray& arr, bool indep=false) std::size_t size = 1; for (const auto x : shape) size *= x; - T* buffer = new T[size]; + std::vector buffer(size); // Read data from attribute - read_dataset(dset, nullptr, H5TypeMap::type_id, buffer, indep); + read_dataset(dset, nullptr, H5TypeMap::type_id, buffer.data(), indep); // Adapt into xarray - arr = xt::adapt(buffer, size, xt::acquire_ownership(), shape); + arr = xt::adapt(buffer, shape); } template @@ -273,13 +275,13 @@ void read_dataset_as_shape(hid_t obj_id, const char* name, std::size_t size = 1; for (const auto x : arr.shape()) size *= x; - T* buffer = new T[size]; + std::vector buffer(size); // Read data from attribute - read_dataset(dset, nullptr, H5TypeMap::type_id, buffer, indep); + read_dataset(dset, nullptr, H5TypeMap::type_id, buffer.data(), indep); // Adapt into xarray - arr = xt::adapt(buffer, size, xt::acquire_ownership(), arr.shape()); + arr = xt::adapt(buffer, arr.shape()); close_dataset(dset); } diff --git a/include/openmc/search.h b/include/openmc/search.h index 81370a3de0..c91ad18fe1 100644 --- a/include/openmc/search.h +++ b/include/openmc/search.h @@ -14,6 +14,7 @@ template typename std::iterator_traits::difference_type lower_bound_index(It first, It last, const T& value) { + if (*first == value) return 0; It index = std::lower_bound(first, last, value) - 1; return (index == last) ? -1 : index - first; } diff --git a/openmc/capi/core.py b/openmc/capi/core.py index 53c3c32bd6..fa6517488d 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -62,6 +62,18 @@ def calculate_volumes(): _dll.openmc_calculate_volumes() +def current_batch(): + """Return the current batch of the simulation. + + Returns + ------- + int + Current batch of the simulation + + """ + return c_int.in_dll(_dll, 'openmc_current_batch').value + + def finalize(): """Finalize simulation and free memory""" _dll.openmc_finalize() @@ -214,6 +226,18 @@ def keff(): return (mean, std_dev) +def master(): + """Return whether processor is master processor or not. + + Returns + ------- + bool + Whether is master processor or not + + """ + return c_bool.in_dll(_dll, 'openmc_master').value + + def next_batch(): """Run next batch. diff --git a/openmc/capi/settings.py b/openmc/capi/settings.py index 1063d6463e..115962395d 100644 --- a/openmc/capi/settings.py +++ b/openmc/capi/settings.py @@ -1,4 +1,5 @@ -from ctypes import c_int, c_int32, c_int64, c_double, c_char_p, POINTER +from ctypes import (c_int, c_int32, c_int64, c_double, c_char_p, c_bool, + POINTER) from . import _dll from .core import _DLLGlobal @@ -17,9 +18,11 @@ _dll.openmc_get_seed.restype = c_int64 class _Settings(object): # Attributes that are accessed through a descriptor batches = _DLLGlobal(c_int32, 'n_batches') + entropy_on = _DLLGlobal(c_bool, 'entropy_on') generations_per_batch = _DLLGlobal(c_int32, 'gen_per_batch') inactive = _DLLGlobal(c_int32, 'n_inactive') particles = _DLLGlobal(c_int64, 'n_particles') + run_CE = _DLLGlobal(c_bool, 'run_CE') verbosity = _DLLGlobal(c_int, 'verbosity') @property diff --git a/scripts/openmc-voxel-to-silovtk b/scripts/openmc-voxel-to-silovtk deleted file mode 100755 index e0cfc0a5ba..0000000000 --- a/scripts/openmc-voxel-to-silovtk +++ /dev/null @@ -1,85 +0,0 @@ -#!/usr/bin/env python3 - -import struct -import sys -from argparse import ArgumentParser - -import numpy as np -import h5py - - -def main(): - # Process command line arguments - parser = ArgumentParser() - parser.add_argument('voxel_file', help='Path to voxel file') - parser.add_argument('-o', '--output', action='store', - default='plot', help='Path to output SILO or VTK file.') - parser.add_argument('-s', '--silo', action='store_true', - default=False, help='Flag to convert to SILO instead of VTK.') - args = parser.parse_args() - - # Read data from voxel file - fh = h5py.File(args.voxel_file, 'r') - dimension = fh.attrs['num_voxels'] - width = fh.attrs['voxel_width'] - lower_left = fh.attrs['lower_left'] - voxel_data = fh['data'].value - - nx, ny, nz = dimension - upper_right = lower_left + width*dimension - - if not args.silo: - import vtk - - grid = vtk.vtkImageData() - grid.SetDimensions(nx+1, ny+1, nz+1) - grid.SetOrigin(*lower_left) - grid.SetSpacing(*width) - - data = vtk.vtkDoubleArray() - data.SetName("id") - data.SetNumberOfTuples(nx*ny*nz) - for x in range(nx): - sys.stdout.write(" {}%\r".format(int(x/nx*100))) - sys.stdout.flush() - for y in range(ny): - for z in range(nz): - i = z*nx*ny + y*nx + x - data.SetValue(i, voxel_data[x, y, z]) - grid.GetCellData().AddArray(data) - - writer = vtk.vtkXMLImageDataWriter() - if vtk.vtkVersion.GetVTKMajorVersion() > 5: - writer.SetInputData(grid) - else: - writer.SetInput(grid) - if not args.output.endswith(".vti"): - args.output += ".vti" - writer.SetFileName(args.output) - writer.Write() - - else: - import silomesh - - if not args.output.endswith(".silo"): - args.output += ".silo" - silomesh.init_silo(args.output) - meshparams = list(map(int, dimension)) + list(map(float, lower_left)) + \ - list(map(float, upper_right)) - silomesh.init_mesh('plot', *meshparams) - silomesh.init_var("id") - for x in range(nx): - sys.stdout.write(" {}%\r".format(int(x/nx*100))) - sys.stdout.flush() - for y in range(ny): - for z in range(nz): - silomesh.set_value(float(voxel_data[x, y, z]), - x + 1, y + 1, z + 1) - print() - silomesh.finalize_var() - silomesh.finalize_mesh() - silomesh.finalize_silo() - - -if __name__ == '__main__': - main() diff --git a/scripts/openmc-voxel-to-vtk b/scripts/openmc-voxel-to-vtk new file mode 100755 index 0000000000..c0d8e9a147 --- /dev/null +++ b/scripts/openmc-voxel-to-vtk @@ -0,0 +1,58 @@ +#!/usr/bin/env python3 + +import struct +import sys +from argparse import ArgumentParser + +import numpy as np +import h5py +import vtk + + +def main(): + # Process command line arguments + parser = ArgumentParser() + parser.add_argument('voxel_file', help='Path to voxel file') + parser.add_argument('-o', '--output', action='store', + default='plot', help='Path to output VTK file.') + args = parser.parse_args() + + # Read data from voxel file + fh = h5py.File(args.voxel_file, 'r') + dimension = fh.attrs['num_voxels'] + width = fh.attrs['voxel_width'] + lower_left = fh.attrs['lower_left'] + voxel_data = fh['data'].value + + nx, ny, nz = dimension + upper_right = lower_left + width*dimension + + grid = vtk.vtkImageData() + grid.SetDimensions(nx+1, ny+1, nz+1) + grid.SetOrigin(*lower_left) + grid.SetSpacing(*width) + + data = vtk.vtkDoubleArray() + data.SetName("id") + data.SetNumberOfTuples(nx*ny*nz) + for x in range(nx): + sys.stdout.write(" {}%\r".format(int(x/nx*100))) + sys.stdout.flush() + for y in range(ny): + for z in range(nz): + i = z*nx*ny + y*nx + x + data.SetValue(i, voxel_data[x, y, z]) + grid.GetCellData().AddArray(data) + + writer = vtk.vtkXMLImageDataWriter() + if vtk.vtkVersion.GetVTKMajorVersion() > 5: + writer.SetInputData(grid) + else: + writer.SetInput(grid) + if not args.output.endswith(".vti"): + args.output += ".vti" + writer.SetFileName(args.output) + writer.Write() + +if __name__ == '__main__': + main() diff --git a/setup.py b/setup.py index cbf564480f..c0168d8805 100755 --- a/setup.py +++ b/setup.py @@ -64,7 +64,7 @@ kwargs = { # Optional dependencies 'extras_require': { 'test': ['pytest', 'pytest-cov'], - 'vtk': ['vtk', 'silomesh'], + 'vtk': ['vtk'], }, } diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index a12ca217e1..b090c6394c 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -273,39 +273,39 @@ contains ! Left surface cmfd % current(1,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + OUT_LEFT) + score_index + 1 + ng*(OUT_LEFT - 1)) cmfd % current(2,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + IN_LEFT) + score_index + 1 + ng*(IN_LEFT - 1)) ! Right surface cmfd % current(3,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + IN_RIGHT) + score_index + 1 + ng*(IN_RIGHT - 1)) cmfd % current(4,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + OUT_RIGHT) + score_index + 1 + ng*(OUT_RIGHT - 1)) ! Back surface cmfd % current(5,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + OUT_BACK) + score_index + 1 + ng*(OUT_BACK - 1)) cmfd % current(6,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + IN_BACK) + score_index + 1 + ng*(IN_BACK - 1)) ! Front surface cmfd % current(7,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + IN_FRONT) + score_index + 1 + ng*(IN_FRONT - 1)) cmfd % current(8,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + OUT_FRONT) + score_index + 1 + ng*(OUT_FRONT - 1)) ! Left surface cmfd % current(9,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + OUT_BOTTOM) + score_index + 1 + ng*(OUT_BOTTOM - 1)) cmfd % current(10,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + IN_BOTTOM) + score_index + 1 + ng*(IN_BOTTOM - 1)) ! Right surface cmfd % current(11,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + IN_TOP) + score_index + 1 + ng*(IN_TOP - 1)) cmfd % current(12,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + OUT_TOP) + score_index + 1 + ng*(OUT_TOP - 1)) else if (ital == 4) then diff --git a/src/cmfd_execute.cpp b/src/cmfd_execute.cpp index 7774cfe1d2..740c1f51f0 100644 --- a/src/cmfd_execute.cpp +++ b/src/cmfd_execute.cpp @@ -25,8 +25,6 @@ cmfd_populate_sourcecounts(int n_energy, const double* energies, auto& m = meshes.at(index_cmfd_mesh); xt::xarray counts = m->count_sites(openmc_work, source_bank, n_energy, energies, outside); - std::cout << counts << "\n"; - // Copy data from the xarray into the source counts array std::copy(counts.begin(), counts.end(), source_counts); } diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90 index 99d482480a..01404d9b0e 100644 --- a/src/cmfd_solver.F90 +++ b/src/cmfd_solver.F90 @@ -146,8 +146,13 @@ contains call loss % assemble() call prod % assemble() if (cmfd_write_matrices) then - call loss % write('loss.dat') - call prod % write('prod.dat') + if (adjoint) then + call loss % write('adj_loss.dat') + call prod % write('adj_prod.dat') + else + call loss % write('loss.dat') + call prod % write('prod.dat') + end if end if ! Set norms to 0 @@ -740,7 +745,7 @@ contains else filename = 'fluxvec.dat' end if - ! TODO: call phi_n % write(filename) + call phi_n % write(filename) end if end subroutine extract_results diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index 40f56d88ef..8872a84f6a 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -675,6 +675,7 @@ contains ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then + allocate(name_(len_trim(name) + 1)) name_ = to_c_string(name) call read_double_c(obj_id, c_loc(name_), buffer_, indep_) else @@ -698,6 +699,7 @@ contains ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then + allocate(name_(len_trim(name) + 1)) name_ = to_c_string(name) call read_double_c(obj_id, c_loc(name_), buffer, indep_) else @@ -720,6 +722,7 @@ contains ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then + allocate(name_(len_trim(name) + 1)) name_ = to_c_string(name) call read_double_c(obj_id, c_loc(name_), buffer, indep_) else @@ -742,6 +745,7 @@ contains ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then + allocate(name_(len_trim(name) + 1)) name_ = to_c_string(name) call read_double_c(obj_id, c_loc(name_), buffer, indep_) else @@ -764,6 +768,7 @@ contains ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then + allocate(name_(len_trim(name) + 1)) name_ = to_c_string(name) call read_double_c(obj_id, c_loc(name_), buffer, indep_) else @@ -888,6 +893,7 @@ contains ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then + allocate(name_(len_trim(name) + 1)) name_ = to_c_string(name) call read_int_c(obj_id, c_loc(name_), buffer_, indep_) else @@ -911,6 +917,7 @@ contains ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then + allocate(name_(len_trim(name) + 1)) name_ = to_c_string(name) call read_int_c(obj_id, c_loc(name_), buffer, indep_) else @@ -933,6 +940,7 @@ contains ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then + allocate(name_(len_trim(name) + 1)) name_ = to_c_string(name) call read_int_c(obj_id, c_loc(name_), buffer, indep_) else @@ -955,6 +963,7 @@ contains ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then + allocate(name_(len_trim(name) + 1)) name_ = to_c_string(name) call read_int_c(obj_id, c_loc(name_), buffer, indep_) else @@ -977,6 +986,7 @@ contains ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then + allocate(name_(len_trim(name) + 1)) name_ = to_c_string(name) call read_int_c(obj_id, c_loc(name_), buffer, indep_) else @@ -1025,6 +1035,7 @@ contains ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then + allocate(name_(len_trim(name) + 1)) name_ = to_c_string(name) call read_llong_c(obj_id, c_loc(name_), buffer_, indep_) else @@ -1120,6 +1131,7 @@ contains do i = 0, dims(1) - 1 n = len_trim(buffer(i+1)) + 1 buffer_(i*m+1 : i*m+n) = to_c_string(buffer(i+1)) + if (n < m) buffer_(i*m+n : i*m+m) = C_NULL_CHAR end do call write_string_c(group_id, 1, dims, m, to_c_string(name), & @@ -1350,7 +1362,7 @@ contains ! Allocate a C char array to get strings n = attribute_typesize(obj_id, to_c_string(name)) m = size(buffer) - allocate(buffer_((n+1)*m)) + allocate(buffer_(n*m)) ! Read attribute call read_attr_string_c(obj_id, to_c_string(name), n, buffer_) @@ -1359,7 +1371,7 @@ contains do i = 1, m buffer(i) = '' do j = 1, n - k = (i-1)*(n+1) + j + k = (i-1)*n + j if (buffer_(k) == C_NULL_CHAR) exit buffer(i)(j:j) = buffer_(k) end do @@ -1425,6 +1437,7 @@ contains ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then + allocate(name_(len_trim(name) + 1)) name_ = to_c_string(name) call read_complex_c(obj_id, c_loc(name_), buffer, indep_) else diff --git a/src/hdf5_interface.cpp b/src/hdf5_interface.cpp index cecf283ac1..d72f7d424d 100644 --- a/src/hdf5_interface.cpp +++ b/src/hdf5_interface.cpp @@ -392,7 +392,7 @@ object_name(hid_t obj_id) // Read and return name H5Iget_name(obj_id, buffer, size); - return {buffer, size}; + return buffer; } @@ -439,7 +439,9 @@ read_attr_string(hid_t obj_id, const char* name, size_t slen, char* buffer) { // Create datatype for a string hid_t datatype = H5Tcopy(H5T_C_S1); - H5Tset_size(datatype, slen + 1); + H5Tset_size(datatype, slen); + // numpy uses null-padding when writing fixed-length strings + H5Tset_strpad(datatype, H5T_STR_NULLPAD); // Read data into buffer read_attr(obj_id, name, datatype, buffer); @@ -503,7 +505,9 @@ read_string(hid_t obj_id, const char* name, size_t slen, char* buffer, bool inde { // Create datatype for a string hid_t datatype = H5Tcopy(H5T_C_S1); - H5Tset_size(datatype, slen + 1); + H5Tset_size(datatype, slen); + // numpy uses null-padding when writing fixed-length strings + H5Tset_strpad(datatype, H5T_STR_NULLPAD); // Read data into buffer read_dataset(obj_id, name, datatype, buffer, indep); diff --git a/src/initialize.cpp b/src/initialize.cpp index 19fc82c859..0c26135dc9 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -115,10 +115,9 @@ parse_command_line(int argc, char* argv[]) // Check what type of file this is hid_t file_id = file_open(argv[i], 'r', true); - size_t len = attribute_typesize(file_id, "filetype"); - read_attr_string(file_id, "filetype", len, buffer); + std::string filetype; + read_attribute(file_id, "filetype", filetype); file_close(file_id); - std::string filetype {buffer}; // Set path and flag for type of run if (filetype == "statepoint") { @@ -141,8 +140,7 @@ parse_command_line(int argc, char* argv[]) // Check file type is a source file file_id = file_open(argv[i+1], 'r', true); - len = attribute_typesize(file_id, "filetype"); - read_attr_string(file_id, "filetype", len, buffer); + read_attribute(file_id, "filetype", filetype); file_close(file_id); if (filetype != "source") { std::string msg {"Second file after restart flag must be a source file"}; diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp index 56f3392ff2..0bb464da07 100644 --- a/src/mgxs_interface.cpp +++ b/src/mgxs_interface.cpp @@ -208,7 +208,7 @@ void get_name_c(int index, int name_len, char* name) { // First blank out our input string - std::string str(name_len, ' '); + std::string str(name_len - 1, ' '); std::strcpy(name, str.c_str()); // Now get the data and copy to the C-string diff --git a/src/simulation.F90 b/src/simulation.F90 index 074d80e09f..356cca1a83 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -87,13 +87,6 @@ contains call initialize_batch() - ! Handle restart runs - if (restart_run .and. current_batch <= restart_batch) then - call replay_batch_history() - status = STATUS_EXIT_NORMAL - return - end if - ! ======================================================================= ! LOOP OVER GENERATIONS GENERATION_LOOP: do current_gen = 1, gen_per_batch @@ -204,10 +197,12 @@ contains ! Reset total starting particle weight used for normalizing tallies total_weight = ZERO - if (n_inactive > 0 .and. current_batch == 1) then + if ((n_inactive > 0 .and. current_batch == 1) .or. & + (restart_run .and. restart_batch < n_inactive .and. current_batch == restart_batch + 1)) then ! Turn on inactive timer call time_inactive % start() - elseif (current_batch == n_inactive + 1) then + elseif ((current_batch == n_inactive + 1) .or. & + (restart_run .and. restart_batch > n_inactive .and. current_batch == restart_batch + 1)) then ! Switch from inactive batch timer to active batch timer call time_inactive % stop() call time_active % start() @@ -386,46 +381,6 @@ contains end subroutine finalize_batch -!=============================================================================== -! REPLAY_BATCH_HISTORY displays keff and entropy for each generation within a -! batch using data read from a state point file -!=============================================================================== - - subroutine replay_batch_history - - ! Write message at beginning - if (current_batch == 1) then - call write_message("Replaying history from state point...", 6) - end if - - if (run_mode == MODE_EIGENVALUE) then - do current_gen = 1, gen_per_batch - call calculate_average_keff() - - ! print out batch keff - if (verbosity >= 7) then - if (current_gen < gen_per_batch) then - if (master) call print_generation() - else - if (master) call print_batch_keff() - end if - end if - end do - end if - - ! Increment n_realizations as would ordinarily be done in finalize_batch - if (reduce_tallies) then - n_realizations = n_realizations + 1 - else - n_realizations = n_realizations + n_procs - end if - - ! Write message at end - if (current_batch == restart_batch) then - call write_message("Resuming simulation...", 6) - end if - - end subroutine replay_batch_history !=============================================================================== ! INITIALIZE_SIMULATION @@ -488,23 +443,11 @@ contains ! file if (restart_run) then call load_state_point() + call write_message("Resuming simulation...", 6) else call initialize_source() end if - ! In restart, set the batch to begin from in order to reproduce the correct - ! average keff (used in sampling the fission bank). Use n_realizations from - ! the statepoint rather than n_inactive in case openmc_reset was called in - ! the previous run. - if (restart_run) then - if (reduce_tallies) then - current_batch = restart_batch - n_realizations - else - current_batch = restart_batch - n_realizations*n_procs - end if - n_realizations = 0 - end if - ! Display header if (master) then if (run_mode == MODE_FIXEDSOURCE) then diff --git a/src/state_point.F90 b/src/state_point.F90 index 526066c529..91d0f10908 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -16,7 +16,7 @@ module state_point use bank_header, only: Bank use cmfd_header use constants - use eigenvalue, only: openmc_get_keff + use eigenvalue, only: openmc_get_keff, k_sum use endf, only: reaction_name use error, only: fatal_error, warning, write_message use hdf5_interface @@ -749,6 +749,20 @@ contains ! Read number of realizations for global tallies call read_dataset(n_realizations, file_id, "n_realizations", indep=.true.) + ! Set k_sum, keff, and current_batch based on whether restart file is part + ! of active cycle or inactive cycle + if (restart_batch > n_inactive) then + do i = n_inactive + 1, restart_batch + k_sum(1) = k_sum(1) + k_generation % data(i) + k_sum(2) = k_sum(2) + k_generation % data(i)**2 + end do + n = gen_per_batch*n_realizations + keff = k_sum(1) / n + else + keff = k_generation % data(n) + end if + current_batch = restart_batch + ! Check to make sure source bank is present if (path_source_point == path_state_point .and. .not. source_present) then call fatal_error("Source bank must be contained in statepoint restart & diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index b0532a40c0..75779d1ad3 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -280,6 +280,7 @@ contains end subroutine tally_allocate_results + function tally_set_filters(this, filter_indices) result(err) class(TallyObject), intent(inout) :: this integer(C_INT32_T), intent(in) :: filter_indices(:) @@ -1032,6 +1033,7 @@ contains end if end function openmc_tally_set_scores + function openmc_tally_set_type(index, type) result(err) bind(C) ! Update the type of a tally that is already allocated integer(C_INT32_T), value, intent(in) :: index diff --git a/src/vector_header.F90 b/src/vector_header.F90 index 1ada0fba74..fd7726c79e 100644 --- a/src/vector_header.F90 +++ b/src/vector_header.F90 @@ -14,7 +14,7 @@ module vector_header procedure :: destroy => vector_destroy procedure :: add_value => vector_add_value procedure :: copy => vector_copy - ! TODO: procedure :: write => vector_write + procedure :: write => vector_write end type Vector contains @@ -88,4 +88,26 @@ contains end subroutine vector_copy +!=============================================================================== +! VECTOR_WRITE writes a vector to file +!=============================================================================== + + subroutine vector_write(self, filename) + + class(Vector), target, intent(inout) :: self ! vector instance + character(*), intent(in) :: filename ! filename to output to + + integer :: unit_ + integer :: i + + open(newunit=unit_, file=filename) + + do i = 1, self % n + write(unit_,*) i, self % data(i) + end do + + close(unit_) + + end subroutine vector_write + end module vector_header diff --git a/tests/regression_tests/plot_voxel/__init__.py b/tests/regression_tests/plot_voxel/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/regression_tests/plot_voxel/geometry.xml b/tests/regression_tests/plot_voxel/geometry.xml new file mode 100644 index 0000000000..83619d9f78 --- /dev/null +++ b/tests/regression_tests/plot_voxel/geometry.xml @@ -0,0 +1,13 @@ + + + + + + + + + + + + + diff --git a/tests/regression_tests/plot_voxel/materials.xml b/tests/regression_tests/plot_voxel/materials.xml new file mode 100644 index 0000000000..90b3542675 --- /dev/null +++ b/tests/regression_tests/plot_voxel/materials.xml @@ -0,0 +1,19 @@ + + + + + + + + + + + + + + + + + + + diff --git a/tests/regression_tests/plot_voxel/plots.xml b/tests/regression_tests/plot_voxel/plots.xml new file mode 100644 index 0000000000..833329b427 --- /dev/null +++ b/tests/regression_tests/plot_voxel/plots.xml @@ -0,0 +1,10 @@ + + + + + 50 50 10 + 0. 0. 0. + 20 20 10 + + + diff --git a/tests/regression_tests/plot_voxel/results_true.dat b/tests/regression_tests/plot_voxel/results_true.dat new file mode 100644 index 0000000000..41193c9bc3 --- /dev/null +++ b/tests/regression_tests/plot_voxel/results_true.dat @@ -0,0 +1 @@ +20a0c3598bb698efa4bba65c5e55d71f4385e5b1bcf0067964e45001c12f5c0a729935a96409c760dbd3ec439ae6c676fce77d6e578b41f8f3e0ac68c9277acb \ No newline at end of file diff --git a/tests/regression_tests/plot_voxel/settings.xml b/tests/regression_tests/plot_voxel/settings.xml new file mode 100644 index 0000000000..adf256d2d4 --- /dev/null +++ b/tests/regression_tests/plot_voxel/settings.xml @@ -0,0 +1,13 @@ + + + + plot + + + 5 4 3 + -10 -10 -10 + 10 10 10 + + 1 + + diff --git a/tests/regression_tests/plot_voxel/test.py b/tests/regression_tests/plot_voxel/test.py new file mode 100644 index 0000000000..d2767e2f18 --- /dev/null +++ b/tests/regression_tests/plot_voxel/test.py @@ -0,0 +1,62 @@ +import glob +import hashlib +import os +from subprocess import check_call + +import h5py +import openmc +import pytest + +from tests.testing_harness import TestHarness +from tests.regression_tests import config + + +vtk = pytest.importorskip('vtk') +class PlotVoxelTestHarness(TestHarness): + """Specialized TestHarness for running OpenMC voxel plot tests.""" + def __init__(self, plot_names): + super().__init__(None) + self._plot_names = plot_names + + def _run_openmc(self): + openmc.plot_geometry(openmc_exec=config['exe']) + + check_call(['../../../scripts/openmc-voxel-to-vtk'] + + glob.glob('plot_4.h5')) + + def _test_output_created(self): + """Make sure *.ppm has been created.""" + for fname in self._plot_names: + assert os.path.exists(fname), 'Plot output file does not exist.' + + def _cleanup(self): + super()._cleanup() + for fname in self._plot_names: + if os.path.exists(fname): + os.remove(fname) + + def _get_results(self): + """Return a string hash of the plot files.""" + outstr = bytes() + + for fname in self._plot_names: + if fname.endswith('.h5'): + # Add voxel data to results + with h5py.File(fname, 'r') as fh: + outstr += fh.attrs['filetype'] + outstr += fh.attrs['num_voxels'].tostring() + outstr += fh.attrs['lower_left'].tostring() + outstr += fh.attrs['voxel_width'].tostring() + outstr += fh['data'].value.tostring() + + # Hash the information and return. + sha512 = hashlib.sha512() + sha512.update(outstr) + outstr = sha512.hexdigest() + + return outstr + + +def test_plot_voxel(): + harness = PlotVoxelTestHarness(('plot_4.h5', 'plot.vti')) + harness.main() diff --git a/tests/unit_tests/test_capi.py b/tests/unit_tests/test_capi.py index 9d3bc106d2..ae2d927191 100644 --- a/tests/unit_tests/test_capi.py +++ b/tests/unit_tests/test_capi.py @@ -375,13 +375,13 @@ def test_restart(capi_init): openmc.capi.next_batch() keff0 = openmc.capi.keff() - # Restart the simulation from the statepoint and the 5 active batches. + # Restart the simulation from the statepoint and the 3 remaining active batches. openmc.capi.simulation_finalize() openmc.capi.hard_reset() openmc.capi.finalize() openmc.capi.init(args=('-r', 'restart_test.h5')) openmc.capi.simulation_init() - for i in range(5): + for i in range(3): openmc.capi.next_batch() keff1 = openmc.capi.keff() openmc.capi.simulation_finalize() diff --git a/tools/ci/travis-install.sh b/tools/ci/travis-install.sh index c10c222793..670f28f33b 100755 --- a/tools/ci/travis-install.sh +++ b/tools/ci/travis-install.sh @@ -23,8 +23,13 @@ fi # Build and install OpenMC executable python tools/ci/travis-install.py -# Install Python API in editable mode -pip install -e .[test] +if [[ $TRAVIS_PYTHON_VERSION == "3.7" ]]; then + # Install Python API in editable mode + pip install -e .[test] +else + # Conditionally install vtk + pip install -e .[test,vtk] +fi # For uploading to coveralls pip install coveralls