crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-25 12:35:29 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Convert all photon physics except bremsstrahlung (not done yet)

Browse source
This commit is contained in:
Paul Romano 2019-01-04 16:27:45 -06:00
parent 6c3f35d019
commit 2e6ac03433
9 changed files with 717 additions and 872 deletions
Download patch file Download diff file Expand all files Collapse all files

4
include/openmc/endf.h
View file

@ -78,6 +78,10 @@ public:
//! \param[in] x independent variable
//! \return Function evaluated at x
double operator()(double x) const;
// Accessors
const std::vector<double>& x() const { return x_; }
const std::vector<double>& y() const { return y_; }
private:
std::size_t n_regions_ {0}; //!< number of interpolation regions
std::vector<int> nbt_; //!< values separating interpolation regions

2
include/openmc/math_functions.h
View file

@ -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);

24
include/openmc/photon.h
View file

@ -2,12 +2,14 @@
#define OPENMC_PHOTON_H
#include "openmc/endf.h"
#include "openmc/particle.h"
#include <hdf5.h>
#include "xtensor/xtensor.hpp"
#include <string>
#include <unordered_map>
#include <utility> // for pair
#include <vector>
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<double, 2> 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<double, double> klein_nishina(double alpha);
void thick_target_bremsstrahlung(Particle& p, double* E_lost);
//==============================================================================
// Global variables
//==============================================================================

4
include/openmc/physics.h
View file

@ -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);

2
src/math_functions.cpp
View file

@ -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

523
src/photon.cpp
View file

@ -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 <array>
#include <cmath>
#include <tuple> // 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<double, 3> 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<int>(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<double, 3> 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<int>(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<int>(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<double, double> 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<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
//==============================================================================

529
src/photon_physics.F90
View file

@ -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
!===============================================================================

188
src/physics.F90
View file

@ -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

313
src/physics.cpp
View file

@ -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<int>(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<int>(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<double, 3> 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<int>(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<int>(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<int>(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)