diff --git a/include/openmc/bremsstrahlung.h b/include/openmc/bremsstrahlung.h new file mode 100644 index 000000000..ac74a2ab4 --- /dev/null +++ b/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 pdf; //!< Bremsstrahlung energy PDF + xt::xtensor cdf; //!< Bremsstrahlung energy CDF + xt::xtensor yield; //!< Photon yield +}; + +class Bremsstrahlung { +public: + // Data + BremsstrahlungData electron; + BremsstrahlungData positron; +}; + +//============================================================================== +// Global variables +//============================================================================== + +namespace data { + +extern xt::xtensor ttb_e_grid; //! energy T of incident electron in [eV] +extern xt::xtensor 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 diff --git a/include/openmc/material.h b/include/openmc/material.h index ad9b9c9f5..8df8fa308 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -1,11 +1,14 @@ #ifndef OPENMC_MATERIAL_H #define OPENMC_MATERIAL_H +#include // for unique_ptr #include #include #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 ttb_; + Material() {}; explicit Material(pugi::xml_node material_node); + +private: + void init_bremsstrahlung(); }; //============================================================================== diff --git a/include/openmc/photon.h b/include/openmc/photon.h index 81675585f..592ae6b97 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -119,8 +119,6 @@ struct ElementMicroXS { std::pair 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 compton_profile_pz; //! Compton profile momentum grid -extern xt::xtensor ttb_e_grid; //! energy T of incident electron in [eV] -extern xt::xtensor ttb_k_grid; //! reduced energy W/T of emitted photon //! Photon interaction data for each element extern std::vector elements; diff --git a/src/bremsstrahlung.cpp b/src/bremsstrahlung.cpp new file mode 100644 index 000000000..1ccf1f351 --- /dev/null +++ b/src/bremsstrahlung.cpp @@ -0,0 +1,109 @@ +#include "openmc/bremsstrahlung.h" + +namespace openmc { + +//============================================================================== +// Global variables +//============================================================================== + +namespace data { + +xt::xtensor ttb_e_grid; +xt::xtensor ttb_k_grid; +std::vector 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(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(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(ParticleType::photon); +// p.create_secondary(p.coord[0].uvw, w, photon, true); +// *E_lost += w; +// } +// } +// } + +} // namespace openmc diff --git a/src/material.cpp b/src/material.cpp index 7f8cfc50b..968094994 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -1,11 +1,19 @@ #include "openmc/material.h" +#include #include #include -#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(); + + // 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({n_e, n_e}); + ttb->cdf = xt::zeros({n_e, n_e}); + ttb->yield = xt::empty({n_e}); + + // Allocate temporary arrays + std::array shape {n_e}; + xt::xtensor stopping_power_collision({n_e}, 0.0); + xt::xtensor stopping_power_radiative({n_e}, 0.0); + xt::xtensor dcs({n_e, n_k}, 0.0); + xt::xtensor f({n_e}, 0.0); + xt::xtensor 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 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 //============================================================================== diff --git a/src/material_header.F90 b/src/material_header.F90 index 934173241..ee3f1f8bf 100644 --- a/src/material_header.F90 +++ b/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 diff --git a/src/photon.cpp b/src/photon.cpp index 7432a2dcd..b324f5ff3 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -768,93 +768,6 @@ std::pair 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(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(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(ParticleType::photon); -// p.create_secondary(p.coord[0].uvw, w, photon, true); -// *E_lost += w; -// } -// } -// } - //============================================================================== // Fortran compatibility //==============================================================================