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

Merge remote-tracking branch 'gitee/develop' into virtual_lattice_0.15.2

Browse source
This commit is contained in:
skywalker_cn 2025-10-21 08:34:59 +00:00
commit ce605fa2d1
81 changed files with 1529 additions and 879 deletions
Download patch file Download diff file Expand all files Collapse all files

35
.github/workflows/ci.yml vendored
View file

@ -75,6 +75,7 @@ jobs:
LIBMESH: ${{ matrix.libmesh }}
NPY_DISABLE_CPU_FEATURES: "AVX512F AVX512_SKX"
OPENBLAS_NUM_THREADS: 1
PYTEST_ADDOPTS: --cov=openmc --cov-report=lcov:coverage-python.lcov
# libfabric complains about fork() as a result of using Python multiprocessing.
# We can work around it with RDMAV_FORK_SAFE=1 in libfabric < 1.13 and with
# FI_EFA_FORK_SAFE=1 in more recent versions.
@ -171,11 +172,37 @@ jobs:
uses: mxschmitt/action-tmate@v3
timeout-minutes: 10
- name: after_success
- name: Generate C++ coverage (gcovr)
shell: bash
run: |
cpp-coveralls -i src -i include -e src/external --exclude-pattern "/usr/*" --dump cpp_cov.json
coveralls --merge=cpp_cov.json --service=github
# Produce LCOV directly from gcov data in the build tree
gcovr \
--root "$GITHUB_WORKSPACE" \
--object-directory "$GITHUB_WORKSPACE/build" \
--filter "$GITHUB_WORKSPACE/src" \
--filter "$GITHUB_WORKSPACE/include" \
--exclude "$GITHUB_WORKSPACE/src/external/.*" \
--exclude "$GITHUB_WORKSPACE/src/include/openmc/external/.*" \
--gcov-ignore-errors source_not_found \
--gcov-ignore-errors output_error \
--gcov-ignore-parse-errors suspicious_hits.warn \
--print-summary \
--lcov -o coverage-cpp.lcov || true
- name: Merge C++ and Python coverage
shell: bash
run: |
# Merge C++ and Python LCOV into a single file for upload
cat coverage-cpp.lcov coverage-python.lcov > coverage.lcov
- name: Upload coverage to Coveralls
if: ${{ hashFiles('coverage.lcov') != '' }}
uses: coverallsapp/github-action@v2
with:
github-token: ${{ secrets.GITHUB_TOKEN }}
parallel: true
flag-name: C++ and Python
path-to-lcov: coverage.lcov
finish:
needs: main
@ -184,5 +211,5 @@ jobs:
- name: Coveralls Finished
uses: coverallsapp/github-action@v2
with:
github-token: ${{ secrets.github_token }}
github-token: ${{ secrets.GITHUB_TOKEN }}
parallel-finished: true

21
docs/source/io_formats/depletion_chain.rst
View file

@ -56,6 +56,27 @@ attributes:
.. _io_chain_reaction:
--------------------
``<source>`` Element
--------------------
The ``<source>`` element represents photon and electron sources associated with
the decay of a nuclide and contains information to construct an
:class:`openmc.stats.Univariate` object that represents this emission as an
energy distribution. This element has the following attributes:
:type:
The type of :class:`openmc.stats.Univariate` source term.
:particle:
The type of particle emitted, e.g., 'photon' or 'electron'
:parameters:
The parameters of the source term, e.g., for a
:class:`openmc.stats.Discrete` source, the energies (in [eV]) at which the
particles are emitted and their relative intensities in [Bq/atom] (in other
words, decay constants).
----------------------
``<reaction>`` Element
----------------------

10
docs/source/io_formats/settings.rst
View file

@ -178,6 +178,16 @@ history-based parallelism.
*Default*: false
--------------------------------
``<free_gas_threshold>`` Element
--------------------------------
The ``<free_gas_threshold>`` element specifies the energy multiplier, expressed
in units of :math:`kT`, that determines when the free gas scattering approach is
used for elastic scattering. Values must be positive.
*Default*: 400.0
-----------------------------------
``<generations_per_batch>`` Element
-----------------------------------

362
docs/source/methods/charged_particles_physics.rst Normal file
View file

@ -0,0 +1,362 @@
.. _methods_charged_particle_physics:
========================
Charged Particle Physics
========================
OpenMC neglects the spatial transport of charged particles (electrons and
positrons), assuming they deposit all their energy locally and produce
bremsstrahlung photons at their birth location. This approximation, called
thick-target bremsstrahlung (TTB) approximation is justified by the fact that
charged particles have much shorter stopping ranges compared to neutrons and
photons, especially in high-density materials.
-----------------------------
Charged Particle Interactions
-----------------------------
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
\begin{aligned}
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),
\end{aligned}
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. The
radiative stopping power for electrons is given by
.. math::
:label: radiative-stopping-power
S_{\text{rad}}(T) = n \frac{Z^2}{\beta^2} T \int_0^1 \chi(Z,T,\kappa)
d\kappa.
To obtain the radiative stopping power for positrons,
:eq:`radiative-stopping-power` is multiplied by :eq:`positron-factor`.
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 and
:math:`S_{\text{rad},i}(T)` is found for element :math:`i` using
:eq:`radiative-stopping-power`. 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. The Bethe
formula can be used to find the collision stopping power of the material:
.. math::
:label: material-collision-stopping-power
S_{\text{col}}(T) = \frac{2 \pi r_e^2 m_e c^2}{\beta^2} N_A \frac{Z}{A_M}
[\ln(T^2/I^2) + \ln(1 + \tau/2) + F(\tau) - \delta_F(T)],
where :math:`N_A` is Avogadro's number, :math:`A_M` is the molar mass,
:math:`\tau = T/m_e`, and :math:`F(\tau)` depends on the particle type. For
electrons,
.. math::
:label: F-electron
F_{-}(\tau) = (1 - \beta^2)[1 + \tau^2/8 - (2\tau + 1) \ln2],
while for positrons
.. math::
:label: F-positron
F_{+}(\tau) = 2\ln2 - (\beta^2/12)[23 + 14/(\tau + 2) + 10/(\tau + 2)^2 +
4/(\tau + 2)^3].
The density effect correction :math:`\delta_F` takes into account the reduction
of the collision stopping power due to the polarization of the material the
charged particle is passing through by the electric field of the particle.
It can be evaluated using the method described by Sternheimer_, where the
equation for :math:`\delta_F` is
.. math::
:label: density-effect-correction
\delta_F(\beta) = \sum_{i=1}^n f_i \ln[(l_i^2 + l^2)/l_i^2] -
l^2(1-\beta^2).
Here, :math:`f_i` is the oscillator strength of the :math:`i`-th transition,
given by :math:`f_i = n_i/Z`, where :math:`n_i` is the number of electrons in
the :math:`i`-th subshell. The frequency :math:`l` is the solution of the
equation
.. math::
:label: density-effect-l
\frac{1}{\beta^2} - 1 = \sum_{i=1}^{n} \frac{f_i}{\bar{\nu}_i^2 + l^2},
where :math:`\bar{v}_i` is defined as
.. math::
:label: density-effect-nubar
\bar{\nu}_i = h\nu_i \rho / h\nu_p.
The plasma energy :math:`h\nu_p` of the medium is given by
.. math::
:label: plasma-frequency
h\nu_p = \sqrt{\frac{(hc)^2 r_e \rho_m N_A Z}{\pi A}},
where :math:`A` is the atomic weight and :math:`\rho_m` is the density of the
material. In :eq:`density-effect-nubar`, :math:`h\nu_i` is the oscillator
energy, and :math:`\rho` is an adjustment factor introduced to give agreement
between the experimental values of the oscillator energies and the mean
excitation energy. The :math:`l_i` in :eq:`density-effect-correction` are
defined as
.. math::
:label: density-effect-li
\begin{aligned}
l_i &= (\bar{\nu}_i^2 + 2/3f_i)^{1/2} ~~~~&\text{for}~~ \bar{\nu}_i > 0 \\
l_n &= f_n^{1/2} ~~~~&\text{for}~~ \bar{\nu}_n = 0,
\end{aligned}
where the second case applies to conduction electrons. For a conductor,
:math:`f_n` is given by :math:`n_c/Z`, where :math:`n_c` is the effective
number of conduction electrons, and :math:`v_n = 0`. The adjustment factor
:math:`\rho` is determined using the equation for the mean excitation energy:
.. math::
:label: mean-excitation-energy
\ln I = \sum_{i=1}^{n-1} f_i \ln[(h\nu_i\rho)^2 + 2/3f_i(h\nu_p)^2]^{1/2} +
f_n \ln (h\nu_pf_n^{1/2}).
.. _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.
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.
.. _Kaltiaisenaho: https://aaltodoc.aalto.fi/bitstream/handle/123456789/21004/master_Kaltiaisenaho_Toni_2016.pdf
.. _Salvat: https://doi.org/10.1787/32da5043-en
.. _Sternheimer: https://doi.org/10.1103/PhysRevB.26.6067

3
docs/source/methods/index.rst
View file

@ -14,6 +14,7 @@ Theory and Methodology
random_numbers
neutron_physics
photon_physics
charged_particles_physics
tallies
eigenvalue
depletion
@ -21,4 +22,4 @@ Theory and Methodology
parallelization
cmfd
variance_reduction
random_ray
random_ray

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

@ -667,342 +667,6 @@ and Auger electrons:
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
\begin{aligned}
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),
\end{aligned}
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. The
radiative stopping power for electrons is given by
.. math::
:label: radiative-stopping-power
S_{\text{rad}}(T) = n \frac{Z^2}{\beta^2} T \int_0^1 \chi(Z,T,\kappa)
d\kappa.
To obtain the radiative stopping power for positrons,
:eq:`radiative-stopping-power` is multiplied by :eq:`positron-factor`.
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 and
:math:`S_{\text{rad},i}(T)` is found for element :math:`i` using
:eq:`radiative-stopping-power`. 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. The Bethe
formula can be used to find the collision stopping power of the material:
.. math::
:label: material-collision-stopping-power
S_{\text{col}}(T) = \frac{2 \pi r_e^2 m_e c^2}{\beta^2} N_A \frac{Z}{A_M}
[\ln(T^2/I^2) + \ln(1 + \tau/2) + F(\tau) - \delta_F(T)],
where :math:`N_A` is Avogadro's number, :math:`A_M` is the molar mass,
:math:`\tau = T/m_e`, and :math:`F(\tau)` depends on the particle type. For
electrons,
.. math::
:label: F-electron
F_{-}(\tau) = (1 - \beta^2)[1 + \tau^2/8 - (2\tau + 1) \ln2],
while for positrons
.. math::
:label: F-positron
F_{+}(\tau) = 2\ln2 - (\beta^2/12)[23 + 14/(\tau + 2) + 10/(\tau + 2)^2 +
4/(\tau + 2)^3].
The density effect correction :math:`\delta_F` takes into account the reduction
of the collision stopping power due to the polarization of the material the
charged particle is passing through by the electric field of the particle.
It can be evaluated using the method described by Sternheimer_, where the
equation for :math:`\delta_F` is
.. math::
:label: density-effect-correction
\delta_F(\beta) = \sum_{i=1}^n f_i \ln[(l_i^2 + l^2)/l_i^2] -
l^2(1-\beta^2).
Here, :math:`f_i` is the oscillator strength of the :math:`i`-th transition,
given by :math:`f_i = n_i/Z`, where :math:`n_i` is the number of electrons in
the :math:`i`-th subshell. The frequency :math:`l` is the solution of the
equation
.. math::
:label: density-effect-l
\frac{1}{\beta^2} - 1 = \sum_{i=1}^{n} \frac{f_i}{\bar{\nu}_i^2 + l^2},
where :math:`\bar{v}_i` is defined as
.. math::
:label: density-effect-nubar
\bar{\nu}_i = h\nu_i \rho / h\nu_p.
The plasma energy :math:`h\nu_p` of the medium is given by
.. math::
:label: plasma-frequency
h\nu_p = \sqrt{\frac{(hc)^2 r_e \rho_m N_A Z}{\pi A}},
where :math:`A` is the atomic weight and :math:`\rho_m` is the density of the
material. In :eq:`density-effect-nubar`, :math:`h\nu_i` is the oscillator
energy, and :math:`\rho` is an adjustment factor introduced to give agreement
between the experimental values of the oscillator energies and the mean
excitation energy. The :math:`l_i` in :eq:`density-effect-correction` are
defined as
.. math::
:label: density-effect-li
\begin{aligned}
l_i &= (\bar{\nu}_i^2 + 2/3f_i)^{1/2} ~~~~&\text{for}~~ \bar{\nu}_i > 0 \\
l_n &= f_n^{1/2} ~~~~&\text{for}~~ \bar{\nu}_n = 0,
\end{aligned}
where the second case applies to conduction electrons. For a conductor,
:math:`f_n` is given by :math:`n_c/Z`, where :math:`n_c` is the effective
number of conduction electrons, and :math:`v_n = 0`. The adjustment factor
:math:`\rho` is determined using the equation for the mean excitation energy:
.. math::
:label: mean-excitation-energy
\ln I = \sum_{i=1}^{n-1} f_i \ln[(h\nu_i\rho)^2 + 2/3f_i(h\nu_p)^2]^{1/2} +
f_n \ln (h\nu_pf_n^{1/2}).
.. _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.
.. _photon_production:
@ -1070,5 +734,3 @@ emitted photon.
.. _Kaltiaisenaho: https://aaltodoc.aalto.fi/bitstream/handle/123456789/21004/master_Kaltiaisenaho_Toni_2016.pdf
.. _Salvat: https://doi.org/10.1787/32da5043-en
.. _Sternheimer: https://doi.org/10.1103/PhysRevB.26.6067

5
include/openmc/constants.h
View file

@ -68,6 +68,11 @@ constexpr double MIN_HITS_PER_BATCH {1.5};
// prevent extremely large adjoint source terms from being generated.
constexpr double ZERO_FLUX_CUTOFF {1e-22};
// The minimum macroscopic cross section value considered non-void for the
// random ray solver. Materials with any group with a cross section below this
// value will be converted to pure void.
constexpr double MINIMUM_MACRO_XS {1e-6};
// ============================================================================
// MATH AND PHYSICAL CONSTANTS

4
include/openmc/nuclide.h
View file

