crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge remote-tracking branch 'upstream/develop' into cpp_tallies

Browse source
This commit is contained in:
Sterling Harper 2018-10-11 19:19:55 -04:00
commit 6a4cea41ed
29 changed files with 983 additions and 224 deletions
Download patch file Download diff file Expand all files Collapse all files

1
.gitignore vendored
View file

@ -69,6 +69,7 @@ scripts/G4EMLOW*/
# Images
*.ppm
*.voxel
*.vti
# PyCharm project configuration files
.idea

693
docs/source/methods/photon_physics.rst
View file

@ -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

24
include/openmc/hdf5_interface.h
View file

@ -1,9 +1,11 @@
#ifndef OPENMC_HDF5_INTERFACE_H
#define OPENMC_HDF5_INTERFACE_H
#include <algorithm> // for min
#include <array>
#include <complex>
#include <cstddef>
#include <cstring> // for strlen
#include <string>
#include <sstream>
#include <vector>
@ -162,13 +164,13 @@ void read_attribute(hid_t obj_id, const char* name, xt::xarray<T>& arr)
std::size_t size = 1;
for (const auto x : shape)
size *= x;
T* buffer = new T[size];
std::vector<T> buffer(size);
// Read data from attribute
read_attr(obj_id, name, H5TypeMap<T>::type_id, buffer);
read_attr(obj_id, name, H5TypeMap<T>::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<std::string>& 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<T>& arr, bool indep=false)
std::size_t size = 1;
for (const auto x : shape)
size *= x;
T* buffer = new T[size];
std::vector<T> buffer(size);
// Read data from attribute
read_dataset(dset, nullptr, H5TypeMap<T>::type_id, buffer, indep);
read_dataset(dset, nullptr, H5TypeMap<T>::type_id, buffer.data(), indep);
// Adapt into xarray
arr = xt::adapt(buffer, size, xt::acquire_ownership(), shape);
arr = xt::adapt(buffer, shape);
}
template <typename T>
@ -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<T> buffer(size);
// Read data from attribute
read_dataset(dset, nullptr, H5TypeMap<T>::type_id, buffer, indep);
read_dataset(dset, nullptr, H5TypeMap<T>::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);
}

1
include/openmc/search.h
View file

@ -14,6 +14,7 @@ template<class It, class T>
typename std::iterator_traits<It>::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;
}

24
openmc/capi/core.py
View file

@ -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.

5
openmc/capi/settings.py
View file

@ -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

85
scripts/openmc-voxel-to-silovtk
View file

@ -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()

58
scripts/openmc-voxel-to-vtk Executable file
View file

@ -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()

2
setup.py
View file

@ -64,7 +64,7 @@ kwargs = {
# Optional dependencies
'extras_require': {
'test': ['pytest', 'pytest-cov'],
'vtk': ['vtk', 'silomesh'],
'vtk': ['vtk'],
},
}

24
src/cmfd_data.F90
View file

@ -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

2
src/cmfd_execute.cpp
View file

@ -25,8 +25,6 @@ cmfd_populate_sourcecounts(int n_energy, const double* energies,
auto& m = meshes.at(index_cmfd_mesh);
xt::xarray<double> 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);
}

11
src/cmfd_solver.F90
View file

@ -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

17
src/hdf5_interface.F90
View file

@ -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

10
src/hdf5_interface.cpp
View file

@ -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);

8
src/initialize.cpp
View file

@ -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"};

2
src/mgxs_interface.cpp
View file

@ -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

67
src/simulation.F90
View file

@ -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

16
src/state_point.F90
View file

@ -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 &

2
src/tallies/tally_header.F90
View file

@ -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

24
src/vector_header.F90
View file

@ -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

0
tests/regression_tests/plot_voxel/__init__.py Normal file
View file

13
tests/regression_tests/plot_voxel/geometry.xml Normal file
View file

@ -0,0 +1,13 @@
<?xml version="1.0"?>
<geometry>
<surface id="1" type="z-cylinder" coeffs="0 0 2" />
<surface id="2" type="z-cylinder" coeffs="0 0 5" />
<surface id="3" type="z-cylinder" coeffs="0 0 10" boundary="vacuum" />
<surface id="4" type="y-plane" coeffs="5" boundary="vacuum" />
<cell id="1" material="1" region=" -1 -4" />
<cell id="2" material="3" region="1 -2 -4" />
<cell id="3" material="2" region="2 -3 -4" />
</geometry>

19
tests/regression_tests/plot_voxel/materials.xml Normal file
View file

@ -0,0 +1,19 @@
<?xml version="1.0"?>
<materials>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />
</material>
<material id="2">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />
</material>
<material id="3">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />
</material>
</materials>

10
tests/regression_tests/plot_voxel/plots.xml Normal file
View file

@ -0,0 +1,10 @@
<?xml version="1.0"?>
<plots>
<plot id="4" type="voxel">
<pixels>50 50 10</pixels>
<origin>0. 0. 0.</origin>
<width>20 20 10</width>
</plot>
</plots>

1
tests/regression_tests/plot_voxel/results_true.dat Normal file
View file

@ -0,0 +1 @@
20a0c3598bb698efa4bba65c5e55d71f4385e5b1bcf0067964e45001c12f5c0a729935a96409c760dbd3ec439ae6c676fce77d6e578b41f8f3e0ac68c9277acb

13
tests/regression_tests/plot_voxel/settings.xml Normal file
View file

@ -0,0 +1,13 @@
<?xml version="1.0"?>
<settings>
<run_mode>plot</run_mode>
<mesh id="1">
<dimension>5 4 3</dimension>
<lower_left>-10 -10 -10</lower_left>
<upper_right>10 10 10</upper_right>
</mesh>
<entropy_mesh>1</entropy_mesh>
</settings>

62
tests/regression_tests/plot_voxel/test.py Normal file
View file

@ -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()

4
tests/unit_tests/test_capi.py
View file

@ -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()

9
tools/ci/travis-install.sh
View file

@ -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