Calculate radiative and collision stopping power instead of using NIST ESTAR data

include/openmc/photon.h

@@ -89,8 +89,6 @@ public:
   double I_; // mean excitation energy
   xt::xtensor<int, 1> n_electrons_;
   xt::xtensor<double, 1> ionization_energy_;
-  // TODO: Calculate stopping powers instead of reading in
-  xt::xtensor<double, 1> stopping_power_collision_;
   xt::xtensor<double, 1> stopping_power_radiative_;
 
   // Bremsstrahlung scaled DCS

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')

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

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<int>(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<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;
@@ -269,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))));
     }
@@ -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<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