diff --git a/docs/source/io_formats/nuclear_data.rst b/docs/source/io_formats/nuclear_data.rst index 5595f2061a..81e4f4bcd7 100644 --- a/docs/source/io_formats/nuclear_data.rst +++ b/docs/source/io_formats/nuclear_data.rst @@ -133,11 +133,17 @@ Incident Photon Data **//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 **//coherent/** @@ -176,13 +182,6 @@ Incident Photon Data :Datasets: - **xs** (*double[]*) -- Total photoionization cross section in [b] -**//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] - **//subshells/** :Attributes: - **designators** (*char[][]*) -- Designator for each shell, e.g. 'M2' diff --git a/docs/source/methods/photon_physics.rst b/docs/source/methods/photon_physics.rst index 028a71de72..23809f6745 100644 --- a/docs/source/methods/photon_physics.rst +++ b/docs/source/methods/photon_physics.rst @@ -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 diff --git a/docs/source/usersguide/cross_sections.rst b/docs/source/usersguide/cross_sections.rst index dc0e9092ad..e6a79815f9 100644 --- a/docs/source/usersguide/cross_sections.rst +++ b/docs/source/usersguide/cross_sections.rst @@ -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 diff --git a/include/openmc/material.h b/include/openmc/material.h index 8a09dd9a24..932945ffce 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -95,6 +95,9 @@ public: std::unique_ptr 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& f, const + std::vector& 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& f, const std::vector& + 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(); diff --git a/include/openmc/photon.h b/include/openmc/photon.h index a2279f5027..7c20a716d6 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -87,7 +87,8 @@ public: // Stopping power data double I_; // mean excitation energy - xt::xtensor stopping_power_collision_; + xt::xtensor n_electrons_; + xt::xtensor ionization_energy_; xt::xtensor stopping_power_radiative_; // Bremsstrahlung scaled DCS diff --git a/openmc/data/density_effect.h5 b/openmc/data/density_effect.h5 new file mode 100644 index 0000000000..77d973517c Binary files /dev/null and b/openmc/data/density_effect.h5 differ diff --git a/openmc/data/photon.py b/openmc/data/photon.py index c73216b84d..01184a744c 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -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): diff --git a/openmc/data/stopping_powers.h5 b/openmc/data/stopping_powers.h5 deleted file mode 100644 index 4bc5b1388a..0000000000 Binary files a/openmc/data/stopping_powers.h5 and /dev/null differ diff --git a/src/material.cpp b/src/material.cpp index e6b2ca458e..da68a536dd 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -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 f; + std::vector 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& f, const + std::vector& 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& f, const std::vector& + 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); diff --git a/src/photon.cpp b/src/photon.cpp index 2dd0928ba4..d977d62ac3 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -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(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 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 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 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({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 diff --git a/src/physics.cpp b/src/physics.cpp index d898ed7881..85b6e2417e 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -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(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 diff --git a/tests/regression_tests/photon_production/inputs_true.dat b/tests/regression_tests/photon_production/inputs_true.dat index 32e3bdcf38..427078c889 100644 --- a/tests/regression_tests/photon_production/inputs_true.dat +++ b/tests/regression_tests/photon_production/inputs_true.dat @@ -26,7 +26,7 @@ - 14.0 1.0 + 14000000.0 1.0 ttb diff --git a/tests/regression_tests/photon_production/results_true.dat b/tests/regression_tests/photon_production/results_true.dat index 6f53e2d334..3d9cf68ede 100644 --- a/tests/regression_tests/photon_production/results_true.dat +++ b/tests/regression_tests/photon_production/results_true.dat @@ -1,3 +1,3 @@ tally 1: -sum = 1.550000E-02 -sum_sq = 2.402500E-04 +sum = 7.938000E-01 +sum_sq = 6.301184E-01 diff --git a/tests/regression_tests/photon_production/test.py b/tests/regression_tests/photon_production/test.py index 68992c867b..1abc184a2c 100644 --- a/tests/regression_tests/photon_production/test.py +++ b/tests/regression_tests/photon_production/test.py @@ -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() diff --git a/tests/regression_tests/photon_source/results_true.dat b/tests/regression_tests/photon_source/results_true.dat index 609f887bd1..1448ad6db8 100644 --- a/tests/regression_tests/photon_source/results_true.dat +++ b/tests/regression_tests/photon_source/results_true.dat @@ -1,3 +1,3 @@ tally 1: -sum = 2.254985E+02 -sum_sq = 5.084955E+04 +sum = 2.275713E+02 +sum_sq = 5.178870E+04 diff --git a/tests/unit_tests/test_data_photon.py b/tests/unit_tests/test_data_photon.py index 5036d08fd0..e2e7372864 100644 --- a/tests/unit_tests/test_data_photon.py +++ b/tests/unit_tests/test_data_photon.py @@ -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')) diff --git a/tools/ci/download-xs.sh b/tools/ci/download-xs.sh index 5fbe93c71d..b2068ceca4 100755 --- a/tools/ci/download-xs.sh +++ b/tools/ci/download-xs.sh @@ -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