@ -164,8 +164,8 @@ namespace data {
// Minimum/maximum transport energy for each particle type. Order corresponds to
// that of the ParticleType enum
extern array<double, 2> energy_min;
extern array<double, 2> energy_max;
extern array<double, 4> energy_min;
extern array<double, 4> energy_max;
//! Minimum temperature in [K] that nuclide data is available at
extern double temperature_min;

7
include/openmc/particle_data.h
View file

@ -154,9 +154,10 @@ struct NuclideMicroXS {
// Energy and temperature last used to evaluate these cross sections. If
// these values have changed, then the cross sections must be re-evaluated.
double last_E {0.0}; //!< Last evaluated energy
double last_sqrtkT {0.0}; //!< Last temperature in sqrt(Boltzmann constant
//!< * temperature (eV))
double last_E {0.0}; //!< Last evaluated energy
double last_sqrtkT {0.0}; //!< Last temperature in sqrt(Boltzmann constant
//!< * temperature (eV))
double ncrystal_xs {-1.0}; //!< NCrystal cross section
};
//==============================================================================

7
include/openmc/physics.h
View file

@ -10,13 +10,6 @@
namespace openmc {
//==============================================================================
// Constants
//==============================================================================
// Monoatomic ideal-gas scattering treatment threshold
constexpr double FREE_GAS_THRESHOLD {400.0};
//==============================================================================
// Non-member functions
//==============================================================================

47
include/openmc/random_ray/flat_source_domain.h
View file

@ -27,8 +27,9 @@ public:
//----------------------------------------------------------------------------
// Methods
virtual void update_neutron_source(double k_eff);
double compute_k_eff(double k_eff_old) const;
virtual void update_single_neutron_source(SourceRegionHandle& srh);
virtual void update_all_neutron_sources();
void compute_k_eff();
virtual void normalize_scalar_flux_and_volumes(
double total_active_distance_per_iteration);
@ -41,7 +42,7 @@ public:
void output_to_vtk() const;
void convert_external_sources();
void count_external_source_regions();
void set_adjoint_sources(const vector<double>& forward_flux);
void set_adjoint_sources();
void flux_swap();
virtual double evaluate_flux_at_point(Position r, int64_t sr, int g) const;
double compute_fixed_source_normalization_factor() const;
@ -54,9 +55,8 @@ public:
bool is_target_void);
void apply_mesh_to_cell_and_children(int32_t i_cell, int32_t mesh_idx,
int32_t target_material_id, bool is_target_void);
void prepare_base_source_regions();
SourceRegionHandle get_subdivided_source_region_handle(
int64_t sr, int mesh_bin, Position r, double dist, Direction u);
SourceRegionKey sr_key, Position r, Direction u);
void finalize_discovered_source_regions();
void apply_transport_stabilization();
int64_t n_source_regions() const
@ -67,6 +67,10 @@ public:
{
return source_regions_.n_source_regions() * negroups_;
}
int64_t lookup_base_source_region_idx(const GeometryState& p) const;
SourceRegionKey lookup_source_region_key(const GeometryState& p) const;
int64_t lookup_mesh_bin(int64_t sr, Position r) const;
int lookup_mesh_idx(int64_t sr) const;
//----------------------------------------------------------------------------
// Static Data members
@ -86,6 +90,7 @@ public:
//----------------------------------------------------------------------------
// Public Data members
double k_eff_ {1.0}; // Eigenvalue
bool mapped_all_tallies_ {false}; // If all source regions have been visited
int64_t n_external_source_regions_ {0}; // Total number of source regions with
@ -110,14 +115,6 @@ public:
// The abstract container holding all source region-specific data
SourceRegionContainer source_regions_;
// Base source region container. When source region subdivision via mesh
// is in use, this container holds the original (non-subdivided) material
// filled cell instance source regions. These are useful as they can be
// initialized with external source and mesh domain information ahead of time.
// Then, dynamically discovered source regions can be initialized by cloning
// their base region.
SourceRegionContainer base_source_regions_;
// Parallel hash map holding all source regions discovered during
// a single iteration. This is a threadsafe data structure that is cleaned
// out after each iteration and stored in the "source_regions_" container.
@ -134,8 +131,17 @@ public:
// Map that relates a SourceRegionKey to the external source index. This map
// is used to check if there are any point sources within a subdivided source
// region at the time it is discovered.
std::unordered_map<SourceRegionKey, int64_t, SourceRegionKey::HashFunctor>
point_source_map_;
std::unordered_map<SourceRegionKey, vector<int>, SourceRegionKey::HashFunctor>
external_point_source_map_;
// Map that relates a base source region index to the external source index.
// This map is used to check if there are any volumetric sources within a
// subdivided source region at the time it is discovered.
std::unordered_map<int64_t, vector<int>> external_volumetric_source_map_;
// Map that relates a base source region index to a mesh index. This map
// is used to check which subdivision mesh is present in a source region.
std::unordered_map<int64_t, int> mesh_map_;
// If transport corrected MGXS data is being used, there may be negative
// in-group scattering cross sections that can result in instability in MOC
@ -147,12 +153,11 @@ protected:
//----------------------------------------------------------------------------
// Methods
void apply_external_source_to_source_region(
Discrete* discrete, double strength_factor, SourceRegionHandle& srh);
void apply_external_source_to_cell_instances(int32_t i_cell,
Discrete* discrete, double strength_factor, int target_material_id,
const vector<int32_t>& instances);
void apply_external_source_to_cell_and_children(int32_t i_cell,
Discrete* discrete, double strength_factor, int32_t target_material_id);
int src_idx, SourceRegionHandle& srh);
void apply_external_source_to_cell_instances(int32_t i_cell, int src_idx,
int target_material_id, const vector<int32_t>& instances);
void apply_external_source_to_cell_and_children(
int32_t i_cell, int src_idx, int32_t target_material_id);
virtual void set_flux_to_flux_plus_source(int64_t sr, double volume, int g);
void set_flux_to_source(int64_t sr, int g);
virtual void set_flux_to_old_flux(int64_t sr, int g);

2
include/openmc/random_ray/linear_source_domain.h
View file

