Start converting TTB data/physics

include/openmc/bremsstrahlung.h

@@ -0,0 +1,44 @@
+#ifndef OPENMC_BREMSSTRAHLUNG_H
+#define OPENMC_BREMSSTRAHLUNG_H
+
+#include "openmc/particle.h"
+
+#include "xtensor/xtensor.hpp"
+
+namespace openmc {
+
+//==============================================================================
+// Bremsstrahlung classes
+//==============================================================================
+
+class BremsstrahlungData {
+public:
+  // Data
+  xt::xtensor<double, 2> pdf; //!< Bremsstrahlung energy PDF
+  xt::xtensor<double, 2> cdf; //!< Bremsstrahlung energy CDF
+  xt::xtensor<double, 1> yield; //!< Photon yield
+};
+
+class Bremsstrahlung {
+public:
+  // Data
+  BremsstrahlungData electron;
+  BremsstrahlungData positron;
+};
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+namespace data {
+
+extern xt::xtensor<double, 1> ttb_e_grid; //! energy T of incident electron in [eV]
+extern xt::xtensor<double, 1> ttb_k_grid; //! reduced energy W/T of emitted photon
+
+} // namespace data
+
+void thick_target_bremsstrahlung(Particle& p, double* E_lost);
+
+} // namespace openmc
+
+#endif // OPENMC_BREMSSTRAHLUNG_H

include/openmc/material.h

@@ -1,11 +1,14 @@
 #ifndef OPENMC_MATERIAL_H
 #define OPENMC_MATERIAL_H
 
+#include <memory> // for unique_ptr
 #include <unordered_map>
 #include <vector>
 
 #include "pugixml.hpp"
 
+#include "openmc/bremsstrahlung.h"
+#include "openmc/particle.h"
 
 namespace openmc {
 
@@ -38,9 +41,14 @@ public:
   //! A negative value indicates no default temperature was specified.
   double temperature_ {-1};
 
+  std::unique_ptr<Bremsstrahlung> ttb_;
+
   Material() {};
 
   explicit Material(pugi::xml_node material_node);
+
+private:
+  void init_bremsstrahlung();
 };
 
 //==============================================================================

include/openmc/photon.h

