diff --git a/include/openmc/endf.h b/include/openmc/endf.h index ae480f8cf..7ec078e93 100644 --- a/include/openmc/endf.h +++ b/include/openmc/endf.h @@ -78,6 +78,10 @@ public: //! \param[in] x independent variable //! \return Function evaluated at x double operator()(double x) const; + + // Accessors + const std::vector& x() const { return x_; } + const std::vector& y() const { return y_; } private: std::size_t n_regions_ {0}; //!< number of interpolation regions std::vector nbt_; //!< values separating interpolation regions diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h index bc5ba9308..f80ed6d3b 100644 --- a/include/openmc/math_functions.h +++ b/include/openmc/math_functions.h @@ -129,7 +129,7 @@ extern "C" void calc_zn_rad(int n, double rho, double zn_rad[]); //! is passed //============================================================================== -extern "C" void rotate_angle_c(double uvw[3], double mu, double* phi); +extern "C" void rotate_angle_c(double uvw[3], double mu, const double* phi); Direction rotate_angle(Direction u, double mu, double* phi); diff --git a/include/openmc/photon.h b/include/openmc/photon.h index bc3f33e14..fce6a3626 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -2,12 +2,14 @@ #define OPENMC_PHOTON_H #include "openmc/endf.h" +#include "openmc/particle.h" #include #include "xtensor/xtensor.hpp" #include #include +#include // for pair #include namespace openmc { @@ -39,6 +41,17 @@ public: // Constructors PhotonInteraction(hid_t group); + // Methods + void compton_scatter(double alpha, bool doppler, double* alpha_out, + double* mu, int* i_shell) const; + + double rayleigh_scatter(double alpha) const; + + void pair_production(double alpha, double* E_electron, double* E_positron, + double* mu_electron, double* mu_positron) const; + + void atomic_relaxation(const ElectronSubshell& shell, Particle& p) const; + // Data members std::string name_; //! Name of element, e.g. "Zr" int Z_; //! Atomic number @@ -76,6 +89,9 @@ public: // Bremsstrahlung scaled DCS xt::xtensor dcs_; + +private: + void compton_doppler(double alpha, double mu, double* E_out, int* i_shell) const; }; //============================================================================== @@ -94,6 +110,14 @@ struct ElementMicroXS { double pair_production; //!< microscopic pair production xs }; +//============================================================================== +// Non-member functions +//============================================================================== + +std::pair klein_nishina(double alpha); + +void thick_target_bremsstrahlung(Particle& p, double* E_lost); + //============================================================================== // Global variables //============================================================================== diff --git a/include/openmc/physics.h b/include/openmc/physics.h index a51ed91b8..8a3d7c4ef 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -22,7 +22,7 @@ void sample_neutron_reaction(Particle* p); //! Samples an element based on the macroscopic cross sections for each nuclide //! within a material and then samples a reaction for that element and calls the //! appropriate routine to process the physics. -extern "C" void sample_photon_reaction(Particle* p); +void sample_photon_reaction(Particle* p); //! Terminates the particle and either deposits all energy locally //! (electron_treatment = ELECTRON_LED) or creates secondary bremsstrahlung @@ -44,7 +44,7 @@ void sample_nuclide(const Particle* p, int mt, int* i_nuclide, int* i_nuc_mat); void create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, Bank* bank_array, int64_t* bank_size, int64_t bank_capacity); -// void sample_element(Particle* p); +extern "C" int sample_element(Particle* p); Reaction* sample_fission(int i_nuclide, double E); diff --git a/src/math_functions.cpp b/src/math_functions.cpp index dcf8d38b1..041b6ff57 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -617,7 +617,7 @@ void calc_zn_rad(int n, double rho, double zn_rad[]) { } -void rotate_angle_c(double uvw[3], double mu, double* phi) { +void rotate_angle_c(double uvw[3], double mu, const double* phi) { // Copy original directional cosines double u0 = uvw[0]; // original cosine in x direction double v0 = uvw[1]; // original cosine in y direction diff --git a/src/photon.cpp b/src/photon.cpp index edb0de4e4..56ff2042d 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -3,6 +3,7 @@ #include "openmc/constants.h" #include "openmc/hdf5_interface.h" #include "openmc/particle.h" +#include "openmc/random_lcg.h" #include "openmc/search.h" #include "openmc/settings.h" @@ -10,6 +11,10 @@ #include "xtensor/xoperation.hpp" #include "xtensor/xview.hpp" +#include +#include +#include // for tie + namespace openmc { //============================================================================== @@ -107,7 +112,7 @@ PhotonInteraction::PhotonInteraction(hid_t group) shells_.emplace_back(); // Create mapping from designator to index - int j = 0; + int j = 1; for (const auto& subshell : SUBSHELLS) { if (designator == subshell) { shell_map_[j] = i; @@ -126,8 +131,7 @@ PhotonInteraction::PhotonInteraction(hid_t group) // Read subshell cross section dset = open_dataset(tgroup, "xs"); read_attribute(dset, "threshold_idx", j); - // TODO: off-by-one - shell.threshold = j - 1; + shell.threshold = j; close_dataset(dset); read_dataset(tgroup, "xs", shell.cross_section); @@ -151,6 +155,8 @@ PhotonInteraction::PhotonInteraction(hid_t group) shell.transition_probability = xt::view(matrix, xt::all(), 3); shell.transition_probability /= xt::sum(shell.transition_probability)(); } + } else { + shell.n_transitions = 0; } close_group(tgroup); } @@ -277,6 +283,517 @@ PhotonInteraction::PhotonInteraction(hid_t group) xt::log(pair_production_total_), -500.0); } +void PhotonInteraction::compton_scatter(double alpha, bool doppler, + double* alpha_out, double* mu, int* i_shell) const +{ + double form_factor_xmax = 0.0; + while (true) { + // Sample Klein-Nishina distribution for trial energy and angle + std::tie(*alpha_out, *mu) = klein_nishina(alpha); + + // Note that the parameter used here does not correspond exactly to the + // momentum transfer q in ENDF-102 Eq. (27.2). Rather, this is the + // parameter as defined by Hubbell, where the actual data comes from + double x = MASS_ELECTRON_EV/PLANCK_C*alpha*std::sqrt(0.5*(1.0 - *mu)); + + // Calculate S(x, Z) and S(x_max, Z) + double form_factor_x = incoherent_form_factor_(x); + if (form_factor_xmax == 0.0) { + form_factor_xmax = incoherent_form_factor_(MASS_ELECTRON_EV/PLANCK_C*alpha); + } + + // Perform rejection on form factor + if (prn() < form_factor_x / form_factor_xmax) { + if (doppler) { + double E_out; + this->compton_doppler(alpha, *mu, &E_out, i_shell); + *alpha_out = E_out/MASS_ELECTRON_EV; + } else { + *i_shell = -1; + } + break; + } + } +} + +void PhotonInteraction::compton_doppler(double alpha, double mu, + double* E_out, int* i_shell) const +{ + auto n = data::compton_profile_pz.size(); + + int shell; // index for shell + while (true) { + // Sample electron shell + double rn = prn(); + double c = 0.0; + for (shell = 0; shell < electron_pdf_.size(); ++shell) { + c += electron_pdf_(shell); + if (rn < c) break; + } + + // Determine binding energy of shell + double E_b = binding_energy_(shell); + + // Determine p_z,max + double E = alpha*MASS_ELECTRON_EV; + if (E < E_b) { + *E_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON_EV; + break; + } + + double pz_max = -FINE_STRUCTURE*(E_b - (E - E_b)*alpha*(1.0 - mu)) / + std::sqrt(2.0*E*(E - E_b)*(1.0 - mu) + E_b*E_b); + if (pz_max < 0.0) { + *E_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON_EV; + break; + } + + // Determine profile cdf value corresponding to p_z,max + double c_max; + if (pz_max > data::compton_profile_pz(n - 1)) { + c_max = profile_cdf_(shell, n - 1); + } else { + int i = lower_bound_index(data::compton_profile_pz.cbegin(), + data::compton_profile_pz.cend(), pz_max); + double pz_l = data::compton_profile_pz(i); + double pz_r = data::compton_profile_pz(i + 1); + double p_l = profile_pdf_(shell, i); + double p_r = profile_pdf_(shell, i + 1); + double c_l = profile_cdf_(shell, i); + if (pz_l == pz_r) { + c_max = c_l; + } else if (p_l == p_r) { + c_max = c_l + (pz_max - pz_l)*p_l; + } else { + double m = (p_l - p_r)/(pz_l - pz_r); + c_max = c_l + (std::pow((m*(pz_max - pz_l) + p_l), 2) - p_l*p_l)/(2.0*m); + } + } + + // Sample value on bounded cdf + c = prn()*c_max; + + // Determine pz corresponding to sampled cdf value + auto cdf_shell = xt::view(profile_cdf_, shell, xt::all()); + int i = lower_bound_index(cdf_shell.cbegin(), cdf_shell.cend(), c); + double pz_l = data::compton_profile_pz(i); + double pz_r = data::compton_profile_pz(i + 1); + double p_l = profile_pdf_(shell, i); + double p_r = profile_pdf_(shell, i + 1); + double c_l = profile_cdf_(shell, i); + double pz; + if (pz_l == pz_r) { + pz = pz_l; + } else if (p_l == p_r) { + pz = pz_l + (c - c_l)/p_l; + } else { + double m = (p_l - p_r)/(pz_l - pz_r); + pz = pz_l + (std::sqrt(p_l*p_l + 2.0*m*(c - c_l)) - p_l)/m; + } + + // Determine outgoing photon energy corresponding to electron momentum + // (solve Eq. 39 in LA-UR-04-0487 for E') + double momentum_sq = std::pow((pz/FINE_STRUCTURE), 2); + double f = 1.0 + alpha*(1.0 - mu); + double a = momentum_sq - f*f; + double b = 2.0*E*(f - momentum_sq*mu); + c = E*E*(momentum_sq - 1.0); + + double quad = b*b - 4.0*a*c; + if (quad < 0) { + *E_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON_EV; + break; + } + quad = std::sqrt(quad); + double E_out1 = -(b + quad)/(2.0*a); + double E_out2 = -(b - quad)/(2.0*a); + + // Determine solution to quadratic equation that is positive + if (E_out1 > 0.0) { + if (E_out2 > 0.0) { + // If both are positive, pick one at random + *E_out = prn() < 0.5 ? E_out1 : E_out2; + } else { + *E_out = E_out1; + } + } else { + if (E_out2 > 0.0) { + *E_out = E_out2; + } else { + // No positive solution -- resample + continue; + } + } + if (*E_out < E - E_b) break; + } + + *i_shell = shell; +} + +double PhotonInteraction::rayleigh_scatter(double alpha) const +{ + double mu; + while (true) { + // Determine maximum value of x^2 + double x2_max = std::pow(MASS_ELECTRON_EV/PLANCK_C*alpha, 2); + + // Determine F(x^2_max, Z) + double F_max = coherent_int_form_factor_(x2_max); + + // Sample cumulative distribution + double F = prn()*F_max; + + // Determine x^2 corresponding to F + const auto& x {coherent_int_form_factor_.x()}; + const auto& y {coherent_int_form_factor_.y()}; + int i = lower_bound_index(y.cbegin(), y.cend(), F); + double r = (F - y[i]) / (y[i+1] - y[i]); + double x2 = x[i] + r*(x[i+1] - x[i]); + + // Calculate mu + mu = 1.0 - 2.0*x2/x2_max; + + if (prn() < 0.5*(1.0 + mu*mu)) break; + } + return mu; +} + +void PhotonInteraction::pair_production(double alpha, double* E_electron, + double* E_positron, double* mu_electron, double* mu_positron) const +{ + constexpr double r[] { + 122.81, 73.167, 69.228, 67.301, 64.696, 61.228, + 57.524, 54.033, 50.787, 47.851, 46.373, 45.401, + 44.503, 43.815, 43.074, 42.321, 41.586, 40.953, + 40.524, 40.256, 39.756, 39.144, 38.462, 37.778, + 37.174, 36.663, 35.986, 35.317, 34.688, 34.197, + 33.786, 33.422, 33.068, 32.740, 32.438, 32.143, + 31.884, 31.622, 31.438, 31.142, 30.950, 30.758, + 30.561, 30.285, 30.097, 29.832, 29.581, 29.411, + 29.247, 29.085, 28.930, 28.721, 28.580, 28.442, + 28.312, 28.139, 27.973, 27.819, 27.675, 27.496, + 27.285, 27.093, 26.911, 26.705, 26.516, 26.304, + 26.108, 25.929, 25.730, 25.577, 25.403, 25.245, + 25.100, 24.941, 24.790, 24.655, 24.506, 24.391, + 24.262, 24.145, 24.039, 23.922, 23.813, 23.712, + 23.621, 23.523, 23.430, 23.331, 23.238, 23.139, + 23.048, 22.967, 22.833, 22.694, 22.624, 22.545, + 22.446, 22.358, 22.264}; + + // The reduced screening radius r is the ratio of the screening radius to + // the Compton wavelength of the electron, where the screening radius is + // obtained under the assumption that the Coulomb field of the nucleus is + // exponentially screened by atomic electrons. This allows us to use a + // simplified atomic form factor and analytical approximations of the + // screening functions in the pair production DCS instead of computing the + // screening functions numerically. The reduced screening radii above for + // Z = 1-99 come from 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). + + // Compute the minimum and maximum values of the electron reduced energy, + // i.e. the fraction of the photon energy that is given to the electron + double e_min = 1.0/alpha; + double e_max = 1.0 - 1.0/alpha; + + // Compute the high-energy Coulomb correction + double a = Z_ / FINE_STRUCTURE; + double c = a*a*(1.0/(1.0 + a*a) + 0.202059 + a*a*(-0.03693 + a*a*(0.00835 + + a*a*(-0.00201 + a*a*(0.00049 + a*a*(-0.00012 + a*a*0.00003)))))); + + // The analytical approximation of the DCS underestimates the cross section + // at low energies. The correction factor f compensates for this. + double q = std::sqrt(2.0/alpha); + double f = q*(-0.1774 - 12.10*a + 11.18*a*a) + + q*q*(8.523 + 73.26*a - 44.41*a*a) + + q*q*q*(-13.52 - 121.1*a + 96.41*a*a) + + q*q*q*q*(8.946 + 62.05*a - 63.41*a*a); + + // Calculate phi_1(1/2) and phi_2(1/2). The unnormalized PDF for the reduced + // energy is given by p = 2*(1/2 - e)^2*phi_1(e) + phi_2(e), where phi_1 and + // phi_2 are non-negative and maximum at e = 1/2. + double b = 2.0*r[Z_]/alpha; + double t1 = 2.0*std::log(1.0 + b*b); + double t2 = b*std::atan(1.0/b); + double t3 = b*b*(4.0 - 4.0*t2 - 3.0*std::log(1.0 + 1.0/(b*b))); + double t4 = 4.0*std::log(r[Z_]) - 4.0*c + f; + double phi1_max = 7.0/3.0 - t1 - 6.0*t2 - t3 + t4; + double phi2_max = 11.0/6.0 - t1 - 3.0*t2 + 0.5*t3 + t4; + + // To aid sampling, the unnormalized PDF can be expressed as + // p = u_1*U_1(e)*pi_1(e) + u_2*U_2(e)*pi_2(e), where pi_1 and pi_2 are + // normalized PDFs on the interval (e_min, e_max) from which values of e can + // be sampled using the inverse transform method, and + // U_1 = phi_1(e)/phi_1(1/2) and U_2 = phi_2(e)/phi_2(1/2) are valid + // rejection functions. The reduced energy can now be sampled using a + // combination of the composition and rejection methods. + double u1 = 2.0/3.0*std::pow(0.5 - 1.0/alpha, 2)*phi1_max; + double u2 = phi2_max; + double e; + while (true) { + double rn = prn(); + + // Sample the index i in (1, 2) using the point probabilities + // p(1) = u_1/(u_1 + u_2) and p(2) = u_2/(u_1 + u_2) + int i; + if (prn() < u1/(u1 + u2)) { + i = 1; + + // Sample e from pi_1 using the inverse transform method + e = rn >= 0.5 ? + 0.5 + (0.5 - 1.0/alpha)*std::pow(2.0*rn - 1.0, 1.0/3.0) : + 0.5 - (0.5 - 1.0/alpha)*std::pow(1.0 - 2.0*rn, 1.0/3.0); + } else { + i = 2; + + // Sample e from pi_2 using the inverse transform method + e = 1.0/alpha + (0.5 - 1.0/alpha)*2.0*rn; + } + + // Calculate phi_i(e) and deliver e if rn <= U_i(e) + b = r[Z_]/(2.0*alpha*e*(1.0 - e)); + t1 = 2.0*std::log(1.0 + b*b); + t2 = b*std::atan(1.0/b); + t3 = b*b*(4.0 - 4.0*t2 - 3.0*std::log(1.0 + 1.0/(b*b))); + if (i == 1) { + double phi1 = 7.0/3.0 - t1 - 6.0*t2 - t3 + t4; + if (prn() <= phi1/phi1_max) break; + } else { + double phi2 = 11.0/6.0 - t1 - 3.0*t2 + 0.5*t3 + t4; + if (prn() <= phi2/phi2_max) break; + } + } + + // Compute the kinetic energy of the electron and the positron + *E_electron = (alpha*e - 1.0)*MASS_ELECTRON_EV; + *E_positron = (alpha*(1.0 - e) - 1.0)*MASS_ELECTRON_EV; + + // Sample the scattering angle of the electron. The cosine of the polar + // angle of the direction relative to the incident photon is sampled from + // p(mu) = C/(1 - beta*mu)^2 using the inverse transform method. + double beta = std::sqrt(*E_electron*(*E_electron + 2.0*MASS_ELECTRON_EV)) + / (*E_electron + MASS_ELECTRON_EV) ; + double rn = 2.0*prn() - 1.0; + *mu_electron = (rn + beta)/(rn*beta + 1.0); + + // Sample the scattering angle of the positron + beta = std::sqrt(*E_positron*(*E_positron + 2.0*MASS_ELECTRON_EV)) + / (*E_positron + MASS_ELECTRON_EV); + rn = 2.0*prn() - 1.0; + *mu_positron = (rn + beta)/(rn*beta + 1.0); +} + +void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particle& p) const +{ + // If no transitions, assume fluorescent photon from captured free electron + if (shell.n_transitions == 0) { + double mu = 2.0*prn() - 1.0; + double phi = 2.0*PI*prn(); + std::array uvw; + uvw[0] = mu; + uvw[1] = std::sqrt(1.0 - mu*mu)*std::cos(phi); + uvw[2] = std::sqrt(1.0 - mu*mu)*std::sin(phi); + double E = shell.binding_energy; + int photon = static_cast(ParticleType::photon); + p.create_secondary(uvw.data(), E, photon, true); + return; + } + + // Sample transition + double rn = prn(); + double c = 0.0; + int i_transition; + for (i_transition = 0; i_transition < shell.n_transitions; ++i_transition) { + c += shell.transition_probability(i_transition); + if (rn < c) break; + } + + // Get primary and secondary subshell designators + int primary = shell.transition_subshells(i_transition, 0); + int secondary = shell.transition_subshells(i_transition, 1); + + // Sample angle isotropically + double mu = 2.0*prn() - 1.0; + double phi = 2.0*PI*prn(); + std::array uvw; + uvw[0] = mu; + uvw[1] = std::sqrt(1.0 - mu*mu)*std::cos(phi); + uvw[2] = std::sqrt(1.0 - mu*mu)*std::sin(phi); + + // Get the transition energy + double E = shell.transition_energy(i_transition); + + if (secondary != 0) { + // Non-radiative transition -- Auger/Coster-Kronig effect + + // Create auger electron + int electron = static_cast(ParticleType::electron); + p.create_secondary(uvw.data(), E, electron, true); + + // Fill hole left by emitted auger electron + int i_hole = shell_map_.at(secondary); + const auto& hole = shells_[i_hole]; + this->atomic_relaxation(hole, p); + } else { + // Radiative transition -- get X-ray energy + + // Create fluorescent photon + int photon = static_cast(ParticleType::photon); + p.create_secondary(uvw.data(), E, photon, true); + } + + // Fill hole created by electron transitioning to the photoelectron hole + int i_hole = shell_map_.at(primary); + const auto& hole = shells_[i_hole]; + this->atomic_relaxation(hole, p); +} + +//============================================================================== +// Non-member functions +//============================================================================== + +std::pair klein_nishina(double alpha) +{ + double alpha_out, mu; + double beta = 1.0 + 2.0*alpha; + if (alpha < 3.0) { + // Kahn's rejection method + double t = beta/(beta + 8.0); + double x; + while (true) { + if (prn() < t) { + // Left branch of flow chart + double r = 2.0*prn(); + x = 1.0 + alpha*r; + if (prn() < 4.0/x*(1.0 - 1.0/x)) { + mu = 1 - r; + break; + } + } else { + // Right branch of flow chart + x = beta/(1.0 + 2.0*alpha*prn()); + mu = 1.0 + (1.0 - x)/alpha; + if (prn() < 0.5*(mu*mu + 1.0/x)) break; + } + } + alpha_out = alpha/x; + + } else { + // Koblinger's direct method + double gamma = 1.0 - std::pow(beta, -2); + double s = prn()*(4.0/alpha + 0.5*gamma + + (1.0 - (1.0 + beta)/(alpha*alpha))*std::log(beta)); + if (s <= 2.0/alpha) { + // For first term, x = 1 + 2ar + // Therefore, a' = a/(1 + 2ar) + alpha_out = alpha/(1.0 + 2.0*alpha*prn()); + } else if (s <= 4.0/alpha) { + // For third term, x = beta/(1 + 2ar) + // Therefore, a' = a(1 + 2ar)/beta + alpha_out = alpha*(1.0 + 2.0*alpha*prn())/beta; + } else if (s <= 4.0/alpha + 0.5*gamma) { + // For fourth term, x = 1/sqrt(1 - gamma*r) + // Therefore, a' = a*sqrt(1 - gamma*r) + alpha_out = alpha*std::sqrt(1.0 - gamma*prn()); + } else { + // For third term, x = beta^r + // Therefore, a' = a/beta^r + alpha_out = alpha/std::pow(beta, prn()); + } + + // Calculate cosine of scattering angle based on basic relation + mu = 1.0 + 1.0/alpha - 1.0/alpha_out; + } + 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 //============================================================================== diff --git a/src/photon_physics.F90 b/src/photon_physics.F90 index 1a9e4edea..157545074 100644 --- a/src/photon_physics.F90 +++ b/src/photon_physics.F90 @@ -3,539 +3,12 @@ module photon_physics use algorithm, only: binary_search use constants use particle_header - use photon_header, only: PhotonInteraction, BremsstrahlungData, & - compton_profile_pz, ttb_e_grid, ttb + use photon_header, only: BremsstrahlungData, ttb_e_grid, ttb use random_lcg, only: prn use settings contains -!=============================================================================== -! KLEIN_NISHINA -!=============================================================================== - - subroutine klein_nishina(alpha, alpha_out, mu) - real(8), intent(in) :: alpha - real(8), intent(out) :: alpha_out - real(8), intent(out) :: mu - - real(8) :: beta ! 1 + 2a - real(8) :: t ! (1 + 2a)/(9 + 2a) - real(8) :: r, s, x - real(8) :: gamma - - beta = ONE + TWO*alpha - if (alpha < THREE) then - ! Kahn's rejection method - t = beta/(beta + 8.0_8) - do - if (prn() < t) then - ! Left branch of flow chart - r = TWO*prn() - x = ONE + alpha*r - if (prn() < FOUR/x*(ONE - ONE/x)) then - mu = 1 - r - exit - end if - else - ! Right branch of flow chart - x = beta/(ONE + TWO*alpha*prn()) - mu = ONE + (ONE - x)/alpha - if (prn() < HALF*(mu**2 + ONE/x)) exit - end if - end do - alpha_out = alpha/x - - else - ! Koblinger's direct method - gamma = ONE - beta**(-2) - s = prn()*(FOUR/alpha + HALF*gamma + & - (ONE - (ONE + beta)/alpha**2)*log(beta)) - if (s <= 2./alpha) then - ! For first term, x = 1 + 2ar - ! Therefore, a' = a/(1 + 2ar) - alpha_out = alpha/(ONE + TWO*alpha*prn()) - elseif (s <= FOUR/alpha) then - ! For third term, x = beta/(1 + 2ar) - ! Therefore, a' = a(1 + 2ar)/beta - alpha_out = alpha*(ONE + TWO*alpha*prn())/beta - elseif (s <= FOUR/alpha + HALF*gamma) then - ! For fourth term, x = 1/sqrt(1 - gamma*r) - ! Therefore, a' = a*sqrt(1 - gamma*r) - alpha_out = alpha*sqrt(ONE - gamma*prn()) - else - ! For third term, x = beta^r - ! Therefore, a' = a/beta^r - alpha_out = alpha/beta**prn() - end if - - ! Calculate cosine of scattering angle based on basic relation - mu = ONE + ONE/alpha - ONE/alpha_out - end if - - end subroutine klein_nishina - -!=============================================================================== -! COMPTON_SCATTER -!=============================================================================== - - subroutine compton_scatter(el, alpha, alpha_out, mu, i_shell, use_doppler) - type(PhotonInteraction), intent(in) :: el - real(8), intent(in) :: alpha - real(8), intent(out) :: alpha_out - real(8), intent(out) :: mu - integer, intent(out) :: i_shell - logical, intent(in), optional :: use_doppler - - real(8) :: x - real(8) :: form_factor_xmax - real(8) :: form_factor_x - real(8) :: e_out - logical :: use_doppler_ - - if (present(use_doppler)) then - use_doppler_ = use_doppler - else - use_doppler_ = .false. - end if - - form_factor_xmax = ZERO - do - ! Sample Klein-Nishina distribution for trial energy and angle - call klein_nishina(alpha, alpha_out, mu) - - ! Note that the parameter used here does not correspond exactly to the - ! momentum transfer q in ENDF-102 Eq. (27.2). Rather, this is the - ! parameter as defined by Hubbell, where the actual data comes from - x = MASS_ELECTRON_EV/PLANCK_C*alpha*sqrt(HALF*(ONE - mu)) - - ! Calculate S(x, Z) and S(x_max, Z) - form_factor_x = el % incoherent_form_factor % evaluate(x) - if (form_factor_xmax == ZERO) then - form_factor_xmax = el % incoherent_form_factor % evaluate(& - MASS_ELECTRON_EV/PLANCK_C*alpha) - end if - - ! Perform rejection on form factor - if (prn() < form_factor_x / form_factor_xmax) then - if (use_doppler_) then - call compton_doppler(el, alpha, mu, e_out, i_shell) - alpha_out = e_out/MASS_ELECTRON_EV - else - i_shell = 0 - end if - exit - end if - end do - - end subroutine compton_scatter - -!=============================================================================== -! COMPTON_DOPPLER -!=============================================================================== - - subroutine compton_doppler(el, alpha, mu, e_out, i_shell) - type(PhotonInteraction), intent(in) :: el - real(8), intent(in) :: alpha - real(8), intent(in) :: mu - real(8), intent(out) :: e_out - integer, intent(out) :: i_shell - - integer :: i - integer :: n - real(8) :: rn, m - real(8) :: c, c_l, c_max - real(8) :: pz_l, pz_r, pz, pz_max - real(8) :: p_l, p_r - real(8) :: e, e_b - real(8) :: e_out1, e_out2 - real(8) :: a, b, quad - real(8) :: f - real(8) :: momentum_sq - - n = size(compton_profile_pz) - - do - ! Sample electron shell - rn = prn() - c = ZERO - do i_shell = 1, size(el % electron_pdf) - c = c + el % electron_pdf(i_shell) - if (rn < c) exit - end do - - ! Determine binding energy of shell - e_b = el % binding_energy(i_shell) - - ! Determine p_z,max - e = alpha*MASS_ELECTRON_EV - if (e < e_b) then - e_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON_EV - exit - end if - - pz_max = -FINE_STRUCTURE*(e_b - (e - e_b)*alpha*(ONE - mu)) / & - sqrt(TWO*e*(e - e_b)*(ONE - mu) + e_b**2) - if (pz_max < ZERO) then - e_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON_EV - exit - end if - - ! Determine profile cdf value corresponding to p_z,max - if (pz_max > compton_profile_pz(n)) then - c_max = el % profile_cdf(n, i_shell) - else - i = binary_search(compton_profile_pz, n, pz_max) - pz_l = compton_profile_pz(i) - pz_r = compton_profile_pz(i + 1) - p_l = el % profile_pdf(i, i_shell) - p_r = el % profile_pdf(i + 1, i_shell) - c_l = el % profile_cdf(i, i_shell) - if (pz_l == pz_r) then - c_max = c_l - elseif (p_l == p_r) then - c_max = c_l + (pz_max - pz_l)*p_l - else - m = (p_l - p_r)/(pz_l - pz_r) - c_max = c_l + ((m*(pz_max - pz_l) + p_l)**2 - p_l**2)/(TWO*m) - end if - end if - - ! Sample value on bounded cdf - c = prn()*c_max - - ! Determine pz corresponding to sampled cdf value - i = binary_search(el % profile_cdf(:, i_shell), n, c) - pz_l = compton_profile_pz(i) - pz_r = compton_profile_pz(i + 1) - p_l = el % profile_pdf(i, i_shell) - p_r = el % profile_pdf(i + 1, i_shell) - c_l = el % profile_cdf(i, i_shell) - if (pz_l == pz_r) then - pz = pz_l - elseif (p_l == p_r) then - pz = pz_l + (c - c_l)/p_l - else - m = (p_l - p_r)/(pz_l - pz_r) - pz = pz_l + (sqrt(p_l**2 + TWO*m*(c - c_l)) - p_l)/m - end if - - ! Determine outgoing photon energy corresponding to electron momentum - ! (solve Eq. 39 in LA-UR-04-0487 for E') - momentum_sq = (pz/FINE_STRUCTURE)**2 - f = ONE + alpha*(ONE - mu) - a = momentum_sq - f*f - b = TWO*e*(f - momentum_sq*mu) - c = e**2*(momentum_sq - ONE) - - quad = b**2 - FOUR*a*c - if (quad < 0) then - e_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON_EV - exit - end if - quad = sqrt(quad) - e_out1 = -(b + quad)/(TWO*a) - e_out2 = -(b - quad)/(TWO*a) - - ! Determine solution to quadratic equation that is positive - if (e_out1 > ZERO) then - if (e_out2 > ZERO) then - ! If both are positive, pick one at random - if (prn() < HALF) then - e_out = e_out1 - else - e_out = e_out2 - end if - else - e_out = e_out1 - end if - else - if (e_out2 > ZERO) then - e_out = e_out2 - else - ! No positive solution -- resample - cycle - end if - end if - if (e_out < e - e_b) exit - end do - - end subroutine compton_doppler - -!=============================================================================== -! RAYLEIGH_SCATTER -!=============================================================================== - - subroutine rayleigh_scatter(el, alpha, mu) - type(PhotonInteraction), intent(in) :: el - real(8), intent(in) :: alpha - real(8), intent(out) :: mu - - integer :: i - real(8) :: F - real(8) :: F_max - real(8) :: x2 - real(8) :: x2_max - real(8) :: r - - do - ! Determine maximum value of x^2 - x2_max = (MASS_ELECTRON_EV/PLANCK_C*alpha)**2 - - ! Determine F(x^2_max, Z) - F_max = el % coherent_int_form_factor % evaluate(x2_max) - - ! Sample cumulative distribution - F = prn()*F_max - - ! Determine x^2 corresponding to F - i = binary_search(el%coherent_int_form_factor%y, & - size(el%coherent_int_form_factor%y), F) - r = (F - el%coherent_int_form_factor%y(i)) / & - (el%coherent_int_form_factor%y(i+1) - el%coherent_int_form_factor%y(i)) - x2 = el%coherent_int_form_factor%x(i) + r*(el%coherent_int_form_factor%x(i+1) - & - el%coherent_int_form_factor%x(i)) - - ! Calculate mu - mu = ONE - TWO*x2/x2_max - - if (prn() < HALF*(ONE + mu**2)) exit - end do - - end subroutine rayleigh_scatter - -!=============================================================================== -! ATOMIC_RELAXATION -!=============================================================================== - - recursive subroutine atomic_relaxation(p, elm, i_shell) - type(Particle), intent(inout) :: p - type(PhotonInteraction), intent(in) :: elm - integer, intent(in) :: i_shell - - integer :: i_hole - integer :: i_transition - integer :: primary - integer :: secondary - real(8) :: c - real(8) :: rn - real(8) :: E - real(8) :: mu - real(8) :: phi - real(8) :: uvw(3) - - ! If no transitions, assume fluorescent photon from captured free electron - if (elm % shells(i_shell) % n_transitions == 0) then - mu = TWO*prn() - ONE - phi = TWO*PI*prn() - uvw(1) = mu - uvw(2) = sqrt(ONE - mu*mu)*cos(phi) - uvw(3) = sqrt(ONE - mu*mu)*sin(phi) - E = elm % shells(i_shell) % binding_energy - call particle_create_secondary(p, uvw, E, PHOTON, run_ce=.true._C_BOOL) - return - end if - - ! Sample transition - rn = prn() - c = ZERO - do i_transition = 1, elm % shells(i_shell) % n_transitions - c = c + elm % shells(i_shell) % & - transition_probability(i_transition) - if (rn < c) exit - end do - - ! Get primary and secondary subshell designators - primary = elm % shells(i_shell) % transition_subshells(1, i_transition) - secondary = elm % shells(i_shell) % transition_subshells(2, i_transition) - - ! Sample angle isotropically - mu = TWO*prn() - ONE - phi = TWO*PI*prn() - uvw(1) = mu - uvw(2) = sqrt(ONE - mu*mu)*cos(phi) - uvw(3) = sqrt(ONE - mu*mu)*sin(phi) - - ! Get the transition energy - E = elm % shells(i_shell) % transition_energy(i_transition) - - if (secondary /= 0) then - ! Non-radiative transition -- Auger/Coster-Kronig effect - - ! Create auger electron - call particle_create_secondary(p, uvw, E, ELECTRON, run_ce=.true._C_BOOL) - - ! Fill hole left by emitted auger electron - i_hole = elm % shell_dict % get(secondary) - call atomic_relaxation(p, elm, i_hole) - else - ! Radiative transition -- get X-ray energy - - ! Create fluorescent photon - call particle_create_secondary(p, uvw, E, PHOTON, run_ce=.true._C_BOOL) - - end if - - ! Fill hole created by electron transitioning to the photoelectron hole - i_hole = elm % shell_dict % get(primary) - call atomic_relaxation(p, elm, i_hole) - - end subroutine atomic_relaxation - -!=============================================================================== -! PAIR_PRODUCTION samples the kinetic energy and direction of the electron and -! positron created when a photon is absorbed near an atomic nucleus. The -! simulation procedure follows the semiempirical model outlined in 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). -!=============================================================================== - - subroutine pair_production(elm, alpha, E_electron, E_positron, mu_electron, & - mu_positron) - type(PhotonInteraction), intent(in) :: elm - real(8), intent(in) :: alpha - real(8), intent(out) :: E_electron - real(8), intent(out) :: E_positron - real(8), intent(out) :: mu_electron - real(8), intent(out) :: mu_positron - - integer :: i - real(8) :: f - real(8) :: c - real(8) :: a - real(8) :: b - real(8) :: q - real(8) :: rn - real(8) :: beta - real(8) :: e, e_min, e_max - real(8) :: t1, t2, t3, t4 - real(8) :: u1, u2 - real(8) :: phi1, phi2 - real(8) :: phi1_max, phi2_max - real(8), parameter :: r(99) = (/ & - 122.81_8, 73.167_8, 69.228_8, 67.301_8, 64.696_8, 61.228_8, & - 57.524_8, 54.033_8, 50.787_8, 47.851_8, 46.373_8, 45.401_8, & - 44.503_8, 43.815_8, 43.074_8, 42.321_8, 41.586_8, 40.953_8, & - 40.524_8, 40.256_8, 39.756_8, 39.144_8, 38.462_8, 37.778_8, & - 37.174_8, 36.663_8, 35.986_8, 35.317_8, 34.688_8, 34.197_8, & - 33.786_8, 33.422_8, 33.068_8, 32.740_8, 32.438_8, 32.143_8, & - 31.884_8, 31.622_8, 31.438_8, 31.142_8, 30.950_8, 30.758_8, & - 30.561_8, 30.285_8, 30.097_8, 29.832_8, 29.581_8, 29.411_8, & - 29.247_8, 29.085_8, 28.930_8, 28.721_8, 28.580_8, 28.442_8, & - 28.312_8, 28.139_8, 27.973_8, 27.819_8, 27.675_8, 27.496_8, & - 27.285_8, 27.093_8, 26.911_8, 26.705_8, 26.516_8, 26.304_8, & - 26.108_8, 25.929_8, 25.730_8, 25.577_8, 25.403_8, 25.245_8, & - 25.100_8, 24.941_8, 24.790_8, 24.655_8, 24.506_8, 24.391_8, & - 24.262_8, 24.145_8, 24.039_8, 23.922_8, 23.813_8, 23.712_8, & - 23.621_8, 23.523_8, 23.430_8, 23.331_8, 23.238_8, 23.139_8, & - 23.048_8, 22.967_8, 22.833_8, 22.694_8, 22.624_8, 22.545_8, & - 22.446_8, 22.358_8, 22.264_8 /) - - ! The reduced screening radius r is the ratio of the screening radius to - ! the Compton wavelength of the electron, where the screening radius is - ! obtained under the assumption that the Coulomb field of the nucleus is - ! exponentially screened by atomic electrons. This allows us to use a - ! simplified atomic form factor and analytical approximations of the - ! screening functions in the pair production DCS instead of computing the - ! screening functions numerically. The reduced screening radii above for - ! Z = 1-99 come from 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). - - ! Compute the minimum and maximum values of the electron reduced energy, - ! i.e. the fraction of the photon energy that is given to the electron - e_min = ONE/alpha - e_max = ONE - ONE/alpha - - ! Compute the high-energy Coulomb correction - a = elm % Z / FINE_STRUCTURE - c = a**2*(ONE/(ONE + a**2) + 0.202059_8 - 0.03693_8*a**2 + 0.00835_8*a**4 & - - 0.00201_8*a**6 + 0.00049_8*a**8 - 0.00012_8*a**10 + 0.00003_8*a**12) - - ! The analytical approximation of the DCS underestimates the cross section - ! at low energies. The correction factor f compensates for this. - q = sqrt(TWO/alpha) - f = q*(-0.1774_8 - 12.10_8*a + 11.18_8*a**2) & - + q**2*(8.523_8 + 73.26_8*a - 44.41_8*a**2) & - + q**3*(-13.52_8 - 121.1_8*a + 96.41_8*a**2) & - + q**4*(8.946_8 + 62.05_8*a - 63.41_8*a**2) - - ! Calculate phi_1(1/2) and phi_2(1/2). The unnormalized PDF for the reduced - ! energy is given by p = 2*(1/2 - e)^2*phi_1(e) + phi_2(e), where phi_1 and - ! phi_2 are non-negative and maximum at e = 1/2. - b = TWO*r(elm % Z)/alpha - t1 = TWO*log(ONE + b**2) - t2 = b*atan(ONE/b) - t3 = b**2*(FOUR - FOUR*t2 - THREE*log(ONE + ONE/b**2)) - t4 = FOUR*log(r(elm % Z)) - FOUR*c + f - phi1_max = 7.0_8/THREE - t1 - 6.0_8*t2 - t3 + t4 - phi2_max = 11.0_8/6.0_8 - t1 - THREE*t2 + HALF*t3 + t4 - - ! To aid sampling, the unnormalized PDF can be expressed as - ! p = u_1*U_1(e)*pi_1(e) + u_2*U_2(e)*pi_2(e), where pi_1 and pi_2 are - ! normalized PDFs on the interval (e_min, e_max) from which values of e can - ! be sampled using the inverse transform method, and - ! U_1 = phi_1(e)/phi_1(1/2) and U_2 = phi_2(e)/phi_2(1/2) are valid - ! rejection functions. The reduced energy can now be sampled using a - ! combination of the composition and rejection methods. - u1 = TWO/THREE*(HALF - ONE/alpha)**2*phi1_max - u2 = phi2_max - do - rn = prn() - - ! Sample the index i in (1, 2) using the point probabilities - ! p(1) = u_1/(u_1 + u_2) and p(2) = u_2/(u_1 + u_2) - if (prn() < u1/(u1 + u2)) then - i = 1 - - ! Sample e from pi_1 using the inverse transform method - if (rn >= HALF) then - e = HALF + (HALF - ONE/alpha)*(TWO*rn - ONE)**(ONE/THREE) - else - e = HALF - (HALF - ONE/alpha)*(ONE - TWO*rn)**(ONE/THREE) - end if - else - i = 2 - - ! Sample e from pi_2 using the inverse transform method - e = ONE/alpha + (HALF - ONE/alpha)*TWO*rn - end if - - ! Calculate phi_i(e) and deliver e if rn <= U_i(e) - b = r(elm % Z)/(TWO*alpha*e*(ONE - e)) - t1 = TWO*log(ONE + b**2) - t2 = b*atan(ONE/b) - t3 = b**2*(FOUR - FOUR*t2 - THREE*log(ONE + ONE/b**2)) - if (i == 1) then - phi1 = 7.0_8/THREE - t1 - 6.0_8*t2 - t3 + t4 - if (prn() <= phi1/phi1_max) exit - else - phi2 = 11.0_8/6.0_8 - t1 - THREE*t2 + HALF*t3 + t4 - if (prn() <= phi2/phi2_max) exit - end if - end do - - ! Compute the kinetic energy of the electron and the positron - E_electron = (alpha*e - ONE)*MASS_ELECTRON_EV - E_positron = (alpha*(ONE - e) - ONE)*MASS_ELECTRON_EV - - ! Sample the scattering angle of the electron. The cosine of the polar - ! angle of the direction relative to the incident photon is sampled from - ! p(mu) = C/(1 - beta*mu)^2 using the inverse transform method. - beta = sqrt(E_electron*(E_electron + TWO*MASS_ELECTRON_EV)) & - / (E_electron + MASS_ELECTRON_EV) - rn = TWO*prn() - ONE - mu_electron = (rn + beta)/(rn*beta + ONE) - - ! Sample the scattering angle of the positron - beta = sqrt(E_positron*(E_positron + TWO*MASS_ELECTRON_EV)) & - / (E_positron + MASS_ELECTRON_EV) - rn = TWO*prn() - ONE - mu_positron = (rn + beta)/(rn*beta + ONE) - - end subroutine pair_production - !=============================================================================== ! THICK_TARGET_BREMSSTRAHLUNG !=============================================================================== diff --git a/src/physics.F90 b/src/physics.F90 index d77fb04df..c6749a6a1 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -1,21 +1,16 @@ module physics use constants - use error, only: fatal_error, warning, write_message + use error, only: fatal_error use material_header, only: Material, materials use math use message_passing use nuclide_header use particle_header use photon_header - use photon_physics, only: rayleigh_scatter, compton_scatter, & - atomic_relaxation, pair_production, & - thick_target_bremsstrahlung - use physics_common use random_lcg, only: prn use settings use simulation_header - use tally_header implicit none @@ -28,190 +23,13 @@ module physics contains -!=============================================================================== -! SAMPLE_PHOTON_REACTION samples an element based on the macroscopic cross -! sections for each nuclide within a material and then samples a reaction for -! that element and calls the appropriate routine to process the physics. -!=============================================================================== - - subroutine sample_photon_reaction(p) bind(C) - type(Particle), intent(inout) :: p - - integer :: i_shell ! index in subshells - integer :: i_grid ! index on energy grid - integer :: i_element ! index in nuclides array - integer :: i_start ! threshold index - real(8) :: prob ! cumulative probability - real(8) :: cutoff ! sampled total cross section - real(8) :: f ! interpolation factor - real(8) :: xs ! photoionization cross section - real(8) :: r ! random number - real(8) :: prob_after - real(8) :: alpha ! photon energy divided by electron rest mass - real(8) :: alpha_out ! outgoing photon energy over electron rest mass - real(8) :: mu ! scattering cosine - real(8) :: mu_electron ! electron scattering cosine - real(8) :: mu_positron ! positron scattering cosine - real(8) :: phi ! azimuthal angle - real(8) :: uvw(3) ! new direction - real(8) :: rel_vel ! relative velocity of electron - real(8) :: e_b ! binding energy of electron - real(8) :: E_electron ! electron energy - real(8) :: E_positron ! positron energy - - ! Kill photon if below energy cutoff -- an extra check is made here because - ! photons with energy below the cutoff may have been produced by neutrons - ! reactions or atomic relaxation - if (p % E < energy_cutoff(PHOTON)) then - p % E = ZERO - p % alive = .false. - return - end if - - ! Sample element within material - i_element = sample_element(p) - p % event_nuclide = i_element - - ! Calculate photon energy over electron rest mass equivalent - alpha = p % E/MASS_ELECTRON_EV - - ! For tallying purposes, this routine might be called directly. In that - ! case, we need to sample a reaction via the cutoff variable - prob = ZERO - cutoff = prn() * micro_photon_xs(i_element) % total - - associate (elm => elements(i_element)) - ! Coherent (Rayleigh) scattering - prob = prob + micro_photon_xs(i_element) % coherent - if (prob > cutoff) then - call rayleigh_scatter(elm, alpha, mu) - p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu) - p % event_MT = COHERENT - return - end if - - ! Incoherent (Compton) scattering - prob = prob + micro_photon_xs(i_element) % incoherent - if (prob > cutoff) then - call compton_scatter(elm, alpha, alpha_out, mu, i_shell, .true.) - - ! Determine binding energy of shell. The binding energy is zero if - ! doppler broadening is not used. - if (i_shell == 0) then - e_b = ZERO - else - e_b = elm % binding_energy(i_shell) - end if - - ! Create Compton electron - E_electron = (alpha - alpha_out)*MASS_ELECTRON_EV - e_b - mu_electron = (alpha - alpha_out*mu) & - / sqrt(alpha**2 + alpha_out**2 - TWO*alpha*alpha_out*mu) - phi = TWO*PI*prn() - uvw = rotate_angle(p % coord(1) % uvw, mu_electron, phi) - call particle_create_secondary(p, uvw, E_electron, ELECTRON, .true._C_BOOL) - - ! TODO: Compton subshell data does not match atomic relaxation data - ! Allow electrons to fill orbital and produce auger electrons - ! and fluorescent photons - if (i_shell > 0) then - call atomic_relaxation(p, elm, i_shell) - end if - - phi = phi + PI - p % E = alpha_out*MASS_ELECTRON_EV - p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu, phi) - p % event_MT = INCOHERENT - return - end if - - ! Photoelectric effect - prob_after = prob + micro_photon_xs(i_element) % photoelectric - if (prob_after > cutoff) then - do i_shell = 1, size(elm % shells) - ! Get grid index and interpolation factor - i_grid = micro_photon_xs(i_element) % index_grid - f = micro_photon_xs(i_element) % interp_factor - - ! Check threshold of reaction - i_start = elm % shells(i_shell) % threshold - if (i_grid <= i_start) cycle - - ! Evaluation subshell photoionization cross section - xs = exp(elm % shells(i_shell) % cross_section(i_grid - i_start) + & - f*(elm % shells(i_shell) % cross_section(i_grid + 1 - i_start) - & - elm % shells(i_shell) % cross_section(i_grid - i_start))) - - prob = prob + xs - if (prob > cutoff) then - E_electron = p % E - elm % shells(i_shell) % binding_energy - - ! Sample mu using non-relativistic Sauter distribution. - ! See Eqns 3.19 and 3.20 in "Implementing a photon physics - ! model in Serpent 2" by Toni Kaltiaisenaho - SAMPLE_MU: do - r = prn() - if (FOUR * (ONE - r) * r >= prn()) then - rel_vel = sqrt(E_electron * (E_electron + TWO * MASS_ELECTRON_EV))& - / (E_electron + MASS_ELECTRON_EV) - mu = (TWO * r + rel_vel - ONE) / & - (TWO * rel_vel * r - rel_vel + ONE) - exit SAMPLE_MU - end if - end do SAMPLE_MU - - phi = TWO*PI*prn() - uvw(1) = mu - uvw(2) = sqrt(ONE - mu*mu)*cos(phi) - uvw(3) = sqrt(ONE - mu*mu)*sin(phi) - - ! Create secondary electron - call particle_create_secondary(p, uvw, E_electron, ELECTRON, & - run_CE=.true._C_BOOL) - - ! Allow electrons to fill orbital and produce auger electrons - ! and fluorescent photons - call atomic_relaxation(p, elm, i_shell) - p % event_MT = 533 + elm % shells(i_shell) % index_subshell - p % alive = .false. - p % E = ZERO - - return - end if - end do - end if - prob = prob_after - - ! Pair production - prob = prob + micro_photon_xs(i_element) % pair_production - if (prob > cutoff) then - call pair_production(elm, alpha, E_electron, E_positron, mu_electron, & - mu_positron) - - ! Create secondary electron - uvw = rotate_angle(p % coord(1) % uvw, mu_electron) - call particle_create_secondary(p, uvw, E_electron, ELECTRON, .true._C_BOOL) - - ! Create secondary positron - uvw = rotate_angle(p % coord(1) % uvw, mu_positron) - call particle_create_secondary(p, uvw, E_positron, POSITRON, .true._C_BOOL) - - p % event_MT = PAIR_PROD - p % alive = .false. - p % E = ZERO - end if - - end associate - - end subroutine sample_photon_reaction - !=============================================================================== ! SAMPLE_ELEMENT !=============================================================================== - function sample_element(p) result(i_element) + function sample_element(p) result(i_element) bind(C) type(Particle), intent(in) :: p - integer :: i_element + integer(C_INT) :: i_element integer :: i real(8) :: prob diff --git a/src/physics.cpp b/src/physics.cpp index 8aec309e8..a4bd67066 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -214,154 +214,168 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, Bank* bank_ } } -// TODO: Finish converting photon physics functions +void sample_photon_reaction(Particle* p) +{ + // Kill photon if below energy cutoff -- an extra check is made here because + // photons with energy below the cutoff may have been produced by neutrons + // reactions or atomic relaxation + int photon = static_cast(ParticleType::photon) - 1; + if (p->E < settings::energy_cutoff[photon]) { + p->E = 0.0; + p->alive = false; + return; + } -// void sample_photon_reaction(Particle* p) -// { -// // Kill photon if below energy cutoff -- an extra check is made here because -// // photons with energy below the cutoff may have been produced by neutrons -// // reactions or atomic relaxation -// if (p->E < energy_cutoff(PHOTON)) { -// p->E = 0.0 -// p->alive = false -// return -// } + // Sample element within material + int i_element = sample_element(p); + p->event_nuclide = i_element; + // TODO: off-by-one + const auto& micro {simulation::micro_photon_xs[i_element - 1]}; + const auto& element {data::elements[i_element - 1]}; -// // Sample element within material -// i_element = sample_element(p) -// p->event_nuclide = i_element + // Calculate photon energy over electron rest mass equivalent + double alpha = p->E/MASS_ELECTRON_EV; -// // Calculate photon energy over electron rest mass equivalent -// alpha = p->E/MASS_ELECTRON_EV + // For tallying purposes, this routine might be called directly. In that + // case, we need to sample a reaction via the cutoff variable + double prob = 0.0; + double cutoff = prn() * micro.total; -// // For tallying purposes, this routine might be called directly. In that -// // case, we need to sample a reaction via the cutoff variable -// prob = 0.0 -// cutoff = prn() * micro_photon_xs(i_element) % total + // Coherent (Rayleigh) scattering + prob += micro.coherent; + if (prob > cutoff) { + double mu = element.rayleigh_scatter(alpha); + rotate_angle_c(p->coord[0].uvw, mu, nullptr); + p->event_MT = COHERENT; + return; + } -// associate (elm => elements(i_element)) -// // Coherent (Rayleigh) scattering -// prob = prob + micro_photon_xs(i_element) % coherent -// if (prob > cutoff) { -// rayleigh_scatter(elm, alpha, mu) -// p->coord(1) % uvw = rotate_angle(p->coord(1) % uvw, mu) -// p->event_MT = COHERENT -// return -// } + // Incoherent (Compton) scattering + prob += micro.incoherent; + if (prob > cutoff) { + double alpha_out, mu; + int i_shell; + element.compton_scatter(alpha, true, &alpha_out, &mu, &i_shell); -// // Incoherent (Compton) scattering -// prob = prob + micro_photon_xs(i_element) % incoherent -// if (prob > cutoff) { -// compton_scatter(elm, alpha, alpha_out, mu, i_shell, true) + // Determine binding energy of shell. The binding energy is 0.0 if + // doppler broadening is not used. + double e_b; + if (i_shell == -1) { + e_b = 0.0; + } else { + e_b = element.binding_energy_[i_shell]; + } -// // Determine binding energy of shell. The binding energy is 0.0 if -// // doppler broadening is not used. -// if (i_shell == 0) { -// e_b = 0.0 -// } else { -// e_b = elm % binding_energy(i_shell) -// } + // Create Compton electron + double E_electron = (alpha - alpha_out)*MASS_ELECTRON_EV - e_b; + double mu_electron = (alpha - alpha_out*mu) + / std::sqrt(alpha*alpha + alpha_out*alpha_out - 2.0*alpha*alpha_out*mu); + double phi = 2.0*PI*prn(); + double uvw[3]; + std::copy(p->coord[0].uvw, p->coord[0].uvw + 3, uvw); + rotate_angle_c(uvw, mu_electron, &phi); + int electron = static_cast(ParticleType::electron); + p->create_secondary(uvw, E_electron, electron, true); -// // Create Compton electron -// E_electron = (alpha - alpha_out)*MASS_ELECTRON_EV - e_b -// mu_electron = (alpha - alpha_out*mu) & -// / std::sqrt(alpha**2 + alpha_out**2 - 2.0*alpha*alpha_out*mu) -// phi = 2.0*PI*prn() -// uvw = rotate_angle(p->coord(1) % uvw, mu_electron, phi) -// particle_create_secondary(p, uvw, E_electron, ELECTRON, true) + // TODO: Compton subshell data does not match atomic relaxation data + // Allow electrons to fill orbital and produce auger electrons + // and fluorescent photons + if (i_shell >= 0) { + const auto& shell = element.shells_[i_shell]; + element.atomic_relaxation(shell, *p); + } -// // TODO: Compton subshell data does not match atomic relaxation data -// // Allow electrons to fill orbital and produce auger electrons -// // and fluorescent photons -// if (i_shell > 0) { -// atomic_relaxation(p, elm, i_shell) -// } + phi += PI; + p->E = alpha_out*MASS_ELECTRON_EV; + rotate_angle_c(p->coord[0].uvw, mu, &phi); + p->event_MT = INCOHERENT; + return; + } -// phi = phi + PI -// p->E = alpha_out*MASS_ELECTRON_EV -// p->coord(1) % uvw = rotate_angle(p->coord(1) % uvw, mu, phi) -// p->event_MT = INCOHERENT -// return -// } + // Photoelectric effect + double prob_after = prob + micro.photoelectric; + if (prob_after > cutoff) { + for (const auto& shell : element.shells_) { + // Get grid index and interpolation factor + // TODO: off-by-one + int i_grid = micro.index_grid - 1; + double f = micro.interp_factor; -// // Photoelectric effect -// prob_after = prob + micro_photon_xs(i_element) % photoelectric -// if (prob_after > cutoff) { -// do i_shell = 1, size(elm % shells) -// // Get grid index and interpolation factor -// i_grid = micro_photon_xs(i_element) % index_grid -// f = micro_photon_xs(i_element) % interp_factor + // Check threshold of reaction + int i_start = shell.threshold; + if (i_grid <= i_start) continue; -// // Check threshold of reaction -// i_start = elm % shells(i_shell) % threshold -// if (i_grid <= i_start) cycle + // Evaluation subshell photoionization cross section + double xs = std::exp(shell.cross_section(i_grid - i_start) + + f*(shell.cross_section(i_grid + 1 - i_start) - + shell.cross_section(i_grid - i_start))); -// // Evaluation subshell photoionization cross section -// xs = std::exp(elm % shells(i_shell) % cross_section(i_grid - i_start) + & -// f*(elm % shells(i_shell) % cross_section(i_grid + 1 - i_start) - & -// elm % shells(i_shell) % cross_section(i_grid - i_start))) + prob += xs; + if (prob > cutoff) { + double E_electron = p->E - shell.binding_energy; -// prob = prob + xs -// if (prob > cutoff) { -// E_electron = p->E - elm % shells(i_shell) % binding_energy + // Sample mu using non-relativistic Sauter distribution. + // See Eqns 3.19 and 3.20 in "Implementing a photon physics + // model in Serpent 2" by Toni Kaltiaisenaho + double mu; + while (true) { + double r = prn(); + if (4.0*(1.0 - r)*r >= prn()) { + double rel_vel = std::sqrt(E_electron * (E_electron + + 2.0*MASS_ELECTRON_EV)) / (E_electron + MASS_ELECTRON_EV); + mu = (2.0*r + rel_vel - 1.0) / (2.0*rel_vel*r - rel_vel + 1.0); + break; + } + } -// // Sample mu using non-relativistic Sauter distribution. -// // See Eqns 3.19 and 3.20 in "Implementing a photon physics -// // model in Serpent 2" by Toni Kaltiaisenaho -// SAMPLE_MU: do -// r = prn() -// if (FOUR * (1.0 - r) * r >= prn()) { -// rel_vel = std::sqrt(E_electron * (E_electron + 2.0 * MASS_ELECTRON_EV))& -// / (E_electron + MASS_ELECTRON_EV) -// mu = (2.0 * r + rel_vel - 1.0) / & -// (2.0 * rel_vel * r - rel_vel + 1.0) -// exit SAMPLE_MU -// } -// end do SAMPLE_MU + double phi = 2.0*PI*prn(); + std::array uvw; + uvw[0] = mu; + uvw[1] = std::sqrt(1.0 - mu*mu)*std::cos(phi); + uvw[2] = std::sqrt(1.0 - mu*mu)*std::sin(phi); -// phi = 2.0*PI*prn() -// uvw(1) = mu -// uvw(2) = std::sqrt(1.0 - mu*mu)*std::cos(phi) -// uvw(3) = std::sqrt(1.0 - mu*mu)*std::sin(phi) + // Create secondary electron + int electron = static_cast(ParticleType::electron); + p->create_secondary(uvw.data(), E_electron, electron, true); -// // Create secondary electron -// particle_create_secondary(p, uvw, E_electron, ELECTRON, & -// run_CE=true) + // Allow electrons to fill orbital and produce auger electrons + // and fluorescent photons + element.atomic_relaxation(shell, *p); + p->event_MT = 533 + shell.index_subshell; + p->alive = false; + p->E = 0.0; + return; + } + } + } + prob = prob_after; -// // Allow electrons to fill orbital and produce auger electrons -// // and fluorescent photons -// atomic_relaxation(p, elm, i_shell) -// p->event_MT = 533 + elm % shells(i_shell) % index_subshell -// p->alive = false -// p->E = 0.0 + // Pair production + prob += micro.pair_production; + if (prob > cutoff) { + double E_electron, E_positron; + double mu_electron, mu_positron; + element.pair_production(alpha, &E_electron, &E_positron, + &mu_electron, &mu_positron); -// return -// } -// end do -// } -// prob = prob_after + // Create secondary electron + double uvw[3]; + std::copy(p->coord[0].uvw, p->coord[0].uvw + 3, uvw); + rotate_angle_c(uvw, mu_electron, nullptr); + int electron = static_cast(ParticleType::electron); + p->create_secondary(uvw, E_electron, electron, true); -// // Pair production -// prob = prob + micro_photon_xs(i_element) % pair_production -// if (prob > cutoff) { -// pair_production(elm, alpha, E_electron, E_positron, mu_electron, & -// mu_positron) + // Create secondary positron + std::copy(p->coord[0].uvw, p->coord[0].uvw + 3, uvw); + rotate_angle_c(uvw, mu_positron, nullptr); + int positron = static_cast(ParticleType::positron); + p->create_secondary(uvw, E_positron, positron, true); -// // Create secondary electron -// uvw = rotate_angle(p->coord(1) % uvw, mu_electron) -// particle_create_secondary(p, uvw, E_electron, ELECTRON, true) - -// // Create secondary positron -// uvw = rotate_angle(p->coord(1) % uvw, mu_positron) -// particle_create_secondary(p, uvw, E_positron, POSITRON, true) - -// p->event_MT = PAIR_PROD -// p->alive = false -// p->E = 0.0 -// } - -// end associate -// } + p->event_MT = PAIR_PROD; + p->alive = false; + p->E = 0.0; + } +} void sample_electron_reaction(Particle* p) { @@ -465,36 +479,31 @@ void sample_nuclide(const Particle* p, int mt, int* i_nuclide, int* i_nuc_mat) } } -// TODO: Finish converting photon physics functions - -// void sample_element(Particle* p) +// int sample_element(Particle* p) // { -// associate (mat => materials(p->material)) -// // Sample cumulative distribution function -// cutoff = prn() * simulation::material_xs.total +// // Sample cumulative distribution function +// double cutoff = prn() * simulation::material_xs.total; -// i = 0 -// prob = 0.0 -// do while (prob < cutoff) -// i = i + 1 +// int i = 0; +// double prob = 0.0; +// while (prob < cutoff) { +// // Check to make sure that a nuclide was sampled +// if (i > mat % n_nuclides) { +// p->write_restart(); +// fatal_error("Did not sample any element during collision."); +// } -// // Check to make sure that a nuclide was sampled -// if (i > mat % n_nuclides) { -// particle_write_restart(p) -// fatal_error("Did not sample any element during collision.") -// } +// // Find atom density +// int i_element = mat % element(i); +// double atom_density = mat % atom_density(i); -// // Find atom density -// i_element = mat % element(i) -// atom_density = mat % atom_density(i) +// // Determine microscopic cross section +// double sigma = atom_density * simulation::micro_photon_xs[i_element].total; -// // Determine microscopic cross section -// sigma = atom_density * micro_photon_xs(i_element) % total - -// // Increment probability to compare to cutoff -// prob = prob + sigma -// end do -// end associate +// // Increment probability to compare to cutoff +// prob += sigma; +// ++i; +// } // } Reaction* sample_fission(int i_nuclide, double E)