Add function to calculate the density effect correction

include/openmc/material.h

@@ -95,6 +95,9 @@ public:
   std::unique_ptr<Bremsstrahlung> ttb_;
 
 private:
+  //! Calculate density effect correction
+  double density_effect_correction(double E);
+
   //! Initialize bremsstrahlung data
   void init_bremsstrahlung();
 

include/openmc/photon.h

@@ -87,6 +87,9 @@ public:
 
   // Stopping power data
   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_;
 

openmc/data/photon.py

@@ -104,6 +104,10 @@ _STOPPING_POWERS = {}
 # 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'.
 _BREMSSTRAHLUNG = {}
 
 
@@ -352,12 +356,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
@@ -600,6 +608,17 @@ class IncidentPhoton(EqualityMixin):
 
         # 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 +659,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
 
@@ -770,7 +789,6 @@ class IncidentPhoton(EqualityMixin):
         # 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
@@ -780,13 +798,11 @@ class IncidentPhoton(EqualityMixin):
         # 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):

src/material.cpp

@@ -451,6 +451,161 @@ void Material::init_thermal()
   thermal_tables_ = tables;
 }
 
+double Material::density_effect_correction(double E)
+{
+  // Maximum number of iterations and convergence criterion for the
+  // Newton-Raphson method
+  constexpr int max_iter = 100;
+  constexpr double tol = 1.0e-6;
+
+  // Effective number of conduction electrons in the material
+  double n_conduction = 0.0;
+
+  // Log of the mean excitation energy of the material
+  double log_I = 0.0;
+
+  // Average electron number and average atomic weight
+  double electron_density = 0.0;
+  double mass_density = 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 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
+    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;
+  }
+
+  // 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);
+  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
+  delta += n_conduction * std::log((n_conduction + w_sq) / n_conduction);
+
+  return delta - w_sq * (1.0 - beta_sq);
+}
+
 void Material::init_bremsstrahlung()
 {
   // Create new object

src/photon.cpp

@@ -212,6 +212,11 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element)
     if (data::ttb_k_grid.size() == 1) {
       read_dataset(rgroup, "photon_energy", data::ttb_k_grid);
     }
+
+    // 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);
 
     // Read stopping power data
@@ -219,7 +224,7 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element)
       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_);
+      //read_attribute(rgroup, "I", I_);
       close_group(rgroup);
     }