@ -20,7 +20,7 @@ class LinearSourceDomain : public FlatSourceDomain {
public:
//----------------------------------------------------------------------------
// Methods
void update_neutron_source(double k_eff) override;
void update_single_neutron_source(SourceRegionHandle& srh) override;
void normalize_scalar_flux_and_volumes(
double total_active_distance_per_iteration) override;

1
include/openmc/random_ray/random_ray.h
View file

@ -48,7 +48,6 @@ public:
static double distance_active_; // Active ray length
static unique_ptr<Source> ray_source_; // Starting source for ray sampling
static RandomRaySourceShape source_shape_; // Flag for linear source
static bool mesh_subdivision_enabled_; // Flag for mesh subdivision
static RandomRaySampleMethod sample_method_; // Flag for sampling method
//----------------------------------------------------------------------------

9
include/openmc/random_ray/random_ray_simulation.h
View file

@ -21,11 +21,7 @@ public:
// Methods
void compute_segment_correction_factors();
void apply_fixed_sources_and_mesh_domains();
void prepare_fixed_sources_adjoint(vector<double>& forward_flux,
SourceRegionContainer& forward_source_regions,
SourceRegionContainer& forward_base_source_regions,
std::unordered_map<SourceRegionKey, int64_t, SourceRegionKey::HashFunctor>&
forward_source_region_map);
void prepare_fixed_sources_adjoint();
void simulate();
void output_simulation_results() const;
void instability_check(
@ -45,9 +41,6 @@ private:
// Contains all flat source region data
unique_ptr<FlatSourceDomain> domain_;
// Random ray eigenvalue
double k_eff_ {1.0};
// Tracks the average FSR miss rate for analysis and reporting
double avg_miss_rate_ {0.0};

1
include/openmc/random_ray/source_region.h
View file

@ -308,7 +308,6 @@ public:
//----------------------------------------------------------------------------
// Constructors
SourceRegion(int negroups, bool is_linear);
SourceRegion(const SourceRegionHandle& handle, int64_t parent_sr);
SourceRegion() = default;
//----------------------------------------------------------------------------

2
include/openmc/settings.h
View file

@ -147,6 +147,8 @@ extern std::unordered_set<int>
source_write_surf_id; //!< Surface ids where sources will be written
extern double source_rejection_fraction; //!< Minimum fraction of source sites
//!< that must be accepted
extern double free_gas_threshold; //!< Threshold multiplier for free gas
//!< scattering treatment
extern int
max_history_splits; //!< maximum number of particle splits for weight windows

2
openmc/data/data.py
View file

@ -324,7 +324,7 @@ def atomic_mass(isotope):
# isotopes of their element (e.g. C0), calculate the atomic mass as
# the sum of the atomic mass times the natural abundance of the isotopes
# that make up the element.
for element in ['C', 'Zn', 'Pt', 'Os', 'Tl']:
for element in ['C', 'Zn', 'Pt', 'Os', 'Tl', 'V']:
isotope_zero = element.lower() + '0'
_ATOMIC_MASS[isotope_zero] = 0.
for iso, abundance in isotopes(element):

2
openmc/data/decay.py
View file

@ -591,7 +591,7 @@ def decay_photon_energy(nuclide: str) -> Univariate | None:
openmc.stats.Univariate or None
Distribution of energies in [eV] of photons emitted from decay, or None
if no photon source exists. Note that the probabilities represent
intensities, given as [Bq].
intensities, given as [Bq/atom] (in other words, decay constants).
"""
if not _DECAY_PHOTON_ENERGY:
chain_file = openmc.config.get('chain_file')

13
openmc/deplete/d1s.py
View file

@ -108,14 +108,15 @@ def time_correction_factors(
# Create a 2D array for the time correction factors
h = np.zeros((n_timesteps, n_nuclides))
for i, (dt, rate) in enumerate(zip(timesteps, source_rates)):
# Precompute the exponential terms. Since (1 - exp(-x)) is susceptible to
# roundoff error, use expm1 instead (which computes exp(x) - 1)
g = np.exp(-decay_rate*dt)
one_minus_g = -np.expm1(-decay_rate*dt)
# Precompute all exponential terms with same shape as h
decay_dt = decay_rate[np.newaxis, :] * timesteps[:, np.newaxis]
g = np.exp(-decay_dt)
one_minus_g = -np.expm1(-decay_dt)
# Apply recurrence relation step by step
for i in range(len(timesteps)):
# Eq. (4) in doi:10.1016/j.fusengdes.2019.111399
h[i + 1] = rate*one_minus_g + h[i]*g
h[i + 1] = source_rates[i] * one_minus_g[i] + h[i] * g[i]
return {nuclides[i]: h[:, i] for i in range(n_nuclides)}

2
openmc/plotter.py
View file

@ -501,7 +501,7 @@ def _calculate_cexs_nuclide(this, types, temperature=294., sab_name=None,
elif ncrystal_cfg:
import NCrystal
nc_scatter = NCrystal.createScatter(ncrystal_cfg)
nc_func = nc_scatter.crossSectionNonOriented
nc_func = nc_scatter.xsect
nc_emax = 5 # eV # this should be obtained from NCRYSTAL_MAX_ENERGY
energy_grid = np.union1d(np.geomspace(min(energy_grid),
1.1*nc_emax,

58
openmc/settings.py
View file

@ -84,6 +84,10 @@ class Settings:
history-based parallelism.
.. versionadded:: 0.12
free_gas_threshold : float
Energy multiplier (in units of :math:`kT`) below which the free gas
scattering treatment is applied for elastic scattering. If not
specified, a value of 400.0 is used.
generations_per_batch : int
Number of generations per batch
ifp_n_generation : int
@ -376,6 +380,7 @@ class Settings:
self._seed = None
self._stride = None
self._survival_biasing = None
self._free_gas_threshold = None
# Shannon entropy mesh
self._entropy_mesh = None
@ -1255,6 +1260,17 @@ class Settings:
cv.check_less_than('source_rejection_fraction', source_rejection_fraction, 1)
self._source_rejection_fraction = source_rejection_fraction
@property
def free_gas_threshold(self) -> float | None:
return self._free_gas_threshold
@free_gas_threshold.setter
def free_gas_threshold(self, free_gas_threshold: float | None):
if free_gas_threshold is not None:
cv.check_type('free gas threshold', free_gas_threshold, Real)
cv.check_greater_than('free gas threshold', free_gas_threshold, 0.0)
self._free_gas_threshold = free_gas_threshold
def _create_run_mode_subelement(self, root):
elem = ET.SubElement(root, "run_mode")
elem.text = self._run_mode.value
@ -1641,6 +1657,7 @@ class Settings:
if mesh_memo is not None:
mesh_memo.add(ww.mesh.id)
def _create_weight_windows_on_subelement(self, root):
if self._weight_windows_on is not None:
elem = ET.SubElement(root, "weight_windows_on")
elem.text = str(self._weight_windows_on).lower()
@ -1714,9 +1731,15 @@ class Settings:
domain_elem = ET.SubElement(mesh_elem, 'domain')
domain_elem.set('id', str(domain.id))
domain_elem.set('type', domain.__class__.__name__.lower())
if mesh_memo is not None and mesh.id not in mesh_memo:
# See if a <mesh> element already exists -- if not, add it
path = f"./mesh[@id='{mesh.id}']"
if root.find(path) is None:
root.append(mesh.to_xml_element())
mesh_memo.add(mesh.id)
if mesh_memo is not None:
mesh_memo.add(mesh.id)
elif isinstance(value, bool):
subelement = ET.SubElement(element, key)
subelement.text = str(value).lower()
else:
subelement = ET.SubElement(element, key)
subelement.text = str(value)
@ -1726,6 +1749,11 @@ class Settings:
element = ET.SubElement(root, "source_rejection_fraction")
element.text = str(self._source_rejection_fraction)
def _create_free_gas_threshold_subelement(self, root):
if self._free_gas_threshold is not None:
element = ET.SubElement(root, "free_gas_threshold")
element.text = str(self._free_gas_threshold)
def _eigenvalue_from_xml_element(self, root):
elem = root.find('eigenvalue')
if elem is not None:
@ -2074,10 +2102,16 @@ class Settings:
ww = WeightWindows.from_xml_element(elem, meshes)
self.weight_windows.append(ww)
def _weight_windows_on_from_xml_element(self, root):
text = get_text(root, 'weight_windows_on')
if text is not None:
self.weight_windows_on = text in ('true', '1')
def _weight_windows_file_from_xml_element(self, root):
text = get_text(root, 'weight_windows_file')
if text is not None:
self.weight_windows_file = text
def _weight_window_checkpoints_from_xml_element(self, root):
elem = root.find('weight_window_checkpoints')
if elem is None:
@ -2103,7 +2137,7 @@ class Settings:
if text is not None:
self.max_tracks = int(text)
def _random_ray_from_xml_element(self, root):
def _random_ray_from_xml_element(self, root, meshes=None):
elem = root.find('random_ray')
if elem is not None:
self.random_ray = {}
@ -2130,7 +2164,11 @@ class Settings:
elif child.tag == 'source_region_meshes':
self.random_ray['source_region_meshes'] = []
for mesh_elem in child.findall('mesh'):
mesh = MeshBase.from_xml_element(mesh_elem)
mesh_id = int(get_text(mesh_elem, 'id'))
if meshes and mesh_id in meshes:
mesh = meshes[mesh_id]
else:
mesh = MeshBase.from_xml_element(mesh_elem)
domains = []
for domain_elem in mesh_elem.findall('domain'):
domain_id = int(get_text(domain_elem, "id"))
@ -2154,6 +2192,11 @@ class Settings:
if text is not None:
self.source_rejection_fraction = float(text)
def _free_gas_threshold_from_xml_element(self, root):
text = get_text(root, 'free_gas_threshold')
if text is not None:
self.free_gas_threshold = float(text)
def to_xml_element(self, mesh_memo=None):
"""Create a 'settings' element to be written to an XML file.
@ -2214,6 +2257,7 @@ class Settings:
self._create_log_grid_bins_subelement(element)
self._create_write_initial_source_subelement(element)
self._create_weight_windows_subelement(element, mesh_memo)
self._create_weight_windows_on_subelement(element)
self._create_weight_window_generators_subelement(element, mesh_memo)
self._create_weight_windows_file_element(element)
self._create_weight_window_checkpoints_subelement(element)
@ -2223,6 +2267,7 @@ class Settings:
self._create_random_ray_subelement(element, mesh_memo)
self._create_use_decay_photons_subelement(element)
self._create_source_rejection_fraction_subelement(element)
self._create_free_gas_threshold_subelement(element)
# Clean the indentation in the file to be user-readable
clean_indentation(element)
@ -2324,14 +2369,17 @@ class Settings:
settings._log_grid_bins_from_xml_element(elem)
settings._write_initial_source_from_xml_element(elem)
settings._weight_windows_from_xml_element(elem, meshes)
settings._weight_windows_on_from_xml_element(elem)
settings._weight_windows_file_from_xml_element(elem)
settings._weight_window_generators_from_xml_element(elem, meshes)
settings._weight_window_checkpoints_from_xml_element(elem)
settings._max_history_splits_from_xml_element(elem)
settings._max_tracks_from_xml_element(elem)
settings._max_secondaries_from_xml_element(elem)
settings._random_ray_from_xml_element(elem)
settings._random_ray_from_xml_element(elem, meshes)
settings._use_decay_photons_from_xml_element(elem)
settings._source_rejection_fraction_from_xml_element(elem)
settings._free_gas_threshold_from_xml_element(elem)
return settings

7
openmc/source.py
View file

@ -260,7 +260,7 @@ class IndependentSource(SourceBase):
time distribution of source sites
strength : float
Strength of the source
particle : {'neutron', 'photon'}
particle : {'neutron', 'photon', 'electron', 'positron'}
Source particle type
domains : iterable of openmc.Cell, openmc.Material, or openmc.Universe
Domains to reject based on, i.e., if a sampled spatial location is not
@ -299,7 +299,7 @@ class IndependentSource(SourceBase):
.. versionadded:: 0.14.0
particle : {'neutron', 'photon'}
particle : {'neutron', 'photon', 'electron', 'positron'}
Source particle type
constraints : dict
Constraints on sampled source particles. Valid keys include
@ -404,7 +404,8 @@ class IndependentSource(SourceBase):
@particle.setter
def particle(self, particle):
cv.check_value('source particle', particle, ['neutron', 'photon'])
cv.check_value('source particle', particle,
['neutron', 'photon', 'electron', 'positron'])
self._particle = particle
def populate_xml_element(self, element):

9
openmc/statepoint.py
View file

@ -536,7 +536,7 @@ class StatePoint:
def get_tally(self, scores=[], filters=[], nuclides=[],
name=None, id=None, estimator=None, exact_filters=False,
exact_nuclides=False, exact_scores=False,
multiply_density=None, derivative=None):
multiply_density=None, derivative=None, filter_type=None):
"""Finds and returns a Tally object with certain properties.
This routine searches the list of Tallies and returns the first Tally
@ -580,6 +580,9 @@ class StatePoint:
to the same value as this parameter.
derivative : openmc.TallyDerivative, optional
TallyDerivative object to match.
filter_type : type, optional
If not None, the Tally must have at least one Filter that is an
instance of this type. For example `openmc.MeshFilter`.
Returns
-------
@ -653,6 +656,10 @@ class StatePoint:
if not contains_filters:
continue
if filter_type is not None:
if not any(isinstance(f, filter_type) for f in test_tally.filters):
continue
# Determine if Tally has the queried Nuclide(s)
if nuclides:
if not all(nuclide in test_tally.nuclides for nuclide in nuclides):

76
openmc/stats/univariate.py
View file

@ -295,6 +295,20 @@ class Discrete(Univariate):
"""
return np.sum(self.p)
def mean(self) -> float:
"""Return mean of the discrete distribution
The mean is the weighted average of the discrete values.
.. versionadded:: 0.15.3
Returns
-------
float
Mean of discrete distribution
"""
return np.sum(self.x * self.p) / np.sum(self.p)
def clip(self, tolerance: float = 1e-6, inplace: bool = False) -> Discrete:
r"""Remove low-importance points from discrete distribution.
@ -413,6 +427,18 @@ class Uniform(Univariate):
rng = np.random.RandomState(seed)
return rng.uniform(self.a, self.b, n_samples)
def mean(self) -> float:
"""Return mean of the uniform distribution
.. versionadded:: 0.15.3
Returns
-------
float
Mean of uniform distribution
"""
return 0.5 * (self.a + self.b)
def to_xml_element(self, element_name: str):
"""Return XML representation of the uniform distribution
@ -1123,7 +1149,7 @@ class Tabular(Univariate):
"""
interpolation = get_text(elem, 'interpolation')
params = get_elem_list(elem, "parameters", float)
params = get_elem_list(elem, "parameters", float)
m = (len(params) + 1)//2 # +1 for when len(params) is odd
x = params[:m]
p = params[m:]
@ -1347,6 +1373,30 @@ class Mixture(Univariate):
for p, dist in zip(self.probability, self.distribution)
])
def mean(self) -> float:
"""Return mean of the mixture distribution
The mean is the weighted average of the means of the component
distributions, weighted by probability * integral.
.. versionadded:: 0.15.3
Returns
-------
float
Mean of the mixture distribution
"""
# Weight each component by its probability and integral
weights = [p*dist.integral() for p, dist in
zip(self.probability, self.distribution)]
total_weight = sum(weights)
if total_weight == 0:
return 0.0
return sum([w*dist.mean() for w, dist in
zip(weights, self.distribution)]) / total_weight
def clip(self, tolerance: float = 1e-6, inplace: bool = False) -> Mixture:
r"""Remove low-importance points / distributions
@ -1369,14 +1419,14 @@ class Mixture(Univariate):
Distribution with low-importance points / distributions removed
"""
# Determine integral of original distribution to compare later
original_integral = self.integral()
# Calculate mean * integral for original distribution to compare later.
original_mean_integral = self.mean() * self.integral()
# Determine indices for any distributions that contribute non-negligibly
# to overall intensity
intensities = [prob*dist.integral() for prob, dist in
zip(self.probability, self.distribution)]
indices = _intensity_clip(intensities, tolerance=tolerance)
# to overall mean * integral
mean_integrals = [prob*dist.mean()*dist.integral() for prob, dist in
zip(self.probability, self.distribution)]
indices = _intensity_clip(mean_integrals, tolerance=tolerance)
# Clip mixture of distributions
probability = self.probability[indices]
@ -1397,12 +1447,14 @@ class Mixture(Univariate):
# Create new distribution
new_dist = type(self)(probability, distribution)
# Show warning if integral of new distribution is not within
# tolerance of original
diff = (original_integral - new_dist.integral())/original_integral
# Show warning if mean * integral of new distribution is not within
# tolerance of original. For energy distributions, mean * integral
# represents total energy.
new_mean_integral = new_dist.mean() * new_dist.integral()
diff = (original_mean_integral - new_mean_integral)/original_mean_integral
if diff > tolerance:
warn("Clipping mixture distribution resulted in an integral that is "
f"lower by a fraction of {diff} when tolerance={tolerance}.")
warn("Clipping mixture distribution resulted in a mean*integral "
f"that is lower by a fraction of {diff} when tolerance={tolerance}.")
return new_dist

11
pyproject.toml
View file

@ -48,8 +48,15 @@ docs = [
"sphinxcontrib-svg2pdfconverter",
"sphinx-rtd-theme"
]
test = ["packaging", "pytest", "pytest-cov", "colorama", "openpyxl"]
ci = ["cpp-coveralls", "coveralls"]
test = [
"packaging",
"pytest",
"pytest-cov>=4.0",
"pytest-rerunfailures",
"colorama",
"openpyxl",
]
ci = ["coverage>=7.4", "gcovr>=7.2"]
vtk = ["vtk"]
[project.urls]

5
src/finalize.cpp
View file

@ -85,6 +85,7 @@ int openmc_finalize()
settings::time_cutoff = {INFTY, INFTY, INFTY, INFTY};
settings::entropy_on = false;
settings::event_based = false;
settings::free_gas_threshold = 400.0;
settings::gen_per_batch = 1;
settings::legendre_to_tabular = true;
settings::legendre_to_tabular_points = -1;
@ -155,8 +156,8 @@ int openmc_finalize()
simulation::entropy_mesh = nullptr;
simulation::ufs_mesh = nullptr;
data::energy_max = {INFTY, INFTY};
data::energy_min = {0.0, 0.0};
data::energy_max = {INFTY, INFTY, INFTY, INFTY};
data::energy_min = {0.0, 0.0, 0.0, 0.0};
data::temperature_min = 0.0;
data::temperature_max = INFTY;
model::root_universe = -1;

4
src/nuclide.cpp
View file

@ -31,8 +31,8 @@ namespace openmc {
//==============================================================================
namespace data {
array<double, 2> energy_min {0.0, 0.0};
array<double, 2> energy_max {INFTY, INFTY};
array<double, 4> energy_min {0.0, 0.0, 0.0, 0.0};
array<double, 4> energy_max {INFTY, INFTY, INFTY, INFTY};
double temperature_min {INFTY};
double temperature_max {0.0};
std::unordered_map<std::string, int> nuclide_map;

6
src/particle.cpp
View file

@ -229,7 +229,7 @@ void Particle::event_advance()
{
// Sample a distance to collision
if (type() == ParticleType::electron || type() == ParticleType::positron) {
collision_distance() = 0.0;
collision_distance() = material() == MATERIAL_VOID ? INFINITY : 0.0;
} else if (macro_xs().total == 0.0) {
collision_distance() = INFINITY;
} else {
@ -861,10 +861,12 @@ void Particle::update_neutron_xs(
// If the cache doesn't match, recalculate micro xs
if (this->E() != micro.last_E || this->sqrtkT() != micro.last_sqrtkT ||
i_sab != micro.index_sab || sab_frac != micro.sab_frac) {
i_sab != micro.index_sab || sab_frac != micro.sab_frac ||
ncrystal_xs != micro.ncrystal_xs) {
data::nuclides[i_nuclide]->calculate_xs(i_sab, i_grid, sab_frac, *this);
// If NCrystal is being used, update micro cross section cache
micro.ncrystal_xs = ncrystal_xs;
if (ncrystal_xs >= 0.0) {
data::nuclides[i_nuclide]->calculate_elastic_xs(*this);
ncrystal_update_micro(ncrystal_xs, micro);

2
src/physics.cpp
View file

@ -851,7 +851,7 @@ Direction sample_target_velocity(const Nuclide& nuc, double E, Direction u,
// otherwise, use free gas model
} else {
if (E >= FREE_GAS_THRESHOLD * kT && nuc.awr_ > 1.0) {
if (E >= settings::free_gas_threshold * kT && nuc.awr_ > 1.0) {
return {};
} else {
sampling_method = ResScatMethod::cxs;

373
src/random_ray/flat_source_domain.cpp
View file

@ -53,24 +53,6 @@ FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_)
// Initialize source regions.
bool is_linear = RandomRay::source_shape_ != RandomRaySourceShape::FLAT;
source_regions_ = SourceRegionContainer(negroups_, is_linear);
source_regions_.assign(
base_source_regions, SourceRegion(negroups_, is_linear));
// Initialize materials
int64_t source_region_id = 0;
for (int i = 0; i < model::cells.size(); i++) {
Cell& cell = *model::cells[i];
if (cell.type_ == Fill::MATERIAL) {
for (int j = 0; j < cell.n_instances(); j++) {
source_regions_.material(source_region_id++) = cell.material(j);
}
}
}
// Sanity check
if (source_region_id != base_source_regions) {
fatal_error("Unexpected number of source regions");
}
// Initialize tally volumes
if (volume_normalized_flux_tallies_) {
@ -118,34 +100,24 @@ void FlatSourceDomain::accumulate_iteration_flux()
}
}
// Compute new estimate of scattering + fission sources in each source region
// based on the flux estimate from the previous iteration.
void FlatSourceDomain::update_neutron_source(double k_eff)
void FlatSourceDomain::update_single_neutron_source(SourceRegionHandle& srh)
{
simulation::time_update_src.start();
double inverse_k_eff = 1.0 / k_eff;
// Reset all source regions to zero (important for void regions)
#pragma omp parallel for
for (int64_t se = 0; se < n_source_elements(); se++) {
source_regions_.source(se) = 0.0;
// Reset all source regions to zero (important for void regions)
for (int g = 0; g < negroups_; g++) {
srh.source(g) = 0.0;
}
// Add scattering + fission source
#pragma omp parallel for
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
int material = source_regions_.material(sr);
if (material == MATERIAL_VOID) {
continue;
}
int material = srh.material();
if (material != MATERIAL_VOID) {
double inverse_k_eff = 1.0 / k_eff_;
for (int g_out = 0; g_out < negroups_; g_out++) {
double sigma_t = sigma_t_[material * negroups_ + g_out];
double scatter_source = 0.0;
double fission_source = 0.0;
for (int g_in = 0; g_in < negroups_; g_in++) {
double scalar_flux = source_regions_.scalar_flux_old(sr, g_in);
double scalar_flux = srh.scalar_flux_old(g_in);
double sigma_s =
sigma_s_[material * negroups_ * negroups_ + g_out * negroups_ + g_in];
double nu_sigma_f = nu_sigma_f_[material * negroups_ + g_in];
@ -154,18 +126,30 @@ void FlatSourceDomain::update_neutron_source(double k_eff)
scatter_source += sigma_s * scalar_flux;
fission_source += nu_sigma_f * scalar_flux * chi;
}
source_regions_.source(sr, g_out) =
srh.source(g_out) =
(scatter_source + fission_source * inverse_k_eff) / sigma_t;
}
}
// Add external source if in fixed source mode
if (settings::run_mode == RunMode::FIXED_SOURCE) {
#pragma omp parallel for
for (int64_t se = 0; se < n_source_elements(); se++) {
source_regions_.source(se) += source_regions_.external_source(se);
for (int g = 0; g < negroups_; g++) {
srh.source(g) += srh.external_source(g);
}
}
}
// Compute new estimate of scattering + fission sources in each source region
// based on the flux estimate from the previous iteration.
void FlatSourceDomain::update_all_neutron_sources()
{
simulation::time_update_src.start();
#pragma omp parallel for
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
SourceRegionHandle srh = source_regions_.get_source_region_handle(sr);
update_single_neutron_source(srh);
}
simulation::time_update_src.stop();
}
@ -320,7 +304,7 @@ int64_t FlatSourceDomain::add_source_to_scalar_flux()
// Generates new estimate of k_eff based on the differences between this
// iteration's estimate of the scalar flux and the last iteration's estimate.
double FlatSourceDomain::compute_k_eff(double k_eff_old) const
void FlatSourceDomain::compute_k_eff()
{
double fission_rate_old = 0;
double fission_rate_new = 0;
@ -365,7 +349,7 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const
p[sr] = sr_fission_source_new;
}
double k_eff_new = k_eff_old * (fission_rate_new / fission_rate_old);
double k_eff_new = k_eff_ * (fission_rate_new / fission_rate_old);
double H = 0.0;
// defining an inverse sum for better performance
@ -385,7 +369,7 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const
// Adds entropy value to shared entropy vector in openmc namespace.
simulation::entropy.push_back(H);
return k_eff_new;
k_eff_ = k_eff_new;
}
// This function is responsible for generating a mapping between random
@ -652,7 +636,6 @@ void FlatSourceDomain::random_ray_tally()
"random ray mode.");
break;
}
// Apply score to the appropriate tally bin
Tally& tally {*model::tallies[task.tally_idx]};
#pragma omp atomic
@ -726,21 +709,21 @@ void FlatSourceDomain::output_to_vtk() const
print_plot();
// Outer loop over plots
for (int p = 0; p < model::plots.size(); p++) {
for (int plt = 0; plt < model::plots.size(); plt++) {
// Get handle to OpenMC plot object and extract params
Plot* openmc_plot = dynamic_cast<Plot*>(model::plots[p].get());
Plot* openmc_plot = dynamic_cast<Plot*>(model::plots[plt].get());
// Random ray plots only support voxel plots
if (!openmc_plot) {
warning(fmt::format("Plot {} is invalid plot type -- only voxel plotting "
"is allowed in random ray mode.",
p));
plt));
continue;
} else if (openmc_plot->type_ != Plot::PlotType::voxel) {
warning(fmt::format("Plot {} is invalid plot type -- only voxel plotting "
"is allowed in random ray mode.",
p));
plt));
continue;
}
@ -794,23 +777,11 @@ void FlatSourceDomain::output_to_vtk() const
continue;
}
int i_cell = p.lowest_coord().cell();
int64_t sr = source_region_offsets_[i_cell] + p.cell_instance();
if (RandomRay::mesh_subdivision_enabled_) {
int mesh_idx = base_source_regions_.mesh(sr);
int mesh_bin;
if (mesh_idx == C_NONE) {
mesh_bin = 0;
} else {
mesh_bin = model::meshes[mesh_idx]->get_bin(p.r());
}
SourceRegionKey sr_key {sr, mesh_bin};
auto it = source_region_map_.find(sr_key);
if (it != source_region_map_.end()) {
sr = it->second;
} else {
sr = -1;
}
SourceRegionKey sr_key = lookup_source_region_key(p);
int64_t sr = -1;
auto it = source_region_map_.find(sr_key);
if (it != source_region_map_.end()) {
sr = it->second;
}
voxel_indices[z * Ny * Nx + y * Nx + x] = sr;
@ -967,13 +938,17 @@ void FlatSourceDomain::output_to_vtk() const
}
void FlatSourceDomain::apply_external_source_to_source_region(
Discrete* discrete, double strength_factor, SourceRegionHandle& srh)
int src_idx, SourceRegionHandle& srh)
{
srh.external_source_present() = 1;
auto s = model::external_sources[src_idx].get();
auto is = dynamic_cast<IndependentSource*>(s);
auto discrete = dynamic_cast<Discrete*>(is->energy());
double strength_factor = is->strength();
const auto& discrete_energies = discrete->x();
const auto& discrete_probs = discrete->prob();
srh.external_source_present() = 1;
for (int i = 0; i < discrete_energies.size(); i++) {
int g = data::mg.get_group_index(discrete_energies[i]);
srh.external_source(g) += discrete_probs[i] * strength_factor;
@ -981,8 +956,7 @@ void FlatSourceDomain::apply_external_source_to_source_region(
}
void FlatSourceDomain::apply_external_source_to_cell_instances(int32_t i_cell,
Discrete* discrete, double strength_factor, int target_material_id,
const vector<int32_t>& instances)
int src_idx, int target_material_id, const vector<int32_t>& instances)
{
Cell& cell = *model::cells[i_cell];
@ -1000,16 +974,13 @@ void FlatSourceDomain::apply_external_source_to_cell_instances(int32_t i_cell,
if (target_material_id == C_NONE ||
cell_material_id == target_material_id) {
int64_t source_region = source_region_offsets_[i_cell] + j;
SourceRegionHandle srh =
source_regions_.get_source_region_handle(source_region);
apply_external_source_to_source_region(discrete, strength_factor, srh);
external_volumetric_source_map_[source_region].push_back(src_idx);
}
}
}
void FlatSourceDomain::apply_external_source_to_cell_and_children(
int32_t i_cell, Discrete* discrete, double strength_factor,
int32_t target_material_id)
int32_t i_cell, int src_idx, int32_t target_material_id)
{
Cell& cell = *model::cells[i_cell];
@ -1017,14 +988,14 @@ void FlatSourceDomain::apply_external_source_to_cell_and_children(
vector<int> instances(cell.n_instances());
std::iota(instances.begin(), instances.end(), 0);
apply_external_source_to_cell_instances(
i_cell, discrete, strength_factor, target_material_id, instances);
i_cell, src_idx, target_material_id, instances);
} else if (target_material_id == C_NONE) {
std::unordered_map<int32_t, vector<int32_t>> cell_instance_list =
cell.get_contained_cells(0, nullptr);
for (const auto& pair : cell_instance_list) {
int32_t i_child_cell = pair.first;
apply_external_source_to_cell_instances(i_child_cell, discrete,
strength_factor, target_material_id, pair.second);
apply_external_source_to_cell_instances(
i_child_cell, src_idx, target_material_id, pair.second);
}
}
}
@ -1070,36 +1041,17 @@ void FlatSourceDomain::convert_external_sources()
"point source at {}",
sp->r()));
}
int i_cell = gs.lowest_coord().cell();
int64_t sr = source_region_offsets_[i_cell] + gs.cell_instance();
SourceRegionKey key = lookup_source_region_key(gs);
if (RandomRay::mesh_subdivision_enabled_) {
// If mesh subdivision is enabled, we need to determine which subdivided
// mesh bin the point source coordinate is in as well
int mesh_idx = source_regions_.mesh(sr);
int mesh_bin;
if (mesh_idx == C_NONE) {
mesh_bin = 0;
} else {
mesh_bin = model::meshes[mesh_idx]->get_bin(gs.r());
}
// With the source region and mesh bin known, we can use the
// accompanying SourceRegionKey as a key into a map that stores the
// corresponding external source index for the point source. Notably, we
// do not actually apply the external source to any source regions here,
// as if mesh subdivision is enabled, they haven't actually been
// discovered & initilized yet. When discovered, they will read from the
// point_source_map to determine if there are any point source terms
// that should be applied.
SourceRegionKey key {sr, mesh_bin};
point_source_map_[key] = es;
} else {
// If we are not using mesh subdivision, we can apply the external
// source directly to the source region as we do for volumetric domain
// constraint sources.
SourceRegionHandle srh = source_regions_.get_source_region_handle(sr);
apply_external_source_to_source_region(energy, strength_factor, srh);
}
// With the source region and mesh bin known, we can use the
// accompanying SourceRegionKey as a key into a map that stores the
// corresponding external source index for the point source. Notably, we
// do not actually apply the external source to any source regions here,
// as if mesh subdivision is enabled, they haven't actually been
// discovered & initilized yet. When discovered, they will read from the
// external_source_map to determine if there are any external source
// terms that should be applied.
external_point_source_map_[key].push_back(es);
} else {
// If not a point source, then use the volumetric domain constraints to
@ -1107,42 +1059,25 @@ void FlatSourceDomain::convert_external_sources()
if (is->domain_type() == Source::DomainType::MATERIAL) {
for (int32_t material_id : domain_ids) {
for (int i_cell = 0; i_cell < model::cells.size(); i_cell++) {
apply_external_source_to_cell_and_children(
i_cell, energy, strength_factor, material_id);
apply_external_source_to_cell_and_children(i_cell, es, material_id);
}
}
} else if (is->domain_type() == Source::DomainType::CELL) {
for (int32_t cell_id : domain_ids) {
int32_t i_cell = model::cell_map[cell_id];
apply_external_source_to_cell_and_children(
i_cell, energy, strength_factor, C_NONE);
apply_external_source_to_cell_and_children(i_cell, es, C_NONE);
}
} else if (is->domain_type() == Source::DomainType::UNIVERSE) {
for (int32_t universe_id : domain_ids) {
int32_t i_universe = model::universe_map[universe_id];
Universe& universe = *model::universes[i_universe];
for (int32_t i_cell : universe.cells_) {
apply_external_source_to_cell_and_children(
i_cell, energy, strength_factor, C_NONE);
apply_external_source_to_cell_and_children(i_cell, es, C_NONE);
}
}
}
}
} // End loop over external sources
// Divide the fixed source term by sigma t (to save time when applying each
// iteration)
#pragma omp parallel for
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
int material = source_regions_.material(sr);
if (material == MATERIAL_VOID) {
continue;
}
for (int g = 0; g < negroups_; g++) {
double sigma_t = sigma_t_[material * negroups_ + g];
source_regions_.external_source(sr, g) /= sigma_t;
}
}
}
void FlatSourceDomain::flux_swap()
@ -1159,13 +1094,23 @@ void FlatSourceDomain::flatten_xs()
const int a = 0;
n_materials_ = data::mg.macro_xs_.size();
for (auto& m : data::mg.macro_xs_) {
for (int i = 0; i < n_materials_; i++) {
auto& m = data::mg.macro_xs_[i];
for (int g_out = 0; g_out < negroups_; g_out++) {
if (m.exists_in_model) {
double sigma_t =
m.get_xs(MgxsType::TOTAL, g_out, NULL, NULL, NULL, t, a);
sigma_t_.push_back(sigma_t);
if (sigma_t < MINIMUM_MACRO_XS) {
Material* mat = model::materials[i].get();
warning(fmt::format(
"Material \"{}\" (id: {}) has a group {} total cross section "
"({:.3e}) below the minimum threshold "
"({:.3e}). Material will be treated as pure void.",
mat->name(), mat->id(), g_out, sigma_t, MINIMUM_MACRO_XS));
}
double nu_sigma_f =
m.get_xs(MgxsType::NU_FISSION, g_out, NULL, NULL, NULL, t, a);
nu_sigma_f_.push_back(nu_sigma_f);
@ -1206,7 +1151,7 @@ void FlatSourceDomain::flatten_xs()
}
}
void FlatSourceDomain::set_adjoint_sources(const vector<double>& forward_flux)
void FlatSourceDomain::set_adjoint_sources()
{
// Set the adjoint external source to 1/forward_flux. If the forward flux is
// negative, zero, or extremely close to zero, set the adjoint source to zero,
@ -1220,7 +1165,7 @@ void FlatSourceDomain::set_adjoint_sources(const vector<double>& forward_flux)
double max_flux = 0.0;
#pragma omp parallel for reduction(max : max_flux)
for (int64_t se = 0; se < n_source_elements(); se++) {
double flux = forward_flux[se];
double flux = source_regions_.scalar_flux_final(se);
if (flux > max_flux) {
max_flux = flux;
}
@ -1230,7 +1175,7 @@ void FlatSourceDomain::set_adjoint_sources(const vector<double>& forward_flux)
#pragma omp parallel for
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
for (int g = 0; g < negroups_; g++) {
double flux = forward_flux[sr * negroups_ + g];
double flux = source_regions_.scalar_flux_final(sr, g);
if (flux <= ZERO_FLUX_CUTOFF * max_flux) {
source_regions_.external_source(sr, g) = 0.0;
} else {
@ -1239,6 +1184,7 @@ void FlatSourceDomain::set_adjoint_sources(const vector<double>& forward_flux)
if (flux > 0.0) {
source_regions_.external_source_present(sr) = 1;
}
source_regions_.scalar_flux_final(sr, g) = 0.0;
}
}
@ -1265,7 +1211,6 @@ void FlatSourceDomain::set_adjoint_sources(const vector<double>& forward_flux)
source_regions_.external_source_present(sr) = 0;
}
}
// Divide the fixed source term by sigma t (to save time when applying each
// iteration)
#pragma omp parallel for
@ -1326,13 +1271,14 @@ void FlatSourceDomain::apply_mesh_to_cell_instances(int32_t i_cell,
if ((target_material_id == C_NONE && !is_target_void) ||
cell_material_id == target_material_id) {
int64_t sr = source_region_offsets_[i_cell] + j;
if (source_regions_.mesh(sr) != C_NONE) {
// print out the source region that is broken:
// Check if the key is already present in the mesh_map_
if (mesh_map_.find(sr) != mesh_map_.end()) {
fatal_error(fmt::format("Source region {} already has mesh idx {} "
"applied, but trying to apply mesh idx {}",
sr, source_regions_.mesh(sr), mesh_idx));
sr, mesh_map_[sr], mesh_idx));
}
source_regions_.mesh(sr) = mesh_idx;
// If the SR has not already been assigned, then we can write to it
mesh_map_[sr] = mesh_idx;
}
}
}
@ -1402,18 +1348,9 @@ void FlatSourceDomain::apply_meshes()
}
}
void FlatSourceDomain::prepare_base_source_regions()
{
std::swap(source_regions_, base_source_regions_);
source_regions_.negroups() = base_source_regions_.negroups();
source_regions_.is_linear() = base_source_regions_.is_linear();
}
SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
int64_t sr, int mesh_bin, Position r, double dist, Direction u)
SourceRegionKey sr_key, Position r, Direction u)
{
SourceRegionKey sr_key {sr, mesh_bin};
// Case 1: Check if the source region key is already present in the permanent
// map. This is the most common condition, as any source region visited in a
// previous power iteration will already be present in the permanent map. If
@ -1475,9 +1412,8 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
gs.r() = r + TINY_BIT * u;
gs.u() = {1.0, 0.0, 0.0};
exhaustive_find_cell(gs);
int gs_i_cell = gs.lowest_coord().cell();
int64_t sr_found = source_region_offsets_[gs_i_cell] + gs.cell_instance();
if (sr_found != sr) {
int64_t sr_found = lookup_base_source_region_idx(gs);
if (sr_found != sr_key.base_source_region_id) {
discovered_source_regions_.unlock(sr_key);
SourceRegionHandle handle;
handle.is_numerical_fp_artifact_ = true;
@ -1485,9 +1421,9 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
}
// Sanity check on mesh bin
int mesh_idx = base_source_regions_.mesh(sr);
int mesh_idx = lookup_mesh_idx(sr_key.base_source_region_id);
if (mesh_idx == C_NONE) {
if (mesh_bin != 0) {
if (sr_key.mesh_bin != 0) {
discovered_source_regions_.unlock(sr_key);
SourceRegionHandle handle;
handle.is_numerical_fp_artifact_ = true;
@ -1496,7 +1432,7 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
} else {
Mesh* mesh = model::meshes[mesh_idx].get();
int bin_found = mesh->get_bin(r + TINY_BIT * u);
if (bin_found != mesh_bin) {
if (bin_found != sr_key.mesh_bin) {
discovered_source_regions_.unlock(sr_key);
SourceRegionHandle handle;
handle.is_numerical_fp_artifact_ = true;
@ -1508,26 +1444,60 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
// condition only occurs the first time the source region is discovered
// (typically in the first power iteration). In this case, we need to handle
// creation of the new source region and its storage into the parallel map.
// The new source region is created by copying the base source region, so as
// to inherit material, external source, and some flux properties etc. We
// also pass the base source region id to allow the new source region to
// know which base source region it is derived from.
SourceRegion* sr_ptr = discovered_source_regions_.emplace(
sr_key, {base_source_regions_.get_source_region_handle(sr), sr});
discovered_source_regions_.unlock(sr_key);
// Additionally, we need to determine the source region's material, initialize
// the starting scalar flux guess, and apply any known external sources.
// Call the basic constructor for the source region and store in the parallel
// map.
bool is_linear = RandomRay::source_shape_ != RandomRaySourceShape::FLAT;
SourceRegion* sr_ptr =
discovered_source_regions_.emplace(sr_key, {negroups_, is_linear});
SourceRegionHandle handle {*sr_ptr};
// Check if the new source region contains a point source and apply it if so
auto it2 = point_source_map_.find(sr_key);
if (it2 != point_source_map_.end()) {
int es = it2->second;
auto s = model::external_sources[es].get();
auto is = dynamic_cast<IndependentSource*>(s);
auto energy = dynamic_cast<Discrete*>(is->energy());
double strength_factor = is->strength();
apply_external_source_to_source_region(energy, strength_factor, handle);
int material = handle.material();
if (material != MATERIAL_VOID) {
// Determine the material
int gs_i_cell = gs.lowest_coord().cell();
Cell& cell = *model::cells[gs_i_cell];
int material = cell.material(gs.cell_instance());
// If material total XS is extremely low, just set it to void to avoid
// problems with 1/Sigma_t
for (int g = 0; g < negroups_; g++) {
double sigma_t = sigma_t_[material * negroups_ + g];
if (sigma_t < MINIMUM_MACRO_XS) {
material = MATERIAL_VOID;
break;
}
}
handle.material() = material;
// Store the mesh index (if any) assigned to this source region
handle.mesh() = mesh_idx;
if (settings::run_mode == RunMode::FIXED_SOURCE) {
// Determine if there are any volumetric sources, and apply them.
// Volumetric sources are specifc only to the base SR idx.
auto it_vol =
external_volumetric_source_map_.find(sr_key.base_source_region_id);
if (it_vol != external_volumetric_source_map_.end()) {
const vector<int>& vol_sources = it_vol->second;
for (int src_idx : vol_sources) {
apply_external_source_to_source_region(src_idx, handle);
}
}
// Determine if there are any point sources, and apply them.
// Point sources are specific to the source region key.
auto it_point = external_point_source_map_.find(sr_key);
if (it_point != external_point_source_map_.end()) {
const vector<int>& point_sources = it_point->second;
for (int src_idx : point_sources) {
apply_external_source_to_source_region(src_idx, handle);
}
}
// Divide external source term by sigma_t
if (material != C_NONE) {
for (int g = 0; g < negroups_; g++) {
double sigma_t = sigma_t_[material * negroups_ + g];
handle.external_source(g) /= sigma_t;
@ -1535,6 +1505,21 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
}
}
// Compute the combined source term
update_single_neutron_source(handle);
// Unlock the parallel map. Note: we may be tempted to release
// this lock earlier, and then just use the source region's lock to protect
// the flux/source initialization stages above. However, the rest of the code
// only protects updates to the new flux and volume fields, and assumes that
// the source is constant for the duration of transport. Thus, using just the
// source region's lock by itself would result in other threads potentially
// reading from the source before it is computed, as they won't use the lock
// when only reading from the SR's source. It would be expensive to protect
// those operations, whereas generating the SR is only done once, so we just
// hold the map's bucket lock until the source region is fully initialized.
discovered_source_regions_.unlock(sr_key);
return handle;
}
@ -1620,4 +1605,52 @@ void FlatSourceDomain::apply_transport_stabilization()
}
}
// Determines the base source region index (i.e., a material filled cell
// instance) that corresponds to a particular location in the geometry. Requires
// that the "gs" object passed in has already been initialized and has called
// find_cell etc.
int64_t FlatSourceDomain::lookup_base_source_region_idx(
const GeometryState& gs) const
{
int i_cell = gs.lowest_coord().cell();
int64_t sr = source_region_offsets_[i_cell] + gs.cell_instance();
return sr;
}
// Determines the index of the mesh (if any) that has been applied
// to a particular base source region index.
int FlatSourceDomain::lookup_mesh_idx(int64_t sr) const
{
int mesh_idx = C_NONE;
auto mesh_it = mesh_map_.find(sr);
if (mesh_it != mesh_map_.end()) {
mesh_idx = mesh_it->second;
}
return mesh_idx;
}
// Determines the source region key that corresponds to a particular location in
// the geometry. This takes into account both the base source region index as
// well as the mesh bin if a mesh is applied to this source region for
// subdivision.
SourceRegionKey FlatSourceDomain::lookup_source_region_key(
const GeometryState& gs) const
{
int64_t sr = lookup_base_source_region_idx(gs);
int64_t mesh_bin = lookup_mesh_bin(sr, gs.r());
return SourceRegionKey {sr, mesh_bin};
}
// Determines the mesh bin that corresponds to a particular base source region
// index and position.
int64_t FlatSourceDomain::lookup_mesh_bin(int64_t sr, Position r) const
{
int mesh_idx = lookup_mesh_idx(sr);
int mesh_bin = 0;
if (mesh_idx != C_NONE) {
mesh_bin = model::meshes[mesh_idx]->get_bin(r);
}
return mesh_bin;
}
} // namespace openmc

47
src/random_ray/linear_source_domain.cpp
View file

@ -34,25 +34,18 @@ void LinearSourceDomain::batch_reset()
}
}
void LinearSourceDomain::update_neutron_source(double k_eff)
void LinearSourceDomain::update_single_neutron_source(SourceRegionHandle& srh)
{
simulation::time_update_src.start();
double inverse_k_eff = 1.0 / k_eff;
// Reset all source regions to zero (important for void regions)
#pragma omp parallel for
for (int64_t se = 0; se < n_source_elements(); se++) {
source_regions_.source(se) = 0.0;
// Reset all source regions to zero (important for void regions)
for (int g = 0; g < negroups_; g++) {
srh.source(g) = 0.0;
}
#pragma omp parallel for
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
int material = source_regions_.material(sr);
if (material == MATERIAL_VOID) {
continue;
}
MomentMatrix invM = source_regions_.mom_matrix(sr).inverse();
// Add scattering + fission source
int material = srh.material();
if (material != MATERIAL_VOID) {
double inverse_k_eff = 1.0 / k_eff_;
MomentMatrix invM = srh.mom_matrix().inverse();
for (int g_out = 0; g_out < negroups_; g_out++) {
double sigma_t = sigma_t_[material * negroups_ + g_out];
@ -64,8 +57,8 @@ void LinearSourceDomain::update_neutron_source(double k_eff)
for (int g_in = 0; g_in < negroups_; g_in++) {
// Handles for the flat and linear components of the flux
double flux_flat = source_regions_.scalar_flux_old(sr, g_in);
MomentArray flux_linear = source_regions_.flux_moments_old(sr, g_in);
double flux_flat = srh.scalar_flux_old(g_in);
MomentArray flux_linear = srh.flux_moments_old(g_in);
// Handles for cross sections
double sigma_s =
@ -81,7 +74,7 @@ void LinearSourceDomain::update_neutron_source(double k_eff)
}
// Compute the flat source term
source_regions_.source(sr, g_out) =
srh.source(g_out) =
(scatter_flat + fission_flat * inverse_k_eff) / sigma_t;
// Compute the linear source terms. In the first 10 iterations when the
@ -91,25 +84,21 @@ void LinearSourceDomain::update_neutron_source(double k_eff)
// very small/noisy or have poorly developed spatial moments, so we zero
// the source gradients (effectively making this a flat source region
// temporarily), so as to improve stability.
if (simulation::current_batch > 10 &&
source_regions_.source(sr, g_out) >= 0.0) {
source_regions_.source_gradients(sr, g_out) =
if (simulation::current_batch > 10 && srh.source(g_out) >= 0.0) {
srh.source_gradients(g_out) =
invM * ((scatter_linear + fission_linear * inverse_k_eff) / sigma_t);
} else {
source_regions_.source_gradients(sr, g_out) = {0.0, 0.0, 0.0};
srh.source_gradients(g_out) = {0.0, 0.0, 0.0};
}
}
}
// Add external source if in fixed source mode
if (settings::run_mode == RunMode::FIXED_SOURCE) {
// Add external source to flat source term if in fixed source mode
#pragma omp parallel for
for (int64_t se = 0; se < n_source_elements(); se++) {
source_regions_.source(se) += source_regions_.external_source(se);
for (int g = 0; g < negroups_; g++) {
srh.source(g) += srh.external_source(g);
}
}
simulation::time_update_src.stop();
}
void LinearSourceDomain::normalize_scalar_flux_and_volumes(

121
src/random_ray/random_ray.cpp
View file

@ -237,7 +237,6 @@ double RandomRay::distance_inactive_;
double RandomRay::distance_active_;
unique_ptr<Source> RandomRay::ray_source_;
RandomRaySourceShape RandomRay::source_shape_ {RandomRaySourceShape::FLAT};
bool RandomRay::mesh_subdivision_enabled_ {false};
RandomRaySampleMethod RandomRay::sample_method_ {RandomRaySampleMethod::PRNG};
RandomRay::RandomRay()
@ -336,71 +335,60 @@ void RandomRay::event_advance_ray()
void RandomRay::attenuate_flux(double distance, bool is_active, double offset)
{
// Determine source region index etc.
int i_cell = lowest_coord().cell();
// The base source region is the spatial region index
int64_t sr = domain_->source_region_offsets_[i_cell] + cell_instance();
// Lookup base source region index
int64_t sr = domain_->lookup_base_source_region_idx(*this);
// Perform ray tracing across mesh
if (mesh_subdivision_enabled_) {
// Determine the mesh index for the base source region, if any
int mesh_idx = domain_->base_source_regions_.mesh(sr);
// Determine the mesh index for the base source region, if any
int mesh_idx = domain_->lookup_mesh_idx(sr);
if (mesh_idx == C_NONE) {
// If there's no mesh being applied to this cell, then
// we just attenuate the flux as normal, and set
// the mesh bin to 0
attenuate_flux_inner(distance, is_active, sr, 0, r());
} else {
// If there is a mesh being applied to this cell, then
// we loop over all the bin crossings and attenuate
// separately.
Mesh* mesh = model::meshes[mesh_idx].get();
// We adjust the start and end positions of the ray slightly
// to accomodate for floating point precision issues that tend
// to occur at mesh boundaries that overlap with geometry lattice
// boundaries.
Position start = r() + (offset + TINY_BIT) * u();
Position end = start + (distance - 2.0 * TINY_BIT) * u();
double reduced_distance = (end - start).norm();
// Ray trace through the mesh and record bins and lengths
mesh_bins_.resize(0);
mesh_fractional_lengths_.resize(0);
mesh->bins_crossed(start, end, u(), mesh_bins_, mesh_fractional_lengths_);
// Loop over all mesh bins and attenuate flux
for (int b = 0; b < mesh_bins_.size(); b++) {
double physical_length = reduced_distance * mesh_fractional_lengths_[b];
attenuate_flux_inner(
physical_length, is_active, sr, mesh_bins_[b], start);
start += physical_length * u();
}
}
if (mesh_idx == C_NONE) {
// If there's no mesh being applied to this cell, then
// we just attenuate the flux as normal, and set
// the mesh bin to 0
attenuate_flux_inner(distance, is_active, sr, 0, r());
} else {
attenuate_flux_inner(distance, is_active, sr, C_NONE, r());
// If there is a mesh being applied to this cell, then
// we loop over all the bin crossings and attenuate
// separately.
Mesh* mesh = model::meshes[mesh_idx].get();
// We adjust the start and end positions of the ray slightly
// to accomodate for floating point precision issues that tend
// to occur at mesh boundaries that overlap with geometry lattice
// boundaries.
Position start = r() + (offset + TINY_BIT) * u();
Position end = start + (distance - 2.0 * TINY_BIT) * u();
double reduced_distance = (end - start).norm();
// Ray trace through the mesh and record bins and lengths
mesh_bins_.resize(0);
mesh_fractional_lengths_.resize(0);
mesh->bins_crossed(start, end, u(), mesh_bins_, mesh_fractional_lengths_);
// Loop over all mesh bins and attenuate flux
for (int b = 0; b < mesh_bins_.size(); b++) {
double physical_length = reduced_distance * mesh_fractional_lengths_[b];
attenuate_flux_inner(
physical_length, is_active, sr, mesh_bins_[b], start);
start += physical_length * u();
}
}
}
void RandomRay::attenuate_flux_inner(
double distance, bool is_active, int64_t sr, int mesh_bin, Position r)
{
SourceRegionKey sr_key {sr, mesh_bin};
SourceRegionHandle srh;
if (mesh_subdivision_enabled_) {
srh = domain_->get_subdivided_source_region_handle(
sr, mesh_bin, r, distance, u());
if (srh.is_numerical_fp_artifact_) {
return;
}
} else {
srh = domain_->source_regions_.get_source_region_handle(sr);
srh = domain_->get_subdivided_source_region_handle(sr_key, r, u());
if (srh.is_numerical_fp_artifact_) {
return;
}
switch (source_shape_) {
case RandomRaySourceShape::FLAT:
if (this->material() == MATERIAL_VOID) {
if (srh.material() == MATERIAL_VOID) {
attenuate_flux_flat_source_void(srh, distance, is_active, r);
} else {
attenuate_flux_flat_source(srh, distance, is_active, r);
@ -408,7 +396,7 @@ void RandomRay::attenuate_flux_inner(
break;
case RandomRaySourceShape::LINEAR:
case RandomRaySourceShape::LINEAR_XY:
if (this->material() == MATERIAL_VOID) {
if (srh.material() == MATERIAL_VOID) {
attenuate_flux_linear_source_void(srh, distance, is_active, r);
} else {
attenuate_flux_linear_source(srh, distance, is_active, r);
@ -439,7 +427,7 @@ void RandomRay::attenuate_flux_flat_source(
n_event()++;
// Get material
int material = this->material();
int material = srh.material();
// MOC incoming flux attenuation + source contribution/attenuation equation
for (int g = 0; g < negroups_; g++) {
@ -490,7 +478,7 @@ void RandomRay::attenuate_flux_flat_source_void(
// The number of geometric intersections is counted for reporting purposes
n_event()++;
int material = this->material();
int material = srh.material();
// If ray is in the active phase (not in dead zone), make contributions to
// source region bookkeeping
@ -537,7 +525,7 @@ void RandomRay::attenuate_flux_linear_source(
// The number of geometric intersections is counted for reporting purposes
n_event()++;
int material = this->material();
int material = srh.material();
Position& centroid = srh.centroid();
Position midpoint = r + u() * (distance / 2.0);
@ -810,27 +798,12 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain)
cell_born() = lowest_coord().cell();
}
SourceRegionKey sr_key = domain_->lookup_source_region_key(*this);
SourceRegionHandle srh =
domain_->get_subdivided_source_region_handle(sr_key, r(), u());
// Initialize ray's starting angular flux to starting location's isotropic
// source
int i_cell = lowest_coord().cell();
int64_t sr = domain_->source_region_offsets_[i_cell] + cell_instance();
SourceRegionHandle srh;
if (mesh_subdivision_enabled_) {
int mesh_idx = domain_->base_source_regions_.mesh(sr);
int mesh_bin;
if (mesh_idx == C_NONE) {
mesh_bin = 0;
} else {
Mesh* mesh = model::meshes[mesh_idx].get();
mesh_bin = mesh->get_bin(r());
}
srh =
domain_->get_subdivided_source_region_handle(sr, mesh_bin, r(), 0.0, u());
} else {
srh = domain_->source_regions_.get_source_region_handle(sr);
}
if (!srh.is_numerical_fp_artifact_) {
for (int g = 0; g < negroups_; g++) {
angular_flux_[g] = srh.source(g);

190
src/random_ray/random_ray_simulation.cpp
View file

@ -47,97 +47,82 @@ void openmc_run_random_ray()
if (mpi::master)
validate_random_ray_inputs();
// Declare forward flux so that it can be saved for later adjoint simulation
vector<double> forward_flux;
SourceRegionContainer forward_source_regions;
SourceRegionContainer forward_base_source_regions;
std::unordered_map<SourceRegionKey, int64_t, SourceRegionKey::HashFunctor>
forward_source_region_map;
// Initialize Random Ray Simulation Object
RandomRaySimulation sim;
{
// Initialize Random Ray Simulation Object
RandomRaySimulation sim;
// Initialize fixed sources, if present
sim.apply_fixed_sources_and_mesh_domains();
// Initialize fixed sources, if present
sim.apply_fixed_sources_and_mesh_domains();
// Begin main simulation timer
simulation::time_total.start();
// Begin main simulation timer
simulation::time_total.start();
// Execute random ray simulation
sim.simulate();
// Execute random ray simulation
sim.simulate();
// End main simulation timer
simulation::time_total.stop();
// End main simulation timer
simulation::time_total.stop();
// Normalize and save the final forward flux
sim.domain()->serialize_final_fluxes(forward_flux);
double source_normalization_factor =
sim.domain()->compute_fixed_source_normalization_factor() /
(settings::n_batches - settings::n_inactive);
// Normalize and save the final forward flux
double source_normalization_factor =
sim.domain()->compute_fixed_source_normalization_factor() /
(settings::n_batches - settings::n_inactive);
#pragma omp parallel for
for (uint64_t i = 0; i < forward_flux.size(); i++) {
forward_flux[i] *= source_normalization_factor;
}
forward_source_regions = sim.domain()->source_regions_;
forward_source_region_map = sim.domain()->source_region_map_;
forward_base_source_regions = sim.domain()->base_source_regions_;
// Finalize OpenMC
openmc_simulation_finalize();
// Output all simulation results
sim.output_simulation_results();
for (uint64_t se = 0; se < sim.domain()->n_source_elements(); se++) {
sim.domain()->source_regions_.scalar_flux_final(se) *=
source_normalization_factor;
}
// Finalize OpenMC
openmc_simulation_finalize();
// Output all simulation results
sim.output_simulation_results();
//////////////////////////////////////////////////////////
// Run adjoint simulation (if enabled)
//////////////////////////////////////////////////////////
if (adjoint_needed) {
reset_timers();
// Configure the domain for adjoint simulation
FlatSourceDomain::adjoint_ = true;
if (mpi::master)
header("ADJOINT FLUX SOLVE", 3);
// Initialize OpenMC general data structures
openmc_simulation_init();
// Initialize Random Ray Simulation Object
RandomRaySimulation adjoint_sim;
// Initialize adjoint fixed sources, if present
adjoint_sim.prepare_fixed_sources_adjoint(forward_flux,
forward_source_regions, forward_base_source_regions,
forward_source_region_map);
// Transpose scattering matrix
adjoint_sim.domain()->transpose_scattering_matrix();
// Swap nu_sigma_f and chi
adjoint_sim.domain()->nu_sigma_f_.swap(adjoint_sim.domain()->chi_);
// Begin main simulation timer
simulation::time_total.start();
// Execute random ray simulation
adjoint_sim.simulate();
// End main simulation timer
simulation::time_total.stop();
// Finalize OpenMC
openmc_simulation_finalize();
// Output all simulation results
adjoint_sim.output_simulation_results();
if (!adjoint_needed) {
return;
}
reset_timers();
// Configure the domain for adjoint simulation
FlatSourceDomain::adjoint_ = true;
if (mpi::master)
header("ADJOINT FLUX SOLVE", 3);
// Initialize OpenMC general data structures
openmc_simulation_init();
sim.domain()->k_eff_ = 1.0;
// Initialize adjoint fixed sources, if present
sim.prepare_fixed_sources_adjoint();
// Transpose scattering matrix
sim.domain()->transpose_scattering_matrix();
// Swap nu_sigma_f and chi
sim.domain()->nu_sigma_f_.swap(sim.domain()->chi_);
// Begin main simulation timer
simulation::time_total.start();
// Execute random ray simulation
sim.simulate();
// End main simulation timer
simulation::time_total.stop();
// Finalize OpenMC
openmc_simulation_finalize();
// Output all simulation results
sim.output_simulation_results();
}
// Enforces restrictions on inputs in random ray mode. While there are
@ -348,7 +333,6 @@ void validate_random_ray_inputs()
// when generating weight windows with FW-CADIS and an overlaid mesh.
///////////////////////////////////////////////////////////////////
if (RandomRay::source_shape_ == RandomRaySourceShape::LINEAR &&
RandomRay::mesh_subdivision_enabled_ &&
variance_reduction::weight_windows.size() > 0) {
warning(
"Linear sources may result in negative fluxes in small source regions "
@ -366,7 +350,6 @@ void openmc_reset_random_ray()
FlatSourceDomain::mesh_domain_map_.clear();
RandomRay::ray_source_.reset();
RandomRay::source_shape_ = RandomRaySourceShape::FLAT;
RandomRay::mesh_subdivision_enabled_ = false;
RandomRay::sample_method_ = RandomRaySampleMethod::PRNG;
}
@ -412,20 +395,11 @@ void RandomRaySimulation::apply_fixed_sources_and_mesh_domains()
}
}
void RandomRaySimulation::prepare_fixed_sources_adjoint(
vector<double>& forward_flux, SourceRegionContainer& forward_source_regions,
SourceRegionContainer& forward_base_source_regions,
std::unordered_map<SourceRegionKey, int64_t, SourceRegionKey::HashFunctor>&
forward_source_region_map)
void RandomRaySimulation::prepare_fixed_sources_adjoint()
{
domain_->source_regions_.adjoint_reset();
if (settings::run_mode == RunMode::FIXED_SOURCE) {
if (RandomRay::mesh_subdivision_enabled_) {
domain_->source_regions_ = forward_source_regions;
domain_->source_region_map_ = forward_source_region_map;
domain_->base_source_regions_ = forward_base_source_regions;
domain_->source_regions_.adjoint_reset();
}
domain_->set_adjoint_sources(forward_flux);
domain_->set_adjoint_sources();
}
}
@ -445,22 +419,18 @@ void RandomRaySimulation::simulate()
simulation::total_weight = 1.0;
// Update source term (scattering + fission)
domain_->update_neutron_source(k_eff_);
domain_->update_all_neutron_sources();
// Reset scalar fluxes, iteration volume tallies, and region hit flags to
// zero
// Reset scalar fluxes, iteration volume tallies, and region hit flags
// to zero
domain_->batch_reset();
// At the beginning of the simulation, if mesh subvivision is in use, we
// At the beginning of the simulation, if mesh subdivision is in use, we
// need to swap the main source region container into the base container,
// as the main source region container will be used to hold the true
// subdivided source regions. The base container will therefore only
// contain the external source region information, the mesh indices,
// material properties, and initial guess values for the flux/source.
if (RandomRay::mesh_subdivision_enabled_ &&
simulation::current_batch == 1 && !FlatSourceDomain::adjoint_) {
domain_->prepare_base_source_regions();
}
// Start timer for transport
simulation::time_transport.start();
@ -476,11 +446,9 @@ void RandomRaySimulation::simulate()
simulation::time_transport.stop();
// If using mesh subdivision, add any newly discovered source regions
// to the main source region container.
if (RandomRay::mesh_subdivision_enabled_) {
domain_->finalize_discovered_source_regions();
}
// Add any newly discovered source regions to the main source region
// container.
domain_->finalize_discovered_source_regions();
// Normalize scalar flux and update volumes
domain_->normalize_scalar_flux_and_volumes(
@ -494,10 +462,10 @@ void RandomRaySimulation::simulate()
if (settings::run_mode == RunMode::EIGENVALUE) {
// Compute random ray k-eff
k_eff_ = domain_->compute_k_eff(k_eff_);
domain_->compute_k_eff();
// Store random ray k-eff into OpenMC's native k-eff variable
global_tally_tracklength = k_eff_;
global_tally_tracklength = domain_->k_eff_;
}
// Execute all tallying tasks, if this is an active batch
@ -507,12 +475,6 @@ void RandomRaySimulation::simulate()
// estimate
domain_->accumulate_iteration_flux();
// Generate mapping between source regions and tallies
if (!domain_->mapped_all_tallies_ &&
!RandomRay::mesh_subdivision_enabled_) {
domain_->convert_source_regions_to_tallies(0);
}
// Use above mapping to contribute FSR flux data to appropriate
// tallies
domain_->random_ray_tally();
@ -522,7 +484,7 @@ void RandomRaySimulation::simulate()
domain_->flux_swap();
// Check for any obvious insabilities/nans/infs
instability_check(n_hits, k_eff_, avg_miss_rate_);
instability_check(n_hits, domain_->k_eff_, avg_miss_rate_);
} // End MPI master work
// Finalize the current batch
@ -571,7 +533,7 @@ void RandomRaySimulation::instability_check(
}
if (k_eff > 10.0 || k_eff < 0.01 || !(std::isfinite(k_eff))) {
fatal_error("Instability detected");
fatal_error(fmt::format("Instability detected: k-eff = {:.5f}", k_eff));
}
}
}

28
src/random_ray/source_region.cpp
View file

@ -48,7 +48,7 @@ SourceRegion::SourceRegion(int negroups, bool is_linear)
}
scalar_flux_new_.assign(negroups, 0.0);
source_.resize(negroups);
source_.assign(negroups, 0.0);
scalar_flux_final_.assign(negroups, 0.0);
tally_task_.resize(negroups);
@ -60,25 +60,6 @@ SourceRegion::SourceRegion(int negroups, bool is_linear)
}
}
SourceRegion::SourceRegion(const SourceRegionHandle& handle, int64_t parent_sr)
: SourceRegion(handle.negroups_, handle.is_linear_)
{
material_ = handle.material();
mesh_ = handle.mesh();
parent_sr_ = parent_sr;
for (int g = 0; g < scalar_flux_new_.size(); g++) {
scalar_flux_old_[g] = handle.scalar_flux_old(g);
source_[g] = handle.source(g);
}
if (settings::run_mode == RunMode::FIXED_SOURCE) {
external_source_present_ = handle.external_source_present();
for (int g = 0; g < scalar_flux_new_.size(); g++) {
external_source_[g] = handle.external_source(g);
}
}
}
//==============================================================================
// SourceRegionContainer implementation
//==============================================================================
@ -259,9 +240,12 @@ void SourceRegionContainer::adjoint_reset()
MomentMatrix {0.0, 0.0, 0.0, 0.0, 0.0, 0.0});
std::fill(mom_matrix_t_.begin(), mom_matrix_t_.end(),
MomentMatrix {0.0, 0.0, 0.0, 0.0, 0.0, 0.0});
std::fill(scalar_flux_old_.begin(), scalar_flux_old_.end(), 0.0);
if (settings::run_mode == RunMode::FIXED_SOURCE) {
std::fill(scalar_flux_old_.begin(), scalar_flux_old_.end(), 0.0);
} else {
std::fill(scalar_flux_old_.begin(), scalar_flux_old_.end(), 1.0);
}
std::fill(scalar_flux_new_.begin(), scalar_flux_new_.end(), 0.0);
std::fill(scalar_flux_final_.begin(), scalar_flux_final_.end(), 0.0);
std::fill(source_.begin(), source_.end(), 0.0f);
std::fill(external_source_.begin(), external_source_.end(), 0.0f);
std::fill(source_gradients_.begin(), source_gradients_.end(),

6
src/settings.cpp
View file

@ -126,6 +126,7 @@ SolverType solver_type {SolverType::MONTE_CARLO};
std::unordered_set<int> sourcepoint_batch;
std::unordered_set<int> statepoint_batch;
double source_rejection_fraction {0.05};
double free_gas_threshold {400.0};
std::unordered_set<int> source_write_surf_id;
int64_t ssw_max_particles;
int64_t ssw_max_files;
@ -346,7 +347,6 @@ void get_run_parameters(pugi::xml_node node_base)
}
FlatSourceDomain::mesh_domain_map_[mesh_id].emplace_back(
type, domain_id);
RandomRay::mesh_subdivision_enabled_ = true;
}
}
}
@ -652,6 +652,10 @@ void read_settings_xml(pugi::xml_node root)
std::stod(get_node_value(root, "source_rejection_fraction"));
}
if (check_for_node(root, "free_gas_threshold")) {
free_gas_threshold = std::stod(get_node_value(root, "free_gas_threshold"));
}
// Survival biasing
if (check_for_node(root, "survival_biasing")) {
survival_biasing = get_node_value_bool(root, "survival_biasing");

21
src/simulation.cpp
View file

@ -674,8 +674,9 @@ void calculate_work()
void initialize_data()
{
// Determine minimum/maximum energy for incident neutron/photon data
data::energy_max = {INFTY, INFTY};
data::energy_min = {0.0, 0.0};
data::energy_max = {INFTY, INFTY, INFTY, INFTY};
data::energy_min = {0.0, 0.0, 0.0, 0.0};
for (const auto& nuc : data::nuclides) {
if (nuc->grid_.size() >= 1) {
int neutron = static_cast<int>(ParticleType::neutron);
@ -703,11 +704,21 @@ void initialize_data()
// than the current minimum/maximum
if (data::ttb_e_grid.size() >= 1) {
int photon = static_cast<int>(ParticleType::photon);
int electron = static_cast<int>(ParticleType::electron);
int positron = static_cast<int>(ParticleType::positron);
int n_e = data::ttb_e_grid.size();
const std::vector<int> charged = {electron, positron};
for (auto t : charged) {
data::energy_min[t] = std::exp(data::ttb_e_grid(1));
data::energy_max[t] = std::exp(data::ttb_e_grid(n_e - 1));
}
data::energy_min[photon] =
std::max(data::energy_min[photon], std::exp(data::ttb_e_grid(1)));
data::energy_max[photon] = std::min(
data::energy_max[photon], std::exp(data::ttb_e_grid(n_e - 1)));
std::max(data::energy_min[photon], data::energy_min[electron]);
data::energy_max[photon] =
std::min(data::energy_max[photon], data::energy_max[electron]);
}
}
}

6
src/source.cpp
View file

@ -290,6 +290,12 @@ IndependentSource::IndependentSource(pugi::xml_node node) : Source(node)
} else if (temp_str == "photon") {
particle_ = ParticleType::photon;
settings::photon_transport = true;
} else if (temp_str == "electron") {
particle_ = ParticleType::electron;
settings::photon_transport = true;
} else if (temp_str == "positron") {
particle_ = ParticleType::positron;
settings::photon_transport = true;
} else {
fatal_error(std::string("Unknown source particle type: ") + temp_str);
}

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

32
tests/regression_tests/electron_heating/inputs_true.dat Normal file
View file

@ -0,0 +1,32 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material id="1">
<density value="1.0" units="g/cc"/>
<nuclide name="H1" ao="2.0"/>
<nuclide name="O16" ao="1.0"/>
</material>
</materials>
<geometry>
<cell id="1" material="1" region="-1" universe="1"/>
<surface id="1" type="sphere" boundary="reflective" coeffs="0.0 0.0 0.0 1.0"/>
</geometry>
<settings>
<run_mode>fixed source</run_mode>
<particles>10000</particles>
<batches>1</batches>
<source type="independent" strength="1.0" particle="electron">
<energy type="discrete">
<parameters>10000000.0 1.0</parameters>
</energy>
</source>
<cutoff>
<energy_photon>1000.0</energy_photon>
</cutoff>
</settings>
<tallies>
<tally id="1">
<scores>heating</scores>
</tally>
</tallies>
</model>

3
tests/regression_tests/electron_heating/results_true.dat Normal file
View file

@ -0,0 +1,3 @@
tally 1:
1.000000E+07
1.000000E+14

40
tests/regression_tests/electron_heating/test.py Normal file
View file

@ -0,0 +1,40 @@
import pytest
import openmc
from tests.testing_harness import PyAPITestHarness
@pytest.fixture
def water_model():
# Define materals and geometry
water = openmc.Material()
water.add_nuclide("H1", 2.0)
water.add_nuclide("O16", 1.0)
water.set_density("g/cc", 1.0)
sphere = openmc.Sphere(r=1.0, boundary_type="reflective")
sph = openmc.Cell(fill=water, region=-sphere)
geometry = openmc.Geometry([sph])
source = openmc.IndependentSource(
energy=openmc.stats.delta_function(10.0e6),
particle="electron"
)
# Define settings
settings = openmc.Settings()
settings.particles = 10000
settings.batches = 1
settings.cutoff = {"energy_photon": 1000.0}
settings.run_mode = "fixed source"
settings.source = source
# Define tallies
tally = openmc.Tally()
tally.scores = ["heating"]
tallies = openmc.Tallies([tally])
return openmc.Model(geometry=geometry, settings=settings, tallies=tallies)
def test_electron_heating_calc(water_model):
harness = PyAPITestHarness("statepoint.1.h5", water_model)
harness.main()

4
tests/regression_tests/random_ray_adjoint_fixed_source/inputs_true.dat
View file

@ -212,8 +212,8 @@
<parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<adjoint>True</adjoint>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
<adjoint>true</adjoint>
<volume_estimator>naive</volume_estimator>
</random_ray>
</settings>

4
tests/regression_tests/random_ray_adjoint_k_eff/inputs_true.dat
View file

@ -85,8 +85,8 @@
<parameters>-1.26 -1.26 -1 1.26 1.26 1</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<adjoint>True</adjoint>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
<adjoint>true</adjoint>
</random_ray>
</settings>
<tallies>

2
tests/regression_tests/random_ray_fixed_source_domain/cell/inputs_true.dat
View file

@ -212,7 +212,7 @@
<parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
</random_ray>
</settings>
<tallies>

2
tests/regression_tests/random_ray_fixed_source_domain/material/inputs_true.dat
View file

@ -212,7 +212,7 @@
<parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
</random_ray>
</settings>
<tallies>

2
tests/regression_tests/random_ray_fixed_source_domain/universe/inputs_true.dat
View file

@ -212,7 +212,7 @@
<parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
</random_ray>
</settings>
<tallies>

2
tests/regression_tests/random_ray_fixed_source_linear/linear/inputs_true.dat
View file

@ -212,7 +212,7 @@
<parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
<source_shape>linear</source_shape>
</random_ray>
</settings>

2
tests/regression_tests/random_ray_fixed_source_linear/linear_xy/inputs_true.dat
View file

@ -212,7 +212,7 @@
<parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
<source_shape>linear_xy</source_shape>
</random_ray>
</settings>

2
tests/regression_tests/random_ray_fixed_source_mesh/flat/inputs_true.dat
View file

@ -212,7 +212,7 @@
<parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
<source_region_meshes>
<mesh id="1">
<domain id="1" type="universe"/>

2
tests/regression_tests/random_ray_fixed_source_mesh/linear/inputs_true.dat
View file

@ -212,7 +212,7 @@
<parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
<source_region_meshes>
<mesh id="1">
<domain id="1" type="universe"/>

2
tests/regression_tests/random_ray_fixed_source_normalization/False/inputs_true.dat
View file

@ -212,7 +212,7 @@
<parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
</space>
</source>
<volume_normalized_flux_tallies>False</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>false</volume_normalized_flux_tallies>
</random_ray>
</settings>
<tallies>

2
tests/regression_tests/random_ray_fixed_source_normalization/True/inputs_true.dat
View file

@ -212,7 +212,7 @@
<parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
</random_ray>
</settings>
<tallies>

2
tests/regression_tests/random_ray_fixed_source_subcritical/flat/inputs_true.dat
View file

@ -115,7 +115,7 @@
<parameters>-1.26 -1.26 -1 1.26 1.26 1</parameters>
</space>
</source>
<volume_normalized_flux_tallies>False</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>false</volume_normalized_flux_tallies>
<source_shape>flat</source_shape>
</random_ray>
</settings>

2
tests/regression_tests/random_ray_fixed_source_subcritical/linear_xy/inputs_true.dat
View file

@ -115,7 +115,7 @@
<parameters>-1.26 -1.26 -1 1.26 1.26 1</parameters>
</space>
</source>
<volume_normalized_flux_tallies>False</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>false</volume_normalized_flux_tallies>
<source_shape>linear_xy</source_shape>
</random_ray>
</settings>

2
tests/regression_tests/random_ray_halton_samples/inputs_true.dat
View file

@ -85,7 +85,7 @@
<parameters>-1.26 -1.26 -1 1.26 1.26 1</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
<sample_method>halton</sample_method>
</random_ray>
</settings>

2
tests/regression_tests/random_ray_k_eff/inputs_true.dat
View file

@ -85,7 +85,7 @@
<parameters>-1.26 -1.26 -1 1.26 1.26 1</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
</random_ray>
</settings>
<tallies>

2
tests/regression_tests/random_ray_k_eff_mesh/inputs_true.dat
View file

@ -85,7 +85,7 @@
<parameters>-1.26 -1.26 -1 1.26 1.26 1</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
<source_region_meshes>
<mesh id="2">
<domain id="7" type="universe"/>

2
tests/regression_tests/random_ray_linear/linear/inputs_true.dat
View file

@ -85,7 +85,7 @@
<parameters>-1.26 -1.26 -1 1.26 1.26 1</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
<source_shape>linear</source_shape>
</random_ray>
</settings>

2
tests/regression_tests/random_ray_linear/linear_xy/inputs_true.dat
View file

@ -85,7 +85,7 @@
<parameters>-1.26 -1.26 -1 1.26 1.26 1</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
<source_shape>linear_xy</source_shape>
</random_ray>
</settings>

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

244
tests/regression_tests/random_ray_low_density/inputs_true.dat Normal file
View file

@ -0,0 +1,244 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<cross_sections>mgxs.h5</cross_sections>
<material id="1" name="source">
<density value="1.0" units="macro"/>
<macroscopic name="source"/>
</material>
<material id="2" name="void">
<density value="1.0" units="macro"/>
<macroscopic name="void"/>
</material>
<material id="3" name="absorber">
<density value="1.0" units="macro"/>
<macroscopic name="absorber"/>
</material>
</materials>
<geometry>
<cell id="1" name="infinite source region" material="1" universe="1"/>
<cell id="2" name="infinite void region" material="2" universe="2"/>
<cell id="3" name="infinite absorber region" material="3" universe="3"/>
<cell id="4" fill="4" universe="5"/>
<cell id="5" name="full domain" fill="5" region="1 -2 3 -4 5 -6" universe="6"/>
<lattice id="4">
<pitch>2.5 2.5 2.5</pitch>
<dimension>12 12 12</dimension>
<lower_left>0.0 0.0 0.0</lower_left>
<universes>
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
1 1 2 2 2 2 2 2 2 2 3 3
1 1 2 2 2 2 2 2 2 2 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
1 1 2 2 2 2 2 2 2 2 3 3
1 1 2 2 2 2 2 2 2 2 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 </universes>
</lattice>
<surface id="1" type="x-plane" boundary="reflective" coeffs="0.0"/>
<surface id="2" type="x-plane" boundary="vacuum" coeffs="30.0"/>
<surface id="3" type="y-plane" boundary="reflective" coeffs="0.0"/>
<surface id="4" type="y-plane" boundary="vacuum" coeffs="30.0"/>
<surface id="5" type="z-plane" boundary="reflective" coeffs="0.0"/>
<surface id="6" type="z-plane" boundary="vacuum" coeffs="30.0"/>
</geometry>
<settings>
<run_mode>fixed source</run_mode>
<particles>90</particles>
<batches>10</batches>
<inactive>5</inactive>
<source type="independent" strength="3.14" particle="neutron">
<energy type="discrete">
<parameters>100.0 1.0</parameters>
</energy>
<constraints>
<domain_type>universe</domain_type>
<domain_ids>1</domain_ids>
</constraints>
</source>
<energy_mode>multi-group</energy_mode>
<random_ray>
<distance_active>500.0</distance_active>
<distance_inactive>100.0</distance_inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
</space>
</source>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
</random_ray>
</settings>
<tallies>
<filter id="3" type="material">
<bins>1</bins>
</filter>
<filter id="2" type="material">
<bins>2</bins>
</filter>
<filter id="1" type="material">
<bins>3</bins>
</filter>
<tally id="3" name="Source Tally">
<filters>3</filters>
<scores>flux</scores>
<estimator>tracklength</estimator>
</tally>
<tally id="2" name="Void Tally">
<filters>2</filters>
<scores>flux</scores>
<estimator>tracklength</estimator>
</tally>
<tally id="1" name="Absorber Tally">
<filters>1</filters>
<scores>flux</scores>
<estimator>tracklength</estimator>
</tally>
</tallies>
</model>

9
tests/regression_tests/random_ray_low_density/results_true.dat Normal file
View file

@ -0,0 +1,9 @@
tally 1:
5.973607E-01
7.155477E-02
tally 2:
3.206216E-02
2.063375E-04
tally 3:
2.096415E-03
8.804963E-07

60
tests/regression_tests/random_ray_low_density/test.py Normal file
View file

@ -0,0 +1,60 @@
import os
import numpy as np
import openmc
from openmc.examples import random_ray_three_region_cube
from tests.testing_harness import TolerantPyAPITestHarness
class MGXSTestHarness(TolerantPyAPITestHarness):
def _cleanup(self):
super()._cleanup()
f = 'mgxs.h5'
if os.path.exists(f):
os.remove(f)
def test_random_ray_low_density():
model = random_ray_three_region_cube()
# Rebuild the MGXS library to have a material with very
# low macroscopic cross sections
ebins = [1e-5, 20.0e6]
groups = openmc.mgxs.EnergyGroups(group_edges=ebins)
void_sigma_a = 4.0e-6
void_sigma_s = 3.0e-4
void_mat_data = openmc.XSdata('void', groups)
void_mat_data.order = 0
void_mat_data.set_total([void_sigma_a + void_sigma_s])
void_mat_data.set_absorption([void_sigma_a])
void_mat_data.set_scatter_matrix(
np.rollaxis(np.array([[[void_sigma_s]]]), 0, 3))
absorber_sigma_a = 0.75
absorber_sigma_s = 0.25
absorber_mat_data = openmc.XSdata('absorber', groups)
absorber_mat_data.order = 0
absorber_mat_data.set_total([absorber_sigma_a + absorber_sigma_s])
absorber_mat_data.set_absorption([absorber_sigma_a])
absorber_mat_data.set_scatter_matrix(
np.rollaxis(np.array([[[absorber_sigma_s]]]), 0, 3))
multiplier = 0.0000001
source_sigma_a = void_sigma_a * multiplier
source_sigma_s = void_sigma_s * multiplier
source_mat_data = openmc.XSdata('source', groups)
source_mat_data.order = 0
source_mat_data.set_total([source_sigma_a + source_sigma_s])
source_mat_data.set_absorption([source_sigma_a])
source_mat_data.set_scatter_matrix(
np.rollaxis(np.array([[[source_sigma_s]]]), 0, 3))
mg_cross_sections_file = openmc.MGXSLibrary(groups)
mg_cross_sections_file.add_xsdatas(
[source_mat_data, void_mat_data, absorber_mat_data])
mg_cross_sections_file.export_to_hdf5()
harness = MGXSTestHarness('statepoint.10.h5', model)
harness.main()

2
tests/regression_tests/random_ray_point_source_locator/inputs_true.dat
View file

@ -211,7 +211,7 @@
<parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
<source_region_meshes>
<mesh id="1">
<domain id="6" type="universe"/>

8
tests/regression_tests/random_ray_point_source_locator/results_true.dat
View file

@ -1,9 +1,9 @@
tally 1:
2.633900E+00
2.948207E+00
2.633923E+00
2.948228E+00
tally 2:
1.440463E-01
3.294032E-03
1.440456E-01
3.293984E-03
tally 3:
9.425207E-03
1.089748E-05

2
tests/regression_tests/random_ray_void/flat/inputs_true.dat
View file

@ -212,7 +212,7 @@
<parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
<source_shape>flat</source_shape>
</random_ray>
</settings>

2
tests/regression_tests/random_ray_void/linear/inputs_true.dat
View file

@ -212,7 +212,7 @@
<parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
<source_shape>linear</source_shape>
</random_ray>
</settings>

2
tests/regression_tests/random_ray_volume_estimator/hybrid/inputs_true.dat
View file

@ -212,7 +212,7 @@
<parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
<volume_estimator>hybrid</volume_estimator>
</random_ray>
</settings>

2
tests/regression_tests/random_ray_volume_estimator/naive/inputs_true.dat
View file

@ -212,7 +212,7 @@
<parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
<volume_estimator>naive</volume_estimator>
</random_ray>
</settings>

2
tests/regression_tests/random_ray_volume_estimator/simulation_averaged/inputs_true.dat
View file

@ -212,7 +212,7 @@
<parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
<volume_estimator>simulation_averaged</volume_estimator>
</random_ray>
</settings>

2
tests/regression_tests/random_ray_volume_estimator_linear/hybrid/inputs_true.dat
View file

@ -212,7 +212,7 @@
<parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
<source_shape>linear</source_shape>
<volume_estimator>hybrid</volume_estimator>
</random_ray>

2
tests/regression_tests/random_ray_volume_estimator_linear/naive/inputs_true.dat
View file

@ -212,7 +212,7 @@
<parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
<source_shape>linear</source_shape>
<volume_estimator>naive</volume_estimator>
</random_ray>

2
tests/regression_tests/random_ray_volume_estimator_linear/simulation_averaged/inputs_true.dat
View file

@ -212,7 +212,7 @@
<parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
<source_shape>linear</source_shape>
<volume_estimator>simulation_averaged</volume_estimator>
</random_ray>

2
tests/regression_tests/weightwindows_fw_cadis/inputs_true.dat
View file

@ -227,7 +227,7 @@
<parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
<volume_estimator>naive</volume_estimator>
</random_ray>
</settings>

2
tests/regression_tests/weightwindows_fw_cadis_mesh/flat/inputs_true.dat
View file

@ -227,7 +227,7 @@
<parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
<source_region_meshes>
<mesh id="1">
<domain id="6" type="universe"/>

2
tests/regression_tests/weightwindows_fw_cadis_mesh/linear/inputs_true.dat
View file

@ -227,7 +227,7 @@
<parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
<source_region_meshes>
<mesh id="1">
<domain id="6" type="universe"/>

1
tests/unit_tests/test_deplete_activation.py
View file

@ -49,6 +49,7 @@ ENERGIES = np.logspace(log10(1e-5), log10(2e7), 100)
("flux", {'energies': ENERGIES, 'reactions': ['(n,gamma)']}, 1e-5),
("flux", {'energies': ENERGIES, 'reactions': ['(n,gamma)'], 'nuclides': ['W186', 'H3']}, 1e-2),
])
@pytest.mark.flaky(reruns=1)
def test_activation(run_in_tmpdir, model, reaction_rate_mode, reaction_rate_opts, tolerance):
# Determine (n.gamma) reaction rate using initial run
sp = model.run()

27
tests/unit_tests/test_settings.py
View file

@ -59,16 +59,28 @@ def test_export_to_xml(run_in_tmpdir):
s.electron_treatment = 'led'
s.write_initial_source = True
s.weight_window_checkpoints = {'surface': True, 'collision': False}
source_region_mesh = openmc.RegularMesh()
source_region_mesh.dimension = [2, 2, 2]
source_region_mesh.lower_left = [-2, -2, -2]
source_region_mesh.upper_right = [2, 2, 2]
root_universe = openmc.Universe()
s.random_ray = {
'distance_inactive': 10.0,
'distance_active': 100.0,
'ray_source': openmc.IndependentSource(
space=openmc.stats.Box((-1., -1., -1.), (1., 1., 1.))
)
),
'source_region_meshes': [(source_region_mesh, [root_universe])],
'volume_estimator': 'hybrid',
'source_shape': 'linear',
'volume_normalized_flux_tallies': True,
'adjoint': False,
'sample_method': 'halton'
}
s.max_particle_events = 100
s.max_secondaries = 1_000_000
s.source_rejection_fraction = 0.01
s.free_gas_threshold = 800.0
# Make sure exporting XML works
s.export_to_xml()
@ -145,5 +157,18 @@ def test_export_to_xml(run_in_tmpdir):
assert s.random_ray['distance_active'] == 100.0
assert s.random_ray['ray_source'].space.lower_left == [-1., -1., -1.]
assert s.random_ray['ray_source'].space.upper_right == [1., 1., 1.]
assert 'source_region_meshes' in s.random_ray
assert len(s.random_ray['source_region_meshes']) == 1
mesh_and_domains = s.random_ray['source_region_meshes'][0]
recovered_mesh = mesh_and_domains[0]
assert recovered_mesh.dimension == (2, 2, 2)
assert recovered_mesh.lower_left == [-2., -2., -2.]
assert recovered_mesh.upper_right == [2., 2., 2.]
assert s.random_ray['volume_estimator'] == 'hybrid'
assert s.random_ray['source_shape'] == 'linear'
assert s.random_ray['volume_normalized_flux_tallies']
assert not s.random_ray['adjoint']
assert s.random_ray['sample_method'] == 'halton'
assert s.max_secondaries == 1_000_000
assert s.source_rejection_fraction == 0.01
assert s.free_gas_threshold == 800.0

65
tests/unit_tests/test_statepoint.py Normal file
View file

@ -0,0 +1,65 @@
import openmc
def test_get_tally_filter_type(run_in_tmpdir):
"""Test various ways of retrieving tallies from a StatePoint object."""
mat = openmc.Material()
mat.add_nuclide("H1", 1.0)
mat.set_density("g/cm3", 10.0)
sphere = openmc.Sphere(r=10.0, boundary_type="vacuum")
cell = openmc.Cell(fill=mat, region=-sphere)
geometry = openmc.Geometry([cell])
settings = openmc.Settings()
settings.particles = 10
settings.batches = 2
settings.run_mode = "fixed source"
reg_mesh = openmc.RegularMesh().from_domain(cell)
tally1 = openmc.Tally(tally_id=1)
mesh_filter = openmc.MeshFilter(reg_mesh)
tally1.filters = [mesh_filter]
tally1.scores = ["flux"]
tally2 = openmc.Tally(tally_id=2, name="heating tally")
cell_filter = openmc.CellFilter(cell)
tally2.filters = [cell_filter]
tally2.scores = ["heating"]
tallies = openmc.Tallies([tally1, tally2])
model = openmc.Model(
geometry=geometry, materials=[mat], settings=settings, tallies=tallies
)
sp_filename = model.run()
sp = openmc.StatePoint(sp_filename)
tally_found = sp.get_tally(filter_type=openmc.MeshFilter)
assert tally_found.id == 1
tally_found = sp.get_tally(filter_type=openmc.CellFilter)
assert tally_found.id == 2
tally_found = sp.get_tally(filters=[mesh_filter])
assert tally_found.id == 1
tally_found = sp.get_tally(filters=[cell_filter])
assert tally_found.id == 2
tally_found = sp.get_tally(scores=["heating"])
assert tally_found.id == 2
tally_found = sp.get_tally(name="heating tally")
assert tally_found.id == 2
tally_found = sp.get_tally(name=None)
assert tally_found.id == 1
tally_found = sp.get_tally(id=1)
assert tally_found.id == 1
tally_found = sp.get_tally(id=2)
assert tally_found.id == 2

10
tests/unit_tests/test_stats.py
View file

@ -16,6 +16,7 @@ def assert_sample_mean(samples, expected_mean):
assert np.abs(expected_mean - samples.mean()) < 4*std_dev
@pytest.mark.flaky(reruns=1)
def test_discrete():
x = [0.0, 1.0, 10.0]
p = [0.3, 0.2, 0.5]
@ -104,6 +105,7 @@ def test_clip_discrete():
d.clip(5)
@pytest.mark.flaky(reruns=1)
def test_uniform():
a, b = 10.0, 20.0
d = openmc.stats.Uniform(a, b)
@ -127,6 +129,7 @@ def test_uniform():
assert_sample_mean(samples, exp_mean)
@pytest.mark.flaky(reruns=1)
def test_powerlaw():
a, b, n = 10.0, 100.0, 2.0
d = openmc.stats.PowerLaw(a, b, n)
@ -148,6 +151,7 @@ def test_powerlaw():
assert_sample_mean(samples, exp_mean)
@pytest.mark.flaky(reruns=1)
def test_maxwell():
theta = 1.2895e6
d = openmc.stats.Maxwell(theta)
@ -171,6 +175,7 @@ def test_maxwell():
assert samples_2.mean() != samples.mean()
@pytest.mark.flaky(reruns=1)
def test_watt():
a, b = 0.965e6, 2.29e-6
d = openmc.stats.Watt(a, b)
@ -194,6 +199,7 @@ def test_watt():
assert_sample_mean(samples, exp_mean)
@pytest.mark.flaky(reruns=1)
def test_tabular():
# test linear-linear sampling
x = np.array([0.0, 5.0, 7.0, 10.0])
@ -270,6 +276,7 @@ def test_legendre():
d.to_xml_element('distribution')
@pytest.mark.flaky(reruns=1)
def test_mixture():
d1 = openmc.stats.Uniform(0, 5)
d2 = openmc.stats.Uniform(3, 7)
@ -425,6 +432,7 @@ def test_point():
assert d.xyz == pytest.approx(p)
@pytest.mark.flaky(reruns=1)
def test_normal():
mean = 10.0
std_dev = 2.0
@ -444,6 +452,7 @@ def test_normal():
assert_sample_mean(samples, mean)
@pytest.mark.flaky(reruns=1)
def test_muir():
mean = 10.0
mass = 5.0
@ -463,6 +472,7 @@ def test_muir():
assert_sample_mean(samples, mean)
@pytest.mark.flaky(reruns=1)
def test_combine_distributions():
# Combine two discrete (same data as in test_merge_discrete)
x1 = [0.0, 1.0, 10.0]

2
tools/ci/gha-script.sh
View file

@ -15,7 +15,7 @@ if [[ $EVENT == 'y' ]]; then
fi
# Run unit tests and then regression tests
pytest --cov=openmc -v $args \
pytest -v $args \
tests/test_matplotlib_import.py \
tests/unit_tests \
tests/regression_tests