@@ -119,8 +119,6 @@ struct ElementMicroXS {
 
 std::pair<double, double> klein_nishina(double alpha);
 
-void thick_target_bremsstrahlung(Particle& p, double* E_lost);
-
 //==============================================================================
 // Global variables
 //==============================================================================
@@ -128,8 +126,6 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost);
 namespace data {
 
 extern xt::xtensor<double, 1> compton_profile_pz; //! Compton profile momentum grid
-extern xt::xtensor<double, 1> ttb_e_grid; //! energy T of incident electron in [eV]
-extern xt::xtensor<double, 1> ttb_k_grid; //! reduced energy W/T of emitted photon
 
 //! Photon interaction data for each element
 extern std::vector<PhotonInteraction> elements;

src/bremsstrahlung.cpp

@@ -0,0 +1,109 @@
+#include "openmc/bremsstrahlung.h"
+
+namespace openmc {
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+namespace data {
+
+xt::xtensor<double, 1> ttb_e_grid;
+xt::xtensor<double, 1> ttb_k_grid;
+std::vector<Bremsstrahlung> ttb;
+
+} // namespace data
+
+//==============================================================================
+// Non-member functions
+//==============================================================================
+
+
+// void thick_target_bremsstrahlung(Particle& p, double* E_lost)
+// {
+//   if (p.material == MATERIAL_VOID) return;
+
+//   // TODO: off-by-one
+//   int photon = static_cast<int>(ParticleType::photon) - 1;
+//   if (p.E < settings::energy_cutoff[photon]) return;
+
+//   // // Get bremsstrahlung data for this material and particle type
+//   if (p.type == static_cast<int>(ParticleType::positron)) {
+//     mat => ttb(p.material) % positron;
+//   } else {
+//     mat => ttb(p.material) % electron;
+//   }
+
+//   double e = std::log(p.E);
+//   auto n_e = data::ttb_e_grid
+
+//   // Find the lower bounding index of the incident electron energy
+//   int j = lower_bound_index(data::ttb_e_grid.cbegin(), data::ttb_e_grid.cend(), e);
+//   if (j == n_e - 1) --j;
+
+//   // Get the interpolation bounds
+//   double e_l = data::ttb_e_grid(j);
+//   double e_r = data::ttb_e_grid(j+1);
+//   double y_l = mat % yield(j);
+//   double y_r = mat % yield(j+1);
+
+//   // Calculate the interpolation weight w_j+1 of the bremsstrahlung energy PDF
+//   // interpolated in log energy, which can be interpreted as the probability
+//   // of index j+1
+//   double f = (e - e_l)/(e_r - e_l);
+
+//   // Get the photon number yield for the given energy using linear
+//   // interpolation on a log-log scale
+//   double y = std::exp(y_l + (y_r - y_l)*f);
+
+//   // Sample number of secondary bremsstrahlung photons
+//   int n = y + prn();
+
+//   *E_lost = 0.0;
+//   if (n == 0) return;
+
+//   // Sample index of the tabulated PDF in the energy grid, j or j+1
+//   double c_max;
+//   if (prn() <= f || j == 0) {
+//     int i_e = j + 1;
+
+//     // Interpolate the maximum value of the CDF at the incoming particle
+//     // energy on a log-log scale
+//     double p_l = mat % pdf(i_e-1, i_e);
+//     double p_r = mat % pdf(i_e, i_e);
+//     double c_l = mat % cdf(i_e-1, i_e);
+//     double a = std::log(p_r/p_l)/(e_r - e_l) + 1.0;
+//     c_max = c_l + std::exp(e_l)*p_l/a*(std::exp(a*(e - e_l)) - 1.0);
+//   } else {
+//     int i_e = j;
+
+//     // Maximum value of the CDF
+//     c_max = mat % cdf(i_e, i_e)
+//   }
+
+//   // Sample the energies of the emitted photons
+//   for (int i = 0; i < n; ++i) {
+//     // Generate a random number r and determine the index i for which
+//     // cdf(i) <= r*cdf,max <= cdf(i+1)
+//     double c = prn()*c_max;
+//     int i_w = lower_bound_index(mat % cdf(:i_e,i_e), c)
+
+//     // Sample the photon energy
+//     double w_l = data::ttb_e_grid(i_w);
+//     double w_r = data::ttb_e_grid(i_w+1);
+//     double p_l = mat % pdf(i_w, i_e);
+//     double p_r = mat % pdf(i_w+1, i_e);
+//     double c_l = mat % cdf(i_w, i_e);
+//     double a = std::log(p_r/p_l)/(w_r - w_l) + 1.0;
+//     double w = std::exp(w_l)*std::pow(a*(c - c_l)/(std::exp(w_l)*p_l) + 1.0, 1.0/a);
+
+//     if (w > settings::energy_cutoff[photon]) {
+//       // Create secondary photon
+//       int photon = static_cast<int>(ParticleType::photon);
+//       p.create_secondary(p.coord[0].uvw, w, photon, true);
+//       *E_lost += w;
+//     }
+//   }
+// }
+
+} // namespace openmc

src/material.cpp

@@ -1,11 +1,19 @@
 #include "openmc/material.h"
 
+#include <cmath>
 #include <string>
 #include <sstream>
 
-#include "openmc/error.h"
-#include "openmc/xml_interface.h"
+#include "xtensor/xbuilder.hpp"
+#include "xtensor/xoperation.hpp"
+#include "xtensor/xview.hpp"
 
+#include "openmc/error.h"
+#include "openmc/math_functions.h"
+#include "openmc/nuclide.h"
+#include "openmc/photon.h"
+#include "openmc/search.h"
+#include "openmc/xml_interface.h"
 
 namespace openmc {
 
@@ -41,6 +49,197 @@ Material::Material(pugi::xml_node node)
   }
 }
 
