diff --git a/include/openmc/photon.h b/include/openmc/photon.h index ad621c59ec..7c20a716d6 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -89,8 +89,6 @@ public: double I_; // mean excitation energy xt::xtensor n_electrons_; xt::xtensor ionization_energy_; - // TODO: Calculate stopping powers instead of reading in - xt::xtensor stopping_power_collision_; 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 5e9f518bb0..01184a744c 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -92,22 +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, -# 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'. +# 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 = {} @@ -376,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]) """ @@ -389,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): @@ -585,27 +572,6 @@ 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 @@ -786,15 +752,6 @@ 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') 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 13b80ac977..277a568ea5 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -509,8 +509,7 @@ void Material::collision_stopping_power(double* s_col, bool positron) double r_e = PLANCK_C / (2.0e8 * PI * FINE_STRUCTURE * MASS_ELECTRON_EV); // Constant in expression for collision stopping power - double c = 2.0e24 * PI * r_e * r_e * MASS_ELECTRON_EV * N_AVOGADRO * - electron_density / mass_density; + double c = 2.0e24 * 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) { @@ -530,8 +529,8 @@ void Material::collision_stopping_power(double* s_col, bool positron) double F; if (positron) { double t = tau + 2.0; - F = 2.0 * std::log(2.0) - (beta_sq / 12.0) * (23.0 + 14.0 / t + 10.0 / - (t * t) + 4.0 / (t * t * t)); + 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)); @@ -573,16 +572,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(0), 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]]; @@ -591,7 +585,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_; @@ -600,13 +593,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; @@ -661,12 +649,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 diff --git a/src/photon.cpp b/src/photon.cpp index 49b050c61e..07753ce796 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -219,15 +219,6 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element) read_attribute(rgroup, "I", I_); 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); - } - // Truncate the bremsstrahlung data at the cutoff energy int photon = static_cast(ParticleType::photon); const auto& E {electron_energy}; @@ -240,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; @@ -269,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)))); } @@ -279,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