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

Merge pull request #1186 from amandalund/density-effect-correction

Browse source 
TTB improvements
This commit is contained in:
Paul Romano 2019-03-07 12:56:26 -06:00 • committed by GitHub
commit 6250c526ec
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
17 changed files with 488 additions and 162 deletions
Download patch file Download diff file Expand all files Collapse all files

13
docs/source/io_formats/nuclear_data.rst
View file

@ -133,11 +133,17 @@ Incident Photon Data
**/<element>/bremsstrahlung/**
:Attributes: - **I** (*double*) -- Mean excitation energy in [eV]
:Datasets: - **electron_energy** (*double[]*) -- Incident electron energy in [eV]
- **photon_energy** (*double[]*) -- Outgoing photon energy as
fraction of incident electron energy
- **dcs** (*double[][]*) -- Bremsstrahlung differential cross section
at each incident energy in [mb/eV]
- **ionization_energy** (*double[]*) -- Ionization potential of each
subshell in [eV]
- **num_electrons** (*int[]*) -- Number of electrons per subshell,
with conduction electrons indicated by a negative value
**/<element>/coherent/**
@ -176,13 +182,6 @@ Incident Photon Data
:Datasets: - **xs** (*double[]*) -- Total photoionization cross section in [b]
**/<element>/stopping_powers/**
:Datasets: - **I** (*double*) -- Mean excitation energy in [eV]
- **energy** (*double[]*) -- Energies in [eV]
- **s_collision** (*double[]*) -- Collision stopping power in [eV-cm\ :sup:`2`\ /g]
- **s_radiative** (*double[]*) -- Radiative stopping power in [eV-cm\ :sup:`2`\ /g]
**/<element>/subshells/**
:Attributes: - **designators** (*char[][]*) -- Designator for each shell, e.g. 'M2'

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

@ -728,10 +728,18 @@ the cross section differential in energy loss. The total stopping power
power :math:`S_{\text{rad}}(T)`, which refers to energy loss due to
bremsstrahlung, and the collision stopping power :math:`S_{\text{col}}(T)`,
which refers to the energy loss due to inelastic collisions with bound
electrons in the material that result in ionization and excitation. To obtain
the radiative stopping power for positrons, the radiative stopping power for
electrons is multiplied by :eq:`positron-factor`. Currently, the collision
stopping power for electrons is also used for positrons.
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
@ -754,14 +762,97 @@ power is calculated using Bragg's additivity rule as
S_{\text{rad}}(T) = \sum_i w_i S_{\text{rad},i}(T),
where :math:`w_i` is the mass fraction of the :math:`i`-th element. The
collision stopping power, however, is a function of certain quantities such as
the mean excitation energy :math:`I` and the density effect correction
:math:`\delta_F` that depend on molecular properties. These quantities cannot
simply be summed over constituent elements in a compound, but should instead be
calculated for the material. Currently, we use Bragg's additivity rule to
calculate the collision stopping power as well, but this is not a good
approximation and should be fixed in the future.
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
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,
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:
@ -891,6 +982,55 @@ 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
-----------------
In coupled neutron-photon transport, a source neutron is tracked, and photons
produced from neutron reactions are transported after the neutron's history has
terminated. Since these secondary photons form the photon source for the
problem, it is important to correctly describe their energy and angular
distributions as the accuracy of the calculation relies on the accuracy of this
source. The photon production cross section for a particular reaction :math:`i`
and incident neutron energy :math:`E` is defined as
.. math::
:label: photon-production-xs
\sigma_{\gamma, i}(E) = y_i(E)\sigma_i(E),
where :math:`y_i(E)` is the photon yield corresponding to an incident neutron
reaction having cross section :math:`\sigma_i(E)`.
The yield of photons during neutron transport is determined as the sum of the
photon yields from each individual reaction. In OpenMC, production of photons
is treated in an average sense. That is, the total photon production cross
section is used at a collision site to determine how many photons to produce
rather than the photon production from the reaction that actually took place.
This is partly done for convenience but also because the use of variance
reduction techniques such as implicit capture make it difficult in practice to
directly sample photon production from individual reactions.
In OpenMC, secondary photons are created after a nuclide has been sampled in a
neutron collision. The expected number of photons produced is
.. math::
:label: expected-number-photons
n = w\frac{\sigma_{\gamma}(E)}{\sigma_T(E)},
where :math:`w` is the weight of the neutron, :math:`\sigma_{\gamma}` is the
photon production cross section for the sampled nuclide, and :math:`\sigma_T`
is the total cross section for the nuclide. :math:`\lfloor n \rfloor` photons
are created with an additional photon produced with probability :math:`n -
\lfloor n \rfloor`. Next, a reaction is sampled for each secondary photon. The
probability of sampling the :math:`i`-th reaction is given by
:math:`\sigma_{\gamma, i}(E)/\sum_j\sigma_{\gamma, j}(E)`, where
:math:`\sum_j\sigma_{\gamma, j} = \sigma_{\gamma}` is the total photon
production cross section. The secondary angle and energy distributions
associated with the reaction are used to sample the angle and energy of the
emitted photon.
.. _Koblinger: https://doi.org/10.13182/NSE75-A26663
.. _anomalous scattering: http://pd.chem.ucl.ac.uk/pdnn/diff1/anomscat.htm
@ -906,3 +1046,5 @@ angles.
.. _Kaltiaisenaho: https://aaltodoc.aalto.fi/bitstream/handle/123456789/21004/master_Kaltiaisenaho_Toni_2016.pdf
.. _Salvat: http://www.oecd-nea.org/globalsearch/download.php?doc=77434
.. _Sternheimer: https://doi.org/10.1103/PhysRevB.26.6067

10
docs/source/usersguide/cross_sections.rst
View file

@ -192,11 +192,11 @@ Photon Cross Sections
Photon interaction data is needed to run OpenMC with photon transport enabled.
Some of this data, namely bremsstrahlung cross sections from `Seltzer and
Berger`_, stopping powers from the `NIST ESTAR database`_, and Compton profiles
calculated by `Biggs et al.`_ and available in the Geant4 G4EMLOW data file, is
distributed with OpenMC. The rest is available from the NNDC_, which provides
ENDF data from the photo-atomic and atomic relaxation sublibraries of the
ENDF/B-VII.1 library.
Berger`_, mean excitation energy from the `NIST ESTAR database`_, and Compton
profiles calculated by `Biggs et al.`_ and available in the Geant4 G4EMLOW data
file, is distributed with OpenMC. The rest is available from the NNDC_, which
provides ENDF data from the photo-atomic and atomic relaxation sublibraries of
the ENDF/B-VII.1 library.
Most of the pregenerated HDF5 libraries available at https://openmc.mcs.anl.gov
already have photon interaction data included. If you are building a data

13
include/openmc/material.h
View file

@ -95,6 +95,9 @@ public:
std::unique_ptr<Bremsstrahlung> ttb_;
private:
//! Calculate the collision stopping power
void collision_stopping_power(double* s_col, bool positron);
//! Initialize bremsstrahlung data
void init_bremsstrahlung();
@ -109,6 +112,16 @@ private:
// Non-member functions
//==============================================================================
//! Calculate Sternheimer adjustment factor
double sternheimer_adjustment(const std::vector<double>& f, const
std::vector<double>& e_b_sq, double e_p_sq, double n_conduction, double
log_I, double tol, int max_iter);
//! Calculate density effect correction
double density_effect(const std::vector<double>& f, const std::vector<double>&
e_b_sq, double e_p_sq, double n_conduction, double rho, double E, double tol,
int max_iter);
//! Read material data from materials.xml
void read_materials_xml();

3
include/openmc/photon.h
View file

@ -87,7 +87,8 @@ public:
// Stopping power data
double I_; // mean excitation energy
xt::xtensor<double, 1> stopping_power_collision_;
xt::xtensor<int, 1> n_electrons_;
xt::xtensor<double, 1> ionization_energy_;
xt::xtensor<double, 1> stopping_power_radiative_;
// Bremsstrahlung scaled DCS

BIN
openmc/data/density_effect.h5 Normal file
View file

Binary file not shown.

91
openmc/data/photon.py
View file

@ -92,18 +92,15 @@ _REACTION_NAME = {
# is a 2D array with shape (n_shells, n_momentum_values) stored on the key Z
_COMPTON_PROFILES = {}
# Stopping powers are read from a pre-generated HDF5 file when they are first
# needed. The dictionary stores an array of energy values at which the other
# quantities are tabulated with the key 'energy' and for each element has the
# mean excitation energy and arrays containing the collision stopping powers
# and radiative stopping powers stored on the key 'Z'.
_STOPPING_POWERS = {}
# Scaled bremsstrahlung DCSs are read from a data file provided by Selzter and
# Berger when they are first needed. The dictionary stores an array of n
# incident electron kinetic energies with key 'electron_energies', an array of
# k reduced photon energies with key 'photon_energies', and the cross sections
# for each element are in a 2D array with shape (n, k) stored on the key 'Z'.
# It also stores data used for calculating the density effect correction and
# stopping power, namely, the mean excitation energy with the key 'I', number
# of electrons per subshell with the key 'num_electrons', and binding energies
# with the key 'ionization_energy'.
_BREMSSTRAHLUNG = {}
@ -352,12 +349,16 @@ class IncidentPhoton(EqualityMixin):
atomic_relaxation : openmc.data.AtomicRelaxation or None
Atomic relaxation data
bremsstrahlung : dict
Dictionary of bremsstrahlung DCS data with keys 'electron_energy'
(incident electron kinetic energy values in [eV]), 'photon_energy'
(ratio of the energy of the emitted photon to the incident electron
kinetic energy), and 'dcs' (cross section values in [b]). The cross
sections are in scaled form: :math:`(\beta^2/Z^2) E_k (d\sigma/dE_k)`,
where :math:`E_k` is the energy of the emitted photon.
Dictionary of bremsstrahlung data with keys 'I' (mean excitation energy
in [eV]), 'num_electrons' (number of electrons in each subshell),
'ionization_energy' (ionization potential of each subshell),
'electron_energy' (incident electron kinetic energy values in [eV]),
'photon_energy' (ratio of the energy of the emitted photon to the
incident electron kinetic energy), and 'dcs' (cross section values in
[b]). The cross sections are in scaled form: :math:`(\beta^2/Z^2) E_k
(d\sigma/dE_k)`, where :math:`E_k` is the energy of the emitted photon.
A negative number of electrons in a subshell indicates conduction
electrons.
compton_profiles : dict
Dictionary of Compton profile data with keys 'num_electrons' (number of
electrons in each subshell), 'binding_energy' (ionization potential of
@ -368,11 +369,6 @@ class IncidentPhoton(EqualityMixin):
reactions : collections.OrderedDict
Contains the cross sections for each photon reaction. The keys are MT
values and the values are instances of :class:`PhotonReaction`.
stopping_powers : dict
Dictionary of stopping power data with keys 'energy' (in [eV]), 'I' (mean
excitation energy), 's_collision' (collision stopping power in
[eV cm\ :sup:`2`/g]), and 's_radiative' (radiative stopping power in
[eV cm\ :sup:`2`/g])
"""
@ -381,7 +377,6 @@ class IncidentPhoton(EqualityMixin):
self._atomic_relaxation = None
self.reactions = OrderedDict()
self.compton_profiles = {}
self.stopping_powers = {}
self.bremsstrahlung = {}
def __contains__(self, mt):
@ -577,29 +572,19 @@ class IncidentPhoton(EqualityMixin):
data.compton_profiles['binding_energy'] = profile['binding_energy']
data.compton_profiles['J'] = [Tabulated1D(pz, J_k) for J_k in profile['J']]
# Load stopping power data if it has not yet been loaded
if not _STOPPING_POWERS:
filename = os.path.join(os.path.dirname(__file__), 'stopping_powers.h5')
with h5py.File(filename, 'r') as f:
# Units are in MeV; convert to eV
_STOPPING_POWERS['energy'] = f['energy'].value*EV_PER_MEV
for i in range(1, 99):
group = f['{:03}'.format(i)]
# Units are in MeV cm^2/g; convert to eV cm^2/g
_STOPPING_POWERS[i] = {
'I': group.attrs['I'],
's_collision': group['s_collision'].value*EV_PER_MEV,
's_radiative': group['s_radiative'].value*EV_PER_MEV
}
# Add stopping power data
if Z < 99:
data.stopping_powers['energy'] = _STOPPING_POWERS['energy']
data.stopping_powers.update(_STOPPING_POWERS[Z])
# Load bremsstrahlung data if it has not yet been loaded
if not _BREMSSTRAHLUNG:
# Add data used for density effect correction
filename = os.path.join(os.path.dirname(__file__), 'density_effect.h5')
with h5py.File(filename, 'r') as f:
for i in range(1, 101):
group = f['{:03}'.format(i)]
_BREMSSTRAHLUNG[i] = {
'I': group.attrs['I'],
'num_electrons': group['num_electrons'].value,
'ionization_energy': group['ionization_energy'].value
}
filename = os.path.join(os.path.dirname(__file__), 'BREMX.DAT')
brem = open(filename, 'r').read().split()
@ -640,12 +625,12 @@ class IncidentPhoton(EqualityMixin):
# Get scaled DCS values (millibarns) on new energy grid
dcs[:,j] = cs(log_energy)
_BREMSSTRAHLUNG[i] = {'dcs': dcs}
_BREMSSTRAHLUNG[i]['dcs'] = dcs
# Add bremsstrahlung DCS data
data.bremsstrahlung['electron_energy'] = _BREMSSTRAHLUNG['electron_energy']
data.bremsstrahlung['photon_energy'] = _BREMSSTRAHLUNG['photon_energy']
data.bremsstrahlung['dcs'] = _BREMSSTRAHLUNG[Z]['dcs']
data.bremsstrahlung.update(_BREMSSTRAHLUNG[Z])
return data
@ -767,26 +752,14 @@ class IncidentPhoton(EqualityMixin):
J = np.array([Jk.y for Jk in profile['J']])
compton_group.create_dataset('J', data=J)
# Write stopping powers
if self.stopping_powers:
s_group = group.create_group('stopping_powers')
for key, value in self.stopping_powers.items():
if key == 'I':
s_group.attrs[key] = value
else:
s_group.create_dataset(key, data=value)
# Write bremsstrahlung
if self.bremsstrahlung:
brem_group = group.create_group('bremsstrahlung')
brem = self.bremsstrahlung
brem_group.create_dataset('electron_energy',
data=brem['electron_energy'])
brem_group.create_dataset('photon_energy',
data=brem['photon_energy'])
brem_group.create_dataset('dcs', data=brem['dcs'])
for key, value in self.bremsstrahlung.items():
if key == 'I':
brem_group.attrs[key] = value
else:
brem_group.create_dataset(key, data=value)
class PhotonReaction(EqualityMixin):

BIN
openmc/data/stopping_powers.h5
View file

Binary file not shown.

243
src/material.cpp
View file

@ -451,6 +451,101 @@ void Material::init_thermal()
thermal_tables_ = tables;
}
void Material::collision_stopping_power(double* s_col, bool positron)
{
// Average electron number and average atomic weight
double electron_density = 0.0;
double mass_density = 0.0;
// Log of the mean excitation energy of the material
double log_I = 0.0;
// Effective number of conduction electrons in the material
double n_conduction = 0.0;
// Oscillator strength and square of the binding energy for each oscillator
// in material
std::vector<double> f;
std::vector<double> e_b_sq;
for (int i = 0; i < element_.size(); ++i) {
const auto& elm = data::elements[element_[i]];
double awr = data::nuclides[nuclide_[i]]->awr_;
// Get atomic density of nuclide given atom/weight percent
double atom_density = (atom_density_[0] > 0.0) ?
atom_density_[i] : -atom_density_[i] / awr;
electron_density += atom_density * elm.Z_;
mass_density += atom_density * awr * MASS_NEUTRON;
log_I += atom_density * elm.Z_ * std::log(elm.I_);
for (int j = 0; j < elm.n_electrons_.size(); ++j) {
if (elm.n_electrons_[j] < 0) {
n_conduction -= elm.n_electrons_[j] * atom_density;
continue;
}
e_b_sq.push_back(elm.ionization_energy_[j] * elm.ionization_energy_[j]);
f.push_back(elm.n_electrons_[j] * atom_density);
}
}
log_I /= electron_density;
n_conduction /= electron_density;
for (auto& f_i : f) f_i /= electron_density;
// Get density in g/cm^3 if it is given in atom/b-cm
double density = (density_ < 0.0) ? -density_ : mass_density / N_AVOGADRO;
// Calculate the square of the plasma energy
double e_p_sq = PLANCK_C * PLANCK_C * PLANCK_C * N_AVOGADRO *
electron_density * density / (2.0 * PI * PI * FINE_STRUCTURE *
MASS_ELECTRON_EV * mass_density);
// Get the Sternheimer adjustment factor
double rho = sternheimer_adjustment(f, e_b_sq, e_p_sq, n_conduction, log_I,
1.0e-6, 100);
// Classical electron radius in cm
constexpr double CM_PER_ANGSTROM {1.0e-8};
constexpr double r_e = CM_PER_ANGSTROM * PLANCK_C / (2.0 * PI *
FINE_STRUCTURE * MASS_ELECTRON_EV);
// Constant in expression for collision stopping power
constexpr double BARN_PER_CM_SQ {1.0e24};
double c = BARN_PER_CM_SQ * 2.0 * PI * r_e * r_e * MASS_ELECTRON_EV *
electron_density;
// Loop over incident charged particle energies
for (int i = 0; i < data::ttb_e_grid.size(); ++i) {
double E = data::ttb_e_grid(i);
// Get the density effect correction
double delta = density_effect(f, e_b_sq, e_p_sq, n_conduction, rho, E,
1.0e-6, 100);
// Square of the ratio of the speed of light to the velocity of the charged
// particle
double beta_sq = E * (E + 2.0 * MASS_ELECTRON_EV) / ((E + MASS_ELECTRON_EV)
* (E + MASS_ELECTRON_EV));
double tau = E / MASS_ELECTRON_EV;
double F;
if (positron) {
double t = tau + 2.0;
F = std::log(4.0) - (beta_sq / 12.0) * (23.0 + 14.0 / t + 10.0 / (t * t)
+ 4.0 / (t * t * t));
} else {
F = (1.0 - beta_sq) * (1.0 + tau * tau / 8.0 - (2.0 * tau + 1.0) *
std::log(2.0));
}
// Calculate the collision stopping power for this energy
s_col[i] = c / beta_sq * (2.0 * (std::log(E) - log_I) + std::log(1.0 + tau
/ 2.0) + F - delta);
}
}
void Material::init_bremsstrahlung()
{
// Create new object
@ -481,16 +576,11 @@ void Material::init_bremsstrahlung()
double Z_eq_sq = 0.0;
double sum_density = 0.0;
// Calculate the molecular DCS and the molecular total stopping power using
// Get the collision stopping power of the material
this->collision_stopping_power(stopping_power_collision.data(), positron);
// Calculate the molecular DCS and the molecular radiative stopping power using
// Bragg's additivity rule.
// TODO: The collision stopping power cannot be accurately calculated using
// Bragg's additivity rule since the mean excitation energies and the
// density effect corrections cannot simply be summed together. Bragg's
// additivity rule fails especially when a higher-density compound is
// composed of elements that are in lower-density form at normal temperature
// and pressure (at which the NIST stopping powers are given). It will be
// used to approximate the collision stopping powers for now, but should be
// fixed in the future.
for (int i = 0; i < n; ++i) {
// Get pointer to current element
const auto& elm = data::elements[element_[i]];
@ -499,7 +589,6 @@ void Material::init_bremsstrahlung()
// Get atomic density and mass density of nuclide given atom/weight percent
double atom_density = (atom_density_[0] > 0.0) ?
atom_density_[i] : -atom_density_[i] / awr;
double mass_density = atom_density * awr;
// Calculate the "equivalent" atomic number Zeq of the material
Z_eq_sq += atom_density * elm.Z_ * elm.Z_;
@ -508,13 +597,8 @@ void Material::init_bremsstrahlung()
// Accumulate material DCS
dcs += (atom_density * elm.Z_ * elm.Z_) * elm.dcs_;
// Accumulate material collision stopping power
stopping_power_collision += (mass_density * MASS_NEUTRON / N_AVOGADRO)
* elm.stopping_power_collision_;
// Accumulate material radiative stopping power
stopping_power_radiative += (mass_density * MASS_NEUTRON / N_AVOGADRO)
* elm.stopping_power_radiative_;
stopping_power_radiative += atom_density * elm.stopping_power_radiative_;
}
Z_eq_sq /= sum_density;
@ -569,12 +653,13 @@ void Material::init_bremsstrahlung()
// photon energy k
double x = x_l + (k - k_l)*(x_r - x_l)/(k_r - k_l);
// Ratio of the velocity of the charged particle to the speed of light
double beta = std::sqrt(e*(e + 2.0*MASS_ELECTRON_EV)) /
(e + MASS_ELECTRON_EV);
// Square of the ratio of the speed of light to the velocity of the
// charged particle
double beta_sq = e * (e + 2.0 * MASS_ELECTRON_EV) / ((e +
MASS_ELECTRON_EV) * (e + MASS_ELECTRON_EV));
// Compute the integrand of the PDF
f(j) = x / (beta*beta * stopping_power(j) * w);
f(j) = x / (beta_sq * stopping_power(j) * w);
}
// Number of points to integrate
@ -866,6 +951,124 @@ void Material::to_hdf5(hid_t group) const
// Non-method functions
//==============================================================================
double sternheimer_adjustment(const std::vector<double>& f, const
std::vector<double>& e_b_sq, double e_p_sq, double n_conduction, double
log_I, double tol, int max_iter)
{
// Get the total number of oscillators
int n = f.size();
// Calculate the Sternheimer adjustment factor using Newton's method
double rho = 2.0;
int iter;
for (iter = 0; iter < max_iter; ++iter) {
double rho_0 = rho;
// Function to find the root of and its derivative
double g = 0.0;
double gp = 0.0;
for (int i = 0; i < n; ++i) {
// Square of resonance energy of a bound-shell oscillator
double e_r_sq = e_b_sq[i] * rho * rho + 2.0 / 3.0 * f[i] * e_p_sq;
g += f[i] * std::log(e_r_sq);
gp += e_b_sq[i] * f[i] * rho / e_r_sq;
}
// Include conduction electrons
if (n_conduction > 0.0) {
g += n_conduction * std::log(n_conduction * e_p_sq);
}
// Set the next guess: rho_n+1 = rho_n - g(rho_n)/g'(rho_n)
rho -= (g - 2.0 * log_I) / (2.0 * gp);
// If the initial guess is too large, rho can be negative
if (rho < 0.0) rho = rho_0 / 2.0;
// Check for convergence
if (std::abs(rho - rho_0) / rho_0 < tol) break;
}
// Did not converge
if (iter >= max_iter) {
warning("Maximum Newton-Raphson iterations exceeded.");
rho = 1.0e-6;
}
return rho;
}
double density_effect(const std::vector<double>& f, const std::vector<double>&
e_b_sq, double e_p_sq, double n_conduction, double rho, double E, double tol,
int max_iter)
{
// Get the total number of oscillators
int n = f.size();
// Square of the ratio of the speed of light to the velocity of the charged
// particle
double beta_sq = E * (E + 2.0 * MASS_ELECTRON_EV) / ((E + MASS_ELECTRON_EV) *
(E + MASS_ELECTRON_EV));
// For nonmetals, delta = 0 for beta < beta_0, where beta_0 is obtained by
// setting the frequency w = 0.
double beta_0_sq = 0.0;
if (n_conduction == 0.0) {
for (int i = 0; i < n; ++i) {
beta_0_sq += f[i] * e_p_sq / (e_b_sq[i] * rho * rho);
}
beta_0_sq = 1.0 / (1.0 + beta_0_sq);
}
double delta = 0.0;
if (beta_sq < beta_0_sq) return delta;
// Compute the square of the frequency w^2 using Newton's method, with the
// initial guess of w^2 equal to beta^2 * gamma^2
double w_sq = E / MASS_ELECTRON_EV * (E / MASS_ELECTRON_EV + 2);
int iter;
for (iter = 0; iter < max_iter; ++iter) {
double w_sq_0 = w_sq;
// Function to find the root of and its derivative
double g = 0.0;
double gp = 0.0;
for (int i = 0; i < n; ++i) {
double c = e_b_sq[i] * rho * rho / e_p_sq + w_sq;
g += f[i] / c;
gp -= f[i] / (c * c);
}
// Include conduction electrons
g += n_conduction / w_sq;
gp -= n_conduction / (w_sq * w_sq);
// Set the next guess: w_n+1 = w_n - g(w_n)/g'(w_n)
w_sq -= (g + 1.0 - 1.0 / beta_sq) / gp;
// If the initial guess is too large, w can be negative
if (w_sq < 0.0) w_sq = w_sq_0 / 2.0;
// Check for convergence
if (std::abs(w_sq - w_sq_0) / w_sq_0 < tol) break;
}
// Did not converge
if (iter >= max_iter) {
warning("Maximum Newton-Raphson iterations exceeded: setting density "
"effect correction to zero.");
return delta;
}
// Solve for the density effect correction
for (int i = 0; i < n; ++i) {
double l_sq = e_b_sq[i] * rho * rho / e_p_sq + 2.0 / 3.0 * f[i];
delta += f[i] * std::log((l_sq + w_sq)/l_sq);
}
// Include conduction electrons
if (n_conduction > 0.0) {
delta += n_conduction * std::log((n_conduction + w_sq) / n_conduction);
}
return delta - w_sq * (1.0 - beta_sq);
}
void read_materials_xml()
{
write_message("Reading materials XML file...", 5);

53
src/photon.cpp
View file

@ -212,16 +212,12 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element)
if (data::ttb_k_grid.size() == 1) {
read_dataset(rgroup, "photon_energy", data::ttb_k_grid);
}
close_group(rgroup);
// Read stopping power data
if (Z_ < 99) {
rgroup = open_group(group, "stopping_powers");
read_dataset(rgroup, "s_collision", stopping_power_collision_);
read_dataset(rgroup, "s_radiative", stopping_power_radiative_);
read_attribute(rgroup, "I", I_);
close_group(rgroup);
}
// Get data used for density effect correction
read_dataset(rgroup, "num_electrons", n_electrons_);
read_dataset(rgroup, "ionization_energy", ionization_energy_);
read_attribute(rgroup, "I", I_);
close_group(rgroup);
// Truncate the bremsstrahlung data at the cutoff energy
int photon = static_cast<int>(Particle::Type::photon);
@ -235,26 +231,10 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element)
double f = (std::log(cutoff) - std::log(E(i_grid))) /
(std::log(E(i_grid+1)) - std::log(E(i_grid)));
// Interpolate collision stopping power at the cutoff energy and truncate
auto& s_col {stopping_power_collision_};
double y = std::exp(std::log(s_col(i_grid)) + f*(std::log(s_col(i_grid+1)) -
std::log(s_col(i_grid))));
xt::xtensor<double, 1> frst {y};
stopping_power_collision_ = xt::concatenate(xt::xtuple(
frst, xt::view(s_col, xt::range(i_grid+1, n_e))));
// Interpolate radiative stopping power at the cutoff energy and truncate
auto& s_rad {stopping_power_radiative_};
y = std::exp(std::log(s_rad(i_grid)) + f*(std::log(s_rad(i_grid+1)) -
std::log(s_rad(i_grid))));
frst(0) = y;
stopping_power_radiative_ = xt::concatenate(xt::xtuple(
frst, xt::view(s_rad, xt::range(i_grid+1, n_e))));
// Interpolate bremsstrahlung DCS at the cutoff energy and truncate
xt::xtensor<double, 2> dcs({n_e - i_grid, n_k});
for (int i = 0; i < n_k; ++i) {
y = std::exp(std::log(dcs_(i_grid,i)) +
double y = std::exp(std::log(dcs_(i_grid,i)) +
f*(std::log(dcs_(i_grid+1,i)) - std::log(dcs_(i_grid,i))));
auto col_i = xt::view(dcs, xt::all(), i);
col_i(0) = y;
@ -264,7 +244,7 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element)
}
dcs_ = dcs;
frst(0) = cutoff;
xt::xtensor<double, 1> frst {cutoff};
electron_energy = xt::concatenate(xt::xtuple(
frst, xt::view(electron_energy, xt::range(i_grid+1, n_e))));
}
@ -274,6 +254,25 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element)
if (data::ttb_e_grid.size() == 1) {
data::ttb_e_grid = electron_energy;
}
// Calculate the radiative stopping power
stopping_power_radiative_ = xt::empty<double>({data::ttb_e_grid.size()});
for (int i = 0; i < data::ttb_e_grid.size(); ++i) {
// Integrate over reduced photon energy
double c = 0.0;
for (int j = 0; j < data::ttb_k_grid.size() - 1; ++j) {
c += 0.5 * (dcs_(i, j+1) + dcs_(i, j)) * (data::ttb_k_grid(j+1) -
data::ttb_k_grid(j));
}
double e = data::ttb_e_grid(i);
// Square of the ratio of the speed of light to the velocity of the
// charged particle
double beta_sq = e * (e + 2.0 * MASS_ELECTRON_EV) / ((e +
MASS_ELECTRON_EV) * (e + MASS_ELECTRON_EV));
stopping_power_radiative_(i) = Z_ * Z_ / beta_sq * e * c;
}
}
// Take logarithm of energies and cross sections since they are log-log

13
src/physics.cpp
View file

@ -264,12 +264,15 @@ void sample_photon_reaction(Particle* p)
}
// Create Compton electron
double E_electron = (alpha - alpha_out)*MASS_ELECTRON_EV - e_b;
double mu_electron = (alpha - alpha_out*mu)
/ std::sqrt(alpha*alpha + alpha_out*alpha_out - 2.0*alpha*alpha_out*mu);
double phi = 2.0*PI*prn();
Direction u = rotate_angle(p->u(), mu_electron, &phi);
p->create_secondary(u, E_electron, Particle::Type::electron);
double E_electron = (alpha - alpha_out)*MASS_ELECTRON_EV - e_b;
int electron = static_cast<int>(Particle::Type::electron);
if (E_electron >= settings::energy_cutoff[electron]) {
double mu_electron = (alpha - alpha_out*mu)
/ std::sqrt(alpha*alpha + alpha_out*alpha_out - 2.0*alpha*alpha_out*mu);
Direction u = rotate_angle(p->u(), mu_electron, &phi);
p->create_secondary(u, E_electron, Particle::Type::electron);
}
// TODO: Compton subshell data does not match atomic relaxation data
// Allow electrons to fill orbital and produce auger electrons

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

@ -26,7 +26,7 @@
</space>
<angle reference_uvw="1.0 0.0 0.0" type="monodirectional" />
<energy type="discrete">
<parameters>14.0 1.0</parameters>
<parameters>14000000.0 1.0</parameters>
</energy>
</source>
<electron_treatment>ttb</electron_treatment>

4
tests/regression_tests/photon_production/results_true.dat
View file

@ -1,3 +1,3 @@
tally 1:
sum = 1.550000E-02
sum_sq = 2.402500E-04
sum = 7.938000E-01
sum_sq = 6.301184E-01

2
tests/regression_tests/photon_production/test.py
View file

@ -33,7 +33,7 @@ class SourceTestHarness(PyAPITestHarness):
source = openmc.Source()
source.space = openmc.stats.Point((0,0,0))
source.angle = openmc.stats.Monodirectional()
source.energy = openmc.stats.Discrete([14.0], [1.0])
source.energy = openmc.stats.Discrete([14.0e6], [1.0])
source.particle = 'neutron'
settings = openmc.Settings()

4
tests/regression_tests/photon_source/results_true.dat
View file

@ -1,3 +1,3 @@
tally 1:
sum = 2.254985E+02
sum_sq = 5.084955E+04
sum = 2.275713E+02
sum_sq = 5.178870E+04

31
tests/unit_tests/test_data_photon.py
View file

@ -71,10 +71,20 @@ def test_transitions(element):
assert sum(matrix['probability']) == pytest.approx(1.0)
@pytest.mark.parametrize('element', ['H', 'Al', 'Ag'], indirect=True)
def test_bremsstrahlung(element):
@pytest.mark.parametrize(
'element, I, i_shell, ionization_energy, num_electrons', [
('H', 19.2, 0, 13.6, 1),
('O', 95.0, 2, 13.62, 4),
('U', 890.0, 25, 6.033, -3)
],
indirect=['element']
)
def test_bremsstrahlung(element, I, i_shell, ionization_energy, num_electrons):
brems = element.bremsstrahlung
assert isinstance(brems, Mapping)
assert brems['I'] == I
assert brems['num_electrons'][i_shell] == num_electrons
assert brems['ionization_energy'][i_shell] == ionization_energy
assert np.all(np.diff(brems['electron_energy']) > 0.0)
assert np.all(np.diff(brems['photon_energy']) > 0.0)
assert brems['photon_energy'][0] == 0.0
@ -114,23 +124,6 @@ def test_reactions(element, reaction):
reactions[18]
@pytest.mark.parametrize(
'element, I', [
('H', 19.2),
('O', 95.0),
('U', 890.0)
],
indirect=['element']
)
def test_stopping_powers(element, I):
stopping_powers = element.stopping_powers
assert isinstance(stopping_powers, Mapping)
assert stopping_powers['I'] == I
assert np.all(np.diff(stopping_powers['energy']) > 0.0)
assert len(stopping_powers['s_collision']) == 200
assert len(stopping_powers['s_radiative']) == 200
@pytest.mark.parametrize('element', ['Pu'], indirect=True)
def test_export_to_hdf5(tmpdir, element):
filename = str(tmpdir.join('tmp.h5'))

2
tools/ci/download-xs.sh
View file

@ -3,7 +3,7 @@ set -ex
# Download HDF5 data
if [[ ! -e $HOME/nndc_hdf5/cross_sections.xml ]]; then
wget -q -O - https://anl.box.com/shared/static/snrluuy79o2fffpvpng9bsuis5kl811d.xz | tar -C $HOME -xJ
wget -q -O - https://anl.box.com/shared/static/9jmb8v2cx6kx03s6mbr2ai0mnvx5j79p.xz | tar -C $HOME -xJ
fi
# Download ENDF/B-VII.1 distribution