+int* material_element(int index);
+
+void Material::init_bremsstrahlung()
+{
+  // Create new object
+  ttb_ = std::make_unique<Bremsstrahlung>();
+
+  // Get the size of the energy grids
+  auto n_k = data::ttb_k_grid.size();
+  auto n_e = data::ttb_e_grid.size();
+
+  // Get pointers to nuclides, elements, densities
+  int32_t index;
+  openmc_get_material_index(id_, &index);
+  int* nuclide_;
+  double* atom_density_;
+  int n;
+  openmc_material_get_densities(index, &nuclide_, &atom_density_, &n);
+  int* element_ = material_element(index);
+
+  for (int particle = 0; particle < 2; ++particle) {
+    // Loop over logic twice, once for electron, once for positron
+    BremsstrahlungData* ttb = (particle == 0) ? &ttb_->electron : &ttb_->positron;
+    bool positron = (particle == 1);
+
+    // Allocate arrays for TTB data
+    ttb->pdf = xt::zeros<double>({n_e, n_e});
+    ttb->cdf = xt::zeros<double>({n_e, n_e});
+    ttb->yield = xt::empty<double>({n_e});
+
+    // Allocate temporary arrays
+    std::array<std::size_t, 1> shape {n_e};
+    xt::xtensor<double, 1> stopping_power_collision({n_e}, 0.0);
+    xt::xtensor<double, 1> stopping_power_radiative({n_e}, 0.0);
+    xt::xtensor<double, 2> dcs({n_e, n_k}, 0.0);
+    xt::xtensor<double, 1> f({n_e}, 0.0);
+    xt::xtensor<double, 1> z({n_e}, 0.0);
+
+    double Z_eq_sq = 0.0;
+    double sum_density = 0.0;
+
+    // Calculate the molecular DCS and the molecular total 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
+      // TODO: off-by-one
+      const auto& elm = data::elements[element_[i] - 1];
+      // TODO: off-by-one
+      double awr = data::nuclides[nuclide_[i] - 1]->awr_;
+
+      // 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_;
+      sum_density += atom_density;
+
+      // 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_;
+    }
+    Z_eq_sq /= sum_density;
+
+    // Calculate the positron DCS and radiative stopping power. These are
+    // obtained by multiplying the electron DCS and radiative stopping powers by
+    // a factor r, which is a numerical approximation of the ratio of the
+    // radiative stopping powers for positrons and electrons. Source: F. Salvat,
+    // J. M. Fernández-Varea, and J. Sempau, "PENELOPE-2011: A Code System for
+    // Monte Carlo Simulation of Electron and Photon Transport," OECD-NEA,
+    // Issy-les-Moulineaux, France (2011).
+    if (positron) {
+      for (int i = 0; i < n_e; ++i) {
+        double t = std::log(1.0 + 1.0e6*data::ttb_e_grid(i)/(Z_eq_sq*MASS_ELECTRON_EV));
+        double r = 1.0 - std::exp(-1.2359e-1*t + 6.1274e-2*std::pow(t, 2)
+          - 3.1516e-2*std::pow(t, 3) + 7.7446e-3*std::pow(t, 4)
+          - 1.0595e-3*std::pow(t, 5) + 7.0568e-5*std::pow(t, 6)
+          - 1.808e-6*std::pow(t, 7));
+        stopping_power_radiative(i) *= r;
+        auto dcs_i = xt::view(dcs, i, xt::all());
+        dcs_i *= r;
+      }
+    }
+
+    // Total material stopping power
+    xt::xtensor<double, 1> stopping_power = stopping_power_collision +
+      stopping_power_radiative;
+
+    // Loop over photon energies
+    for (int i = 0; i < n_e - 1; ++i) {
+      double w = data::ttb_e_grid(i);
+
+      // Loop over incident particle energies
+      for (int j = 0; j < n_e; ++j) {
+        double e = data::ttb_e_grid(j);
+
+        // Reduced photon energy
+        double k = w / e;
+
+        // Find the lower bounding index of the reduced photon energy
+        int i_k = lower_bound_index(data::ttb_k_grid.cbegin(), data::ttb_k_grid.cend(), k);
+
+        // Get the interpolation bounds
+        double k_l = data::ttb_k_grid(i_k);
+        double k_r = data::ttb_k_grid(i_k + 1);
+        double x_l = dcs(j, i_k);
+        double x_r = dcs(j, i_k+1);
+
+        // Find the value of the DCS using linear interpolation in reduced
+        // 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);
+
+        // Compute the integrand of the PDF
+        f(j) = x / (beta*beta * stopping_power(j) * w);
+      }
+
+      // Number of points to integrate
+      int n = n_e - i;
+
+      // Integrate the PDF using cubic spline integration over the incident
+      // particle energy
+      if (n > 2) {
+        spline_c(n, &data::ttb_e_grid(i), &f(i), &z(i));
+
+        double c = 0.0;
+        for (int j = 0; j < n_e - 1; ++j) {
+          c += spline_integrate_c(n, &data::ttb_e_grid(i), &f(i), &z(i),
+            data::ttb_e_grid(j), data::ttb_e_grid(j+1));
+
+          ttb->pdf(j+1,i) = c;
+        }
+
+      // Integrate the last two points using trapezoidal rule in log-log space
+      } else {
+        double e_l = std::log(data::ttb_e_grid(i));
+        double e_r = std::log(data::ttb_e_grid(i+1));
+        double x_l = std::log(f(i));
+        double x_r = std::log(f(i+1));
+
+        ttb->pdf(i,i+1) = 0.5*(e_r - e_l)*(std::exp(e_l + x_l) + std::exp(e_r + x_r));
+      }
+    }
+
+    // Loop over incident particle energies
+    for (int j = 1; j < n_e; ++j) {
+      // Set last element of PDF to small non-zero value to enable log-log
+      // interpolation
+      ttb->pdf(j,j) = std::exp(-500.0);
+
+      // Loop over photon energies
+      double c = 0.0;
+      for (int i = 0; i < j - 1; ++i) {
+        // Integrate the CDF from the PDF using the trapezoidal rule in log-log
+        // space
+        double w_l = std::log(data::ttb_e_grid(i));
+        double w_r = std::log(data::ttb_e_grid(i+1));
+        double x_l = std::log(ttb->pdf(j,i));
+        double x_r = std::log(ttb->pdf(j,i+1));
+
+        c += 0.5*(w_r - w_l)*(std::exp(w_l + x_l) + std::exp(w_r + x_r));
+        ttb->cdf(j,i+1) = c;
+      }
+
+      // Set photon number yield
+      ttb->yield(j) = c;
+    }
+
+    // Use logarithm of number yield since it is log-log interpolated
+    ttb->yield = xt::where(ttb->yield > 0.0, ttb->yield, -500.0);
+  }
+}
+
 //==============================================================================
 // Non-method functions
 //==============================================================================

src/material_header.F90

@@ -1065,4 +1065,10 @@ contains
     end associate
   end function
 
+  function material_element(i_material) result(ptr) bind(C)
+    integer(C_INT), value :: i_material
+    type(C_PTR) :: ptr
+    ptr = C_LOC(materials(i_material) % element(1))
+  end function
+
 end module material_header

src/photon.cpp

@@ -768,93 +768,6 @@ std::pair<double, double> klein_nishina(double alpha)
   return {alpha_out, mu};
 }
 
-// void thick_target_bremsstrahlung(Particle& p, double* E_lost)
-// {
-//   if (p.material == MATERIAL_VOID) return;
-
-//   // TODO: off-by-one
-//   int photon = static_cast<int>(ParticleType::photon) - 1;
-//   if (p.E < settings::energy_cutoff[photon]) return;
-
-//   // // Get bremsstrahlung data for this material and particle type
-//   if (p.type == static_cast<int>(ParticleType::positron)) {
-//     mat => ttb(p.material) % positron;
-//   } else {
-//     mat => ttb(p.material) % electron;
-//   }
-
-//   double e = std::log(p.E);
-//   auto n_e = data::ttb_e_grid
-
-//   // Find the lower bounding index of the incident electron energy
-//   int j = lower_bound_index(data::ttb_e_grid.cbegin(), data::ttb_e_grid.cend(), e);
-//   if (j == n_e - 1) --j;
-
-//   // Get the interpolation bounds
-//   double e_l = data::ttb_e_grid(j);
-//   double e_r = data::ttb_e_grid(j+1);
-//   double y_l = mat % yield(j);
-//   double y_r = mat % yield(j+1);
-
-//   // Calculate the interpolation weight w_j+1 of the bremsstrahlung energy PDF
-//   // interpolated in log energy, which can be interpreted as the probability
-//   // of index j+1
-//   double f = (e - e_l)/(e_r - e_l);
-
-//   // Get the photon number yield for the given energy using linear
-//   // interpolation on a log-log scale
-//   double y = std::exp(y_l + (y_r - y_l)*f);
-
-//   // Sample number of secondary bremsstrahlung photons
-//   int n = y + prn();
-
-//   *E_lost = 0.0;
-//   if (n == 0) return;
-
-//   // Sample index of the tabulated PDF in the energy grid, j or j+1
-//   double c_max;
-//   if (prn() <= f || j == 0) {
-//     int i_e = j + 1;
-
-//     // Interpolate the maximum value of the CDF at the incoming particle
-//     // energy on a log-log scale
-//     double p_l = mat % pdf(i_e-1, i_e);
-//     double p_r = mat % pdf(i_e, i_e);
-//     double c_l = mat % cdf(i_e-1, i_e);
-//     double a = std::log(p_r/p_l)/(e_r - e_l) + 1.0;
-//     c_max = c_l + std::exp(e_l)*p_l/a*(std::exp(a*(e - e_l)) - 1.0);
-//   } else {
-//     int i_e = j;
-
-//     // Maximum value of the CDF
-//     c_max = mat % cdf(i_e, i_e)
-//   }
-
-//   // Sample the energies of the emitted photons
-//   for (int i = 0; i < n; ++i) {
-//     // Generate a random number r and determine the index i for which
-//     // cdf(i) <= r*cdf,max <= cdf(i+1)
-//     double c = prn()*c_max;
-//     int i_w = lower_bound_index(mat % cdf(:i_e,i_e), c)
-
-//     // Sample the photon energy
-//     double w_l = data::ttb_e_grid(i_w);
-//     double w_r = data::ttb_e_grid(i_w+1);
-//     double p_l = mat % pdf(i_w, i_e);
-//     double p_r = mat % pdf(i_w+1, i_e);
-//     double c_l = mat % cdf(i_w, i_e);
-//     double a = std::log(p_r/p_l)/(w_r - w_l) + 1.0;
-//     double w = std::exp(w_l)*std::pow(a*(c - c_l)/(std::exp(w_l)*p_l) + 1.0, 1.0/a);
-
-//     if (w > settings::energy_cutoff[photon]) {
-//       // Create secondary photon
-//       int photon = static_cast<int>(ParticleType::photon);
-//       p.create_secondary(p.coord[0].uvw, w, photon, true);
-//       *E_lost += w;
-//     }
-//   }
-// }
-
 //==============================================================================
 // Fortran compatibility
 //==============================================================================