mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 14:15:42 -04:00
Merge pull request #1148 from paulromano/cpp-photon
Translate photon physics to C++
This commit is contained in:
Amanda Lund
GitHub
commit
57ed2f8e4b
31 changed files with 1585 additions and 1940 deletions
|
|
@ -17,6 +17,8 @@ addons:
|
|||
- liblapack-dev
|
||||
config:
|
||||
retries: true
|
||||
services:
|
||||
- xvfb
|
||||
cache:
|
||||
directories:
|
||||
- $HOME/nndc_hdf5
|
||||
|
|
@ -31,7 +33,6 @@ env:
|
|||
- OPENMC_ENDF_DATA=$HOME/endf-b-vii.1
|
||||
- LD_LIBRARY_PATH=$HOME/MOAB/lib:$HOME/DAGMC/lib
|
||||
- PATH=$PATH:$HOME/NJOY2016/build
|
||||
- DISPLAY=:99.0
|
||||
- COVERALLS_PARALLEL=true
|
||||
matrix:
|
||||
- OMP=n MPI=n PHDF5=n
|
||||
|
|
|
|||
|
|
@ -331,9 +331,6 @@ add_library(libopenmc SHARED
|
|||
src/particle_header.F90
|
||||
src/particle_restart.F90
|
||||
src/photon_header.F90
|
||||
src/photon_physics.F90
|
||||
src/physics_common.F90
|
||||
src/physics.F90
|
||||
src/pugixml/pugixml_f.F90
|
||||
src/random_lcg.F90
|
||||
src/reaction_header.F90
|
||||
|
|
@ -370,6 +367,7 @@ add_library(libopenmc SHARED
|
|||
src/tallies/trigger.F90
|
||||
src/tallies/trigger_header.F90
|
||||
src/bank.cpp
|
||||
src/bremsstrahlung.cpp
|
||||
src/dagmc.cpp
|
||||
src/cell.cpp
|
||||
src/cmfd_solver.cpp
|
||||
|
|
|
|||
|
|
@ -18,6 +18,11 @@ code:
|
|||
* `Jingang Liang <https://github.com/liangjg>`_
|
||||
* `Colin Josey <https://github.com/cjosey>`_
|
||||
* `Amanda Lund <https://github.com/amandalund>`_
|
||||
* `Guillaume Giudicelli <https://github.com/giudgiud>`_
|
||||
* `Isaac Meyer <https://github.com/icmeyer>`_
|
||||
* `Patrick Shriwise <https://github.com/pshriwise>`_
|
||||
* `Shikhar Kumar <https://github.com/shikhar413>`_
|
||||
* `Andrew Davis <https://github.com/makeclean>`_
|
||||
* `Benoit Forget <http://web.mit.edu/nse/people/faculty/forget.html>`_
|
||||
* `Kord Smith <http://web.mit.edu/nse/people/faculty/smith.html>`_
|
||||
* `Andrew Siegel <http://www.mcs.anl.gov/person/andrew-siegel>`_
|
||||
|
|
|
|||
|
|
@ -57,9 +57,14 @@ Benchmarking
|
|||
Coupling and Multi-physics
|
||||
--------------------------
|
||||
|
||||
- Ze-Long Zhao, Yongwei Yang, and Shuang Hong, "`Application of FLUKA and OpenMC
|
||||
in coupled physics calculation of target and subcritical reactor for ADS
|
||||
<https://doi.org/10.1007/s41365-018-0539-1>`_," *Nucl. Sci. Tech.*, **30**
|
||||
(2019).
|
||||
|
||||
- Jun Chen, Liangzhi Cao, Chuanqi Zhao, and Zhouyu Liu, "`Development of
|
||||
Subchannel Code SUBSC for high-fidelity multi-physics coupling application
|
||||
<https://doi.org/10.1016/j.egypro.2017.08.121>`_", Energy Procedia, **127**,
|
||||
<https://doi.org/10.1016/j.egypro.2017.08.121>`_", *Energy Procedia*, **127**,
|
||||
264-274 (2017).
|
||||
|
||||
- Tianliang Hu, Liangzhu Cao, Hongchun Wu, Xianan Du, and Mingtao He, "`Coupled
|
||||
|
|
|
|||
48
include/openmc/bremsstrahlung.h
Normal file
48
include/openmc/bremsstrahlung.h
Normal file
|
|
@ -0,0 +1,48 @@
|
|||
#ifndef OPENMC_BREMSSTRAHLUNG_H
|
||||
#define OPENMC_BREMSSTRAHLUNG_H
|
||||
|
||||
#include "openmc/particle.h"
|
||||
|
||||
#include "xtensor/xtensor.hpp"
|
||||
|
||||
namespace openmc {
|
||||
|
||||
//==============================================================================
|
||||
// Bremsstrahlung classes
|
||||
//==============================================================================
|
||||
|
||||
class BremsstrahlungData {
|
||||
public:
|
||||
// Data
|
||||
xt::xtensor<double, 2> pdf; //!< Bremsstrahlung energy PDF
|
||||
xt::xtensor<double, 2> cdf; //!< Bremsstrahlung energy CDF
|
||||
xt::xtensor<double, 1> yield; //!< Photon yield
|
||||
};
|
||||
|
||||
class Bremsstrahlung {
|
||||
public:
|
||||
// Data
|
||||
BremsstrahlungData electron;
|
||||
BremsstrahlungData positron;
|
||||
};
|
||||
|
||||
//==============================================================================
|
||||
// Global variables
|
||||
//==============================================================================
|
||||
|
||||
namespace data {
|
||||
|
||||
extern xt::xtensor<double, 1> ttb_e_grid; //! energy T of incident electron in [eV]
|
||||
extern xt::xtensor<double, 1> ttb_k_grid; //! reduced energy W/T of emitted photon
|
||||
|
||||
} // namespace data
|
||||
|
||||
//==============================================================================
|
||||
// Global variables
|
||||
//==============================================================================
|
||||
|
||||
void thick_target_bremsstrahlung(Particle& p, double* E_lost);
|
||||
|
||||
} // namespace openmc
|
||||
|
||||
#endif // OPENMC_BREMSSTRAHLUNG_H
|
||||
|
|
@ -4,8 +4,8 @@
|
|||
#ifndef OPENMC_CONSTANTS_H
|
||||
#define OPENMC_CONSTANTS_H
|
||||
|
||||
#include <cmath>
|
||||
#include <array>
|
||||
#include <cmath>
|
||||
#include <limits>
|
||||
#include <vector>
|
||||
|
||||
|
|
@ -96,12 +96,12 @@ constexpr double N_AVOGADRO {0.6022140857}; // Avogadro's number in 10^24/
|
|||
constexpr double K_BOLTZMANN {8.6173303e-5}; // Boltzmann constant in eV/K
|
||||
|
||||
// Electron subshell labels
|
||||
constexpr char SUBSHELLS[][4] {
|
||||
"K ", "L1 ", "L2 ", "L3 ", "M1 ", "M2 ", "M3 ", "M4 ", "M5 ",
|
||||
"N1 ", "N2 ", "N3 ", "N4 ", "N5 ", "N6 ", "N7 ", "O1 ", "O2 ",
|
||||
"O3 ", "O4 ", "O5 ", "O6 ", "O7 ", "O8 ", "O9 ", "P1 ", "P2 ",
|
||||
"P3 ", "P4 ", "P5 ", "P6 ", "P7 ", "P8 ", "P9 ", "P10", "P11",
|
||||
"Q1 ", "Q2 ", "Q3 "
|
||||
constexpr std::array<const char*, 39> SUBSHELLS = {
|
||||
"K", "L1", "L2", "L3", "M1", "M2", "M3", "M4", "M5",
|
||||
"N1", "N2", "N3", "N4", "N5", "N6", "N7", "O1", "O2",
|
||||
"O3", "O4", "O5", "O6", "O7", "O8", "O9", "P1", "P2",
|
||||
"P3", "P4", "P5", "P6", "P7", "P8", "P9", "P10", "P11",
|
||||
"Q1", "Q2", "Q3"
|
||||
};
|
||||
|
||||
// Void material
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -1,11 +1,14 @@
|
|||
#ifndef OPENMC_MATERIAL_H
|
||||
#define OPENMC_MATERIAL_H
|
||||
|
||||
#include <memory> // for unique_ptr
|
||||
#include <unordered_map>
|
||||
#include <vector>
|
||||
|
||||
#include "pugixml.hpp"
|
||||
|
||||
#include "openmc/bremsstrahlung.h"
|
||||
#include "openmc/particle.h"
|
||||
|
||||
namespace openmc {
|
||||
|
||||
|
|
@ -38,15 +41,21 @@ public:
|
|||
//! A negative value indicates no default temperature was specified.
|
||||
double temperature_ {-1};
|
||||
|
||||
std::unique_ptr<Bremsstrahlung> ttb_;
|
||||
|
||||
Material() {};
|
||||
|
||||
explicit Material(pugi::xml_node material_node);
|
||||
|
||||
//! Initialize bremsstrahlung data
|
||||
void init_bremsstrahlung();
|
||||
};
|
||||
|
||||
//==============================================================================
|
||||
// Fortran compatibility
|
||||
//==============================================================================
|
||||
|
||||
extern "C" int* material_element(int i_material);
|
||||
extern "C" bool material_isotropic(int i_material, int i_nuc_mat);
|
||||
|
||||
} // namespace openmc
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
||||
|
|
@ -167,7 +167,7 @@ extern "C" double watt_spectrum(double a, double b);
|
|||
//!
|
||||
//! Samples from a Normal distribution with a given mean and standard deviation
|
||||
//! The PDF is defined as s(x) = (1/2*sigma*sqrt(2) * e-((mu-x)/2*sigma)^2
|
||||
//! Its sampled according to
|
||||
//! Its sampled according to
|
||||
//! http://www-pdg.lbl.gov/2009/reviews/rpp2009-rev-monte-carlo-techniques.pdf
|
||||
//! section 33.4.4
|
||||
//!
|
||||
|
|
@ -181,7 +181,7 @@ extern "C" double normal_variate(double mean, double std_dev);
|
|||
//==============================================================================
|
||||
//! Samples an energy from the Muir (Gaussian) energy-dependent distribution.
|
||||
//!
|
||||
//! This is another form of the Gaussian distribution but with more easily
|
||||
//! This is another form of the Gaussian distribution but with more easily
|
||||
//! modifiable parameters
|
||||
//! https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-05411-MS
|
||||
//!
|
||||
|
|
@ -193,8 +193,6 @@ extern "C" double normal_variate(double mean, double std_dev);
|
|||
|
||||
extern "C" double muir_spectrum(double e0, double m_rat, double kt);
|
||||
|
||||
|
||||
|
||||
//==============================================================================
|
||||
//! Doppler broadens the windowed multipole curvefit.
|
||||
//!
|
||||
|
|
@ -224,7 +222,7 @@ extern "C" void broaden_wmp_polynomials(double E, double dopp, int n, double fac
|
|||
//! value of x.
|
||||
//==============================================================================
|
||||
|
||||
extern "C" void spline_c(int n, const double x[], const double y[], double z[]);
|
||||
void spline(int n, const double x[], const double y[], double z[]);
|
||||
|
||||
//==============================================================================
|
||||
//! Determine the cubic spline interpolated y-value for a given x-value.
|
||||
|
|
@ -239,8 +237,8 @@ extern "C" void spline_c(int n, const double x[], const double y[], double z[]);
|
|||
//! \return Interpolated value
|
||||
//==============================================================================
|
||||
|
||||
extern "C" double spline_interpolate_c(int n, const double x[], const double y[],
|
||||
const double z[], double xint);
|
||||
double spline_interpolate(int n, const double x[], const double y[],
|
||||
const double z[], double xint);
|
||||
|
||||
//==============================================================================
|
||||
//! Evaluate the definite integral of the interpolating cubic spline between
|
||||
|
|
@ -257,9 +255,8 @@ extern "C" double spline_interpolate_c(int n, const double x[], const double y[]
|
|||
//! \return Integral
|
||||
//==============================================================================
|
||||
|
||||
extern "C" double spline_integrate_c(int n, const double x[], const double y[],
|
||||
const double z[], double xa, double xb);
|
||||
|
||||
double spline_integrate(int n, const double x[], const double y[],
|
||||
const double z[], double xa, double xb);
|
||||
|
||||
//! Evaluate the Faddeeva function
|
||||
//!
|
||||
|
|
|
|||
|
|
@ -1,9 +1,15 @@
|
|||
#ifndef OPENMC_PHOTON_H
|
||||
#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 {
|
||||
|
|
@ -12,26 +18,125 @@ namespace openmc {
|
|||
//! Photon interaction data for a single element
|
||||
//==============================================================================
|
||||
|
||||
class ElectronSubshell {
|
||||
public:
|
||||
// Constructors
|
||||
ElectronSubshell() { };
|
||||
|
||||
int index_subshell; //!< index in SUBSHELLS
|
||||
int threshold;
|
||||
double n_electrons;
|
||||
double binding_energy;
|
||||
xt::xtensor<double, 1> cross_section;
|
||||
|
||||
// Transition data
|
||||
int n_transitions;
|
||||
xt::xtensor<int, 2> transition_subshells;
|
||||
xt::xtensor<double, 1> transition_energy;
|
||||
xt::xtensor<double, 1> transition_probability;
|
||||
};
|
||||
|
||||
class PhotonInteraction {
|
||||
public:
|
||||
// Constructors
|
||||
PhotonInteraction(hid_t group);
|
||||
PhotonInteraction(hid_t group, int i_element);
|
||||
|
||||
// Methods
|
||||
void calculate_xs(double E) const;
|
||||
|
||||
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
|
||||
std::string name_; //!< Name of element, e.g. "Zr"
|
||||
int Z_; //!< Atomic number
|
||||
int i_element_; //!< Index in global elements vector
|
||||
|
||||
// Microscopic cross sections
|
||||
xt::xtensor<double, 1> energy_;
|
||||
xt::xtensor<double, 1> coherent_;
|
||||
xt::xtensor<double, 1> incoherent_;
|
||||
xt::xtensor<double, 1> photoelectric_total_;
|
||||
xt::xtensor<double, 1> pair_production_total_;
|
||||
xt::xtensor<double, 1> pair_production_electron_;
|
||||
xt::xtensor<double, 1> pair_production_nuclear_;
|
||||
|
||||
// Form factors
|
||||
Tabulated1D incoherent_form_factor_;
|
||||
Tabulated1D coherent_int_form_factor_;
|
||||
Tabulated1D coherent_anomalous_real_;
|
||||
Tabulated1D coherent_anomalous_imag_;
|
||||
|
||||
// Photoionization and atomic relaxation data
|
||||
std::unordered_map<int, int> shell_map_; //!< Given a shell designator, e.g. 3, this
|
||||
//!< dictionary gives an index in shells_
|
||||
std::vector<ElectronSubshell> shells_;
|
||||
|
||||
// Compton profile data
|
||||
xt::xtensor<double, 2> profile_pdf_;
|
||||
xt::xtensor<double, 2> profile_cdf_;
|
||||
xt::xtensor<double, 1> binding_energy_;
|
||||
xt::xtensor<double, 1> electron_pdf_;
|
||||
|
||||
// Stopping power data
|
||||
double I_; // mean excitation energy
|
||||
xt::xtensor<double, 1> stopping_power_collision_;
|
||||
xt::xtensor<double, 1> stopping_power_radiative_;
|
||||
|
||||
// Bremsstrahlung scaled DCS
|
||||
xt::xtensor<double, 2> dcs_;
|
||||
|
||||
private:
|
||||
void compton_doppler(double alpha, double mu, double* E_out, int* i_shell) const;
|
||||
};
|
||||
|
||||
//==============================================================================
|
||||
//! Cached microscopic photon cross sections for a particular element at the
|
||||
//! current energy
|
||||
//==============================================================================
|
||||
|
||||
struct ElementMicroXS {
|
||||
int index_grid; //!< index on element energy grid
|
||||
double last_E {0.0}; //!< last evaluated energy in [eV]
|
||||
double interp_factor; //!< interpolation factor on energy grid
|
||||
double total; //!< microscopic total photon xs
|
||||
double coherent; //!< microscopic coherent xs
|
||||
double incoherent; //!< microscopic incoherent xs
|
||||
double photoelectric; //!< microscopic photoelectric xs
|
||||
double pair_production; //!< microscopic pair production xs
|
||||
};
|
||||
|
||||
//==============================================================================
|
||||
// Non-member functions
|
||||
//==============================================================================
|
||||
|
||||
std::pair<double, double> klein_nishina(double alpha);
|
||||
|
||||
//==============================================================================
|
||||
// Global variables
|
||||
//==============================================================================
|
||||
|
||||
namespace data {
|
||||
|
||||
extern xt::xtensor<double, 1> compton_profile_pz; //! Compton profile momentum grid
|
||||
|
||||
//! Photon interaction data for each element
|
||||
extern std::vector<PhotonInteraction> elements;
|
||||
|
||||
} // namespace data
|
||||
|
||||
namespace simulation {
|
||||
extern ElementMicroXS* micro_photon_xs;
|
||||
#pragma omp threadprivate(micro_photon_xs)
|
||||
} // namespace simulation
|
||||
|
||||
} // namespace openmc
|
||||
|
||||
#endif // OPENMC_PHOTON_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);
|
||||
int sample_element(Particle* p);
|
||||
|
||||
Reaction* sample_fission(int i_nuclide, double E);
|
||||
|
||||
|
|
@ -82,8 +82,6 @@ void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p);
|
|||
|
||||
void sample_secondary_photons(Particle* p, int i_nuclide);
|
||||
|
||||
extern "C" void thick_target_bremsstrahlung(Particle* p, double* E_lost);
|
||||
|
||||
} // namespace openmc
|
||||
|
||||
#endif // OPENMC_PHYSICS_H
|
||||
|
|
|
|||
|
|
@ -68,7 +68,7 @@ extern "C" int32_t n_inactive; //!< number of inactive batches
|
|||
extern "C" int32_t gen_per_batch; //!< number of generations per batch
|
||||
extern "C" int64_t n_particles; //!< number of particles per generation
|
||||
|
||||
extern "C" int electron_treatment; //!< how to treat secondary electrons
|
||||
extern int electron_treatment; //!< how to treat secondary electrons
|
||||
extern "C" std::array<double, 4> energy_cutoff; //!< Energy cutoff in [eV] for each particle type
|
||||
extern "C" int legendre_to_tabular_points; //!< number of points to convert Legendres
|
||||
extern "C" int max_order; //!< Maximum Legendre order for multigroup data
|
||||
|
|
|
|||
129
src/bremsstrahlung.cpp
Normal file
129
src/bremsstrahlung.cpp
Normal file
|
|
@ -0,0 +1,129 @@
|
|||
#include "openmc/bremsstrahlung.h"
|
||||
|
||||
#include "openmc/constants.h"
|
||||
#include "openmc/material.h"
|
||||
#include "openmc/random_lcg.h"
|
||||
#include "openmc/search.h"
|
||||
#include "openmc/settings.h"
|
||||
|
||||
#include "xtensor/xmath.hpp"
|
||||
|
||||
namespace openmc {
|
||||
|
||||
//==============================================================================
|
||||
// Global variables
|
||||
//==============================================================================
|
||||
|
||||
namespace data {
|
||||
|
||||
xt::xtensor<double, 1> ttb_e_grid;
|
||||
xt::xtensor<double, 1> ttb_k_grid;
|
||||
std::vector<Bremsstrahlung> ttb;
|
||||
|
||||
} // namespace data
|
||||
|
||||
//==============================================================================
|
||||
// Non-member functions
|
||||
//==============================================================================
|
||||
|
||||
void thick_target_bremsstrahlung(Particle& p, double* E_lost)
|
||||
{
|
||||
if (p.material == MATERIAL_VOID) return;
|
||||
|
||||
// TODO: off-by-one
|
||||
int photon = static_cast<int>(ParticleType::photon) - 1;
|
||||
if (p.E < settings::energy_cutoff[photon]) return;
|
||||
|
||||
// Get bremsstrahlung data for this material and particle type
|
||||
// TODO: off-by-one
|
||||
BremsstrahlungData* mat;
|
||||
if (p.type == static_cast<int>(ParticleType::positron)) {
|
||||
mat = &model::materials[p.material -1]->ttb_->positron;
|
||||
} else {
|
||||
mat = &model::materials[p.material -1]->ttb_->electron;
|
||||
}
|
||||
|
||||
double e = std::log(p.E);
|
||||
auto n_e = data::ttb_e_grid.size();
|
||||
|
||||
// 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;
|
||||
int i_e;
|
||||
if (prn() <= f || j == 0) {
|
||||
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, i_e - 1);
|
||||
double p_r = mat->pdf(i_e, i_e);
|
||||
double c_l = mat->cdf(i_e, i_e - 1);
|
||||
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 {
|
||||
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, 0), &mat->cdf(i_e, 0) + 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_e, i_w);
|
||||
double p_r = mat->pdf(i_e, i_w + 1);
|
||||
double c_l = mat->cdf(i_e, i_w);
|
||||
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
|
||||
//==============================================================================
|
||||
|
||||
extern "C" void set_log_ttb_e_grid()
|
||||
{
|
||||
data::ttb_e_grid = xt::log(data::ttb_e_grid);
|
||||
}
|
||||
|
||||
} // namespace openmc
|
||||
|
|
@ -68,13 +68,8 @@ module constants
|
|||
|
||||
real(8), parameter :: &
|
||||
PI = 3.1415926535898_8, & ! pi
|
||||
SQRT_PI = 1.7724538509055_8, & ! square root of pi
|
||||
MASS_NEUTRON = 1.00866491588_8, & ! mass of a neutron in amu
|
||||
MASS_NEUTRON_EV = 939.5654133e6_8, & ! mass of a neutron in eV/c^2
|
||||
MASS_PROTON = 1.007276466879_8, & ! mass of a proton in amu
|
||||
MASS_ELECTRON_EV = 0.5109989461e6_8, & ! electron mass energy equivalent in eV/c^2
|
||||
FINE_STRUCTURE = 137.035999139_8, & ! inverse fine structure constant
|
||||
PLANCK_C = 1.2398419739062977e4_8,& ! Planck's constant times c in eV-Angstroms
|
||||
AMU = 1.660539040e-27_8, & ! 1 amu in kg
|
||||
C_LIGHT = 2.99792458e8_8, & ! speed of light in m/s
|
||||
N_AVOGADRO = 0.6022140857_8, & ! Avogadro's number in 10^24/mol
|
||||
|
|
@ -83,10 +78,7 @@ module constants
|
|||
ZERO = 0.0_8, &
|
||||
HALF = 0.5_8, &
|
||||
ONE = 1.0_8, &
|
||||
TWO = 2.0_8, &
|
||||
THREE = 3.0_8, &
|
||||
FOUR = 4.0_8
|
||||
complex(8), parameter :: ONEI = (ZERO, ONE)
|
||||
TWO = 2.0_8
|
||||
|
||||
! Electron subshell labels
|
||||
character(3), parameter :: SUBSHELLS(39) = [ &
|
||||
|
|
@ -185,12 +177,6 @@ module constants
|
|||
real(8), parameter :: &
|
||||
MACROSCOPIC_AWR = -TWO
|
||||
|
||||
! Fission neutron emission (nu) type
|
||||
integer, parameter :: &
|
||||
NU_NONE = 0, & ! No nu values (non-fissionable)
|
||||
NU_POLYNOMIAL = 1, & ! Nu values given by polynomial
|
||||
NU_TABULAR = 2 ! Nu values given by tabular distribution
|
||||
|
||||
! Secondary particle emission type
|
||||
integer, parameter :: &
|
||||
EMISSION_PROMPT = 1, & ! Prompt emission of secondary particle
|
||||
|
|
@ -291,11 +277,6 @@ module constants
|
|||
FILTER_ZERNIKE_RADIAL = 21, &
|
||||
FILTER_PARTICLE = 22
|
||||
|
||||
|
||||
! Mesh types
|
||||
integer, parameter :: &
|
||||
MESH_REGULAR = 1
|
||||
|
||||
! Tally surface current directions
|
||||
integer, parameter :: &
|
||||
OUT_LEFT = 1, & ! x min
|
||||
|
|
@ -382,11 +363,6 @@ module constants
|
|||
MODE_PARTICLE = 4, & ! Particle restart mode
|
||||
MODE_VOLUME = 5 ! Volume calculation mode
|
||||
|
||||
! Electron treatments
|
||||
integer, parameter :: &
|
||||
ELECTRON_LED = 1, & ! Local Energy Deposition
|
||||
ELECTRON_TTB = 2 ! Thick Target Bremsstrahlung
|
||||
|
||||
!=============================================================================
|
||||
! DELAYED NEUTRON PRECURSOR CONSTANTS
|
||||
|
||||
|
|
|
|||
|
|
@ -2034,11 +2034,17 @@ contains
|
|||
type(SetChar) :: already_read
|
||||
type(SetChar) :: element_already_read
|
||||
|
||||
interface
|
||||
subroutine photon_from_hdf5(group) bind(C)
|
||||
import HID_T
|
||||
integer(HID_T), value :: group
|
||||
end subroutine
|
||||
|
||||
subroutine read_ce_cross_sections_c() bind(C)
|
||||
end subroutine
|
||||
end interface
|
||||
|
||||
allocate(nuclides(n_nuclides))
|
||||
allocate(elements(n_elements))
|
||||
if (photon_transport .and. electron_treatment == ELECTRON_TTB) then
|
||||
allocate(ttb(n_materials))
|
||||
end if
|
||||
|
||||
! Read cross sections
|
||||
do i = 1, size(materials)
|
||||
|
|
@ -2094,22 +2100,11 @@ contains
|
|||
|
||||
! Read element data from HDF5
|
||||
group_id = open_group(file_id, element)
|
||||
call elements(i_element) % from_hdf5(group_id)
|
||||
|
||||
call photon_from_hdf5(group_id)
|
||||
call close_group(group_id)
|
||||
call file_close(file_id)
|
||||
|
||||
! Determine if minimum/maximum energy for this element is
|
||||
! greater/less than the previous
|
||||
if (size(elements(i_element) % energy) >= 1) then
|
||||
energy_min(PHOTON) = max(energy_min(PHOTON), &
|
||||
exp(elements(i_element) % energy(1)))
|
||||
energy_max(PHOTON) = min(energy_max(PHOTON), &
|
||||
exp(elements(i_element) % energy(size(elements(i_element) &
|
||||
% energy))))
|
||||
call set_particle_energy_bounds(PHOTON, energy_min(PHOTON), &
|
||||
energy_max(PHOTON))
|
||||
end if
|
||||
|
||||
! Add element to set
|
||||
call element_already_read % add(element)
|
||||
end if
|
||||
|
|
@ -2124,38 +2119,9 @@ contains
|
|||
call materials(i) % set_fissionable(.true.)
|
||||
end if
|
||||
end do
|
||||
|
||||
! Generate material bremsstrahlung data for electrons and positrons
|
||||
if (photon_transport .and. electron_treatment == ELECTRON_TTB) then
|
||||
call bremsstrahlung_init(ttb(i) % electron, i, ELECTRON)
|
||||
call bremsstrahlung_init(ttb(i) % positron, i, POSITRON)
|
||||
end if
|
||||
end do
|
||||
|
||||
if (photon_transport .and. electron_treatment == ELECTRON_TTB) then
|
||||
! Deallocate element bremsstrahlung DCS and stopping power data since
|
||||
! only the material bremsstrahlung data is needed
|
||||
do i = 1, size(elements)
|
||||
if (allocated(elements(i) % stopping_power_collision)) &
|
||||
deallocate(elements(i) % stopping_power_collision)
|
||||
if (allocated(elements(i) % stopping_power_radiative)) &
|
||||
deallocate(elements(i) % stopping_power_radiative)
|
||||
if (allocated(elements(i) % dcs)) deallocate(elements(i) % dcs)
|
||||
if (allocated(ttb_k_grid)) deallocate(ttb_k_grid)
|
||||
end do
|
||||
|
||||
! Determine if minimum/maximum energy for bremsstrahlung is greater/less
|
||||
! than the current minimum/maximum
|
||||
if (size(ttb_e_grid) >= 1) then
|
||||
energy_min(PHOTON) = max(energy_min(PHOTON), ttb_e_grid(1))
|
||||
energy_max(PHOTON) = min(energy_max(PHOTON), ttb_e_grid(size(ttb_e_grid)))
|
||||
call set_particle_energy_bounds(PHOTON, energy_min(PHOTON), &
|
||||
energy_max(PHOTON))
|
||||
end if
|
||||
|
||||
! Take logarithm of energies since they are log-log interpolated
|
||||
ttb_e_grid = log(ttb_e_grid)
|
||||
end if
|
||||
call read_ce_cross_sections_c()
|
||||
|
||||
! Set up logarithmic grid for nuclides
|
||||
do i = 1, size(nuclides)
|
||||
|
|
|
|||
234
src/material.cpp
234
src/material.cpp
|
|
@ -1,11 +1,21 @@
|
|||
#include "openmc/material.h"
|
||||
|
||||
#include <algorithm> // for min, max
|
||||
#include <cmath>
|
||||
#include <string>
|
||||
#include <sstream>
|
||||
|
||||
#include "openmc/error.h"
|
||||
#include "openmc/xml_interface.h"
|
||||
#include "xtensor/xbuilder.hpp"
|
||||
#include "xtensor/xoperation.hpp"
|
||||
#include "xtensor/xview.hpp"
|
||||
|
||||
#include "openmc/error.h"
|
||||
#include "openmc/math_functions.h"
|
||||
#include "openmc/nuclide.h"
|
||||
#include "openmc/photon.h"
|
||||
#include "openmc/search.h"
|
||||
#include "openmc/settings.h"
|
||||
#include "openmc/xml_interface.h"
|
||||
|
||||
namespace openmc {
|
||||
|
||||
|
|
@ -41,6 +51,196 @@ Material::Material(pugi::xml_node node)
|
|||
}
|
||||
}
|
||||
|
||||
void Material::init_bremsstrahlung()
|
||||
{
|
||||
// Create new object
|
||||
ttb_ = std::make_unique<Bremsstrahlung>();
|
||||
|
||||
// Get the size of the energy grids
|
||||
auto n_k = data::ttb_k_grid.size();
|
||||
auto n_e = data::ttb_e_grid.size();
|
||||
|
||||
// Get pointers to nuclides, elements, densities
|
||||
int32_t index;
|
||||
openmc_get_material_index(id_, &index);
|
||||
int* nuclide_;
|
||||
double* atom_density_;
|
||||
int n;
|
||||
openmc_material_get_densities(index, &nuclide_, &atom_density_, &n);
|
||||
int* element_ = material_element(index);
|
||||
|
||||
for (int particle = 0; particle < 2; ++particle) {
|
||||
// Loop over logic twice, once for electron, once for positron
|
||||
BremsstrahlungData* ttb = (particle == 0) ? &ttb_->electron : &ttb_->positron;
|
||||
bool positron = (particle == 1);
|
||||
|
||||
// Allocate arrays for TTB data
|
||||
ttb->pdf = xt::zeros<double>({n_e, n_e});
|
||||
ttb->cdf = xt::zeros<double>({n_e, n_e});
|
||||
ttb->yield = xt::empty<double>({n_e});
|
||||
|
||||
// Allocate temporary arrays
|
||||
xt::xtensor<double, 1> stopping_power_collision({n_e}, 0.0);
|
||||
xt::xtensor<double, 1> stopping_power_radiative({n_e}, 0.0);
|
||||
xt::xtensor<double, 2> dcs({n_e, n_k}, 0.0);
|
||||
|
||||
double Z_eq_sq = 0.0;
|
||||
double sum_density = 0.0;
|
||||
|
||||
// Calculate the molecular DCS and the molecular total stopping power using
|
||||
// Bragg's additivity rule.
|
||||
// TODO: The collision stopping power cannot be accurately calculated using
|
||||
// Bragg's additivity rule since the mean excitation energies and the
|
||||
// density effect corrections cannot simply be summed together. Bragg's
|
||||
// additivity rule fails especially when a higher-density compound is
|
||||
// composed of elements that are in lower-density form at normal temperature
|
||||
// and pressure (at which the NIST stopping powers are given). It will be
|
||||
// used to approximate the collision stopping powers for now, but should be
|
||||
// fixed in the future.
|
||||
for (int i = 0; i < n; ++i) {
|
||||
// Get pointer to current element
|
||||
// TODO: off-by-one
|
||||
const auto& elm = data::elements[element_[i] - 1];
|
||||
// TODO: off-by-one
|
||||
double awr = data::nuclides[nuclide_[i] - 1]->awr_;
|
||||
|
||||
// Get atomic density and mass density of nuclide given atom/weight percent
|
||||
double atom_density = (atom_density_[0] > 0.0) ?
|
||||
atom_density_[i] : -atom_density_[i] / awr;
|
||||
double mass_density = atom_density * awr;
|
||||
|
||||
// Calculate the "equivalent" atomic number Zeq of the material
|
||||
Z_eq_sq += atom_density * elm.Z_ * elm.Z_;
|
||||
sum_density += atom_density;
|
||||
|
||||
// Accumulate material DCS
|
||||
dcs += (atom_density * elm.Z_ * elm.Z_) * elm.dcs_;
|
||||
|
||||
// Accumulate material collision stopping power
|
||||
stopping_power_collision += (mass_density * MASS_NEUTRON / N_AVOGADRO)
|
||||
* elm.stopping_power_collision_;
|
||||
|
||||
// Accumulate material radiative stopping power
|
||||
stopping_power_radiative += (mass_density * MASS_NEUTRON / N_AVOGADRO)
|
||||
* elm.stopping_power_radiative_;
|
||||
}
|
||||
Z_eq_sq /= sum_density;
|
||||
|
||||
// Calculate the positron DCS and radiative stopping power. These are
|
||||
// obtained by multiplying the electron DCS and radiative stopping powers by
|
||||
// a factor r, which is a numerical approximation of the ratio of the
|
||||
// radiative stopping powers for positrons and electrons. Source: F. Salvat,
|
||||
// J. M. Fernández-Varea, and J. Sempau, "PENELOPE-2011: A Code System for
|
||||
// Monte Carlo Simulation of Electron and Photon Transport," OECD-NEA,
|
||||
// Issy-les-Moulineaux, France (2011).
|
||||
if (positron) {
|
||||
for (int i = 0; i < n_e; ++i) {
|
||||
double t = std::log(1.0 + 1.0e6*data::ttb_e_grid(i)/(Z_eq_sq*MASS_ELECTRON_EV));
|
||||
double r = 1.0 - std::exp(-1.2359e-1*t + 6.1274e-2*std::pow(t, 2)
|
||||
- 3.1516e-2*std::pow(t, 3) + 7.7446e-3*std::pow(t, 4)
|
||||
- 1.0595e-3*std::pow(t, 5) + 7.0568e-5*std::pow(t, 6)
|
||||
- 1.808e-6*std::pow(t, 7));
|
||||
stopping_power_radiative(i) *= r;
|
||||
auto dcs_i = xt::view(dcs, i, xt::all());
|
||||
dcs_i *= r;
|
||||
}
|
||||
}
|
||||
|
||||
// Total material stopping power
|
||||
xt::xtensor<double, 1> stopping_power = stopping_power_collision +
|
||||
stopping_power_radiative;
|
||||
|
||||
// Loop over photon energies
|
||||
xt::xtensor<double, 1> f({n_e}, 0.0);
|
||||
xt::xtensor<double, 1> z({n_e}, 0.0);
|
||||
for (int i = 0; i < n_e - 1; ++i) {
|
||||
double w = data::ttb_e_grid(i);
|
||||
|
||||
// Loop over incident particle energies
|
||||
for (int j = i; j < n_e; ++j) {
|
||||
double e = data::ttb_e_grid(j);
|
||||
|
||||
// Reduced photon energy
|
||||
double k = w / e;
|
||||
|
||||
// Find the lower bounding index of the reduced photon energy
|
||||
int i_k = lower_bound_index(data::ttb_k_grid.cbegin(),
|
||||
data::ttb_k_grid.cend(), k);
|
||||
|
||||
// Get the interpolation bounds
|
||||
double k_l = data::ttb_k_grid(i_k);
|
||||
double k_r = data::ttb_k_grid(i_k + 1);
|
||||
double x_l = dcs(j, i_k);
|
||||
double x_r = dcs(j, i_k + 1);
|
||||
|
||||
// Find the value of the DCS using linear interpolation in reduced
|
||||
// photon energy k
|
||||
double x = x_l + (k - k_l)*(x_r - x_l)/(k_r - k_l);
|
||||
|
||||
// Ratio of the velocity of the charged particle to the speed of light
|
||||
double beta = std::sqrt(e*(e + 2.0*MASS_ELECTRON_EV)) /
|
||||
(e + MASS_ELECTRON_EV);
|
||||
|
||||
// Compute the integrand of the PDF
|
||||
f(j) = x / (beta*beta * stopping_power(j) * w);
|
||||
}
|
||||
|
||||
// Number of points to integrate
|
||||
int n = n_e - i;
|
||||
|
||||
// Integrate the PDF using cubic spline integration over the incident
|
||||
// particle energy
|
||||
if (n > 2) {
|
||||
spline(n, &data::ttb_e_grid(i), &f(i), &z(i));
|
||||
|
||||
double c = 0.0;
|
||||
for (int j = i; j < n_e - 1; ++j) {
|
||||
c += spline_integrate(n, &data::ttb_e_grid(i), &f(i), &z(i),
|
||||
data::ttb_e_grid(j), data::ttb_e_grid(j+1));
|
||||
|
||||
ttb->pdf(j+1,i) = c;
|
||||
}
|
||||
|
||||
// Integrate the last two points using trapezoidal rule in log-log space
|
||||
} else {
|
||||
double e_l = std::log(data::ttb_e_grid(i));
|
||||
double e_r = std::log(data::ttb_e_grid(i+1));
|
||||
double x_l = std::log(f(i));
|
||||
double x_r = std::log(f(i+1));
|
||||
|
||||
ttb->pdf(i+1,i) = 0.5*(e_r - e_l)*(std::exp(e_l + x_l) + std::exp(e_r + x_r));
|
||||
}
|
||||
}
|
||||
|
||||
// Loop over incident particle energies
|
||||
for (int j = 1; j < n_e; ++j) {
|
||||
// Set last element of PDF to small non-zero value to enable log-log
|
||||
// interpolation
|
||||
ttb->pdf(j,j) = std::exp(-500.0);
|
||||
|
||||
// Loop over photon energies
|
||||
double c = 0.0;
|
||||
for (int i = 0; i < j; ++i) {
|
||||
// Integrate the CDF from the PDF using the trapezoidal rule in log-log
|
||||
// space
|
||||
double w_l = std::log(data::ttb_e_grid(i));
|
||||
double w_r = std::log(data::ttb_e_grid(i+1));
|
||||
double x_l = std::log(ttb->pdf(j,i));
|
||||
double x_r = std::log(ttb->pdf(j,i+1));
|
||||
|
||||
c += 0.5*(w_r - w_l)*(std::exp(w_l + x_l) + std::exp(w_r + x_r));
|
||||
ttb->cdf(j,i+1) = c;
|
||||
}
|
||||
|
||||
// Set photon number yield
|
||||
ttb->yield(j) = c;
|
||||
}
|
||||
|
||||
// Use logarithm of number yield since it is log-log interpolated
|
||||
ttb->yield = xt::where(ttb->yield > 0.0, xt::log(ttb->yield), -500.0);
|
||||
}
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// Non-method functions
|
||||
//==============================================================================
|
||||
|
|
@ -68,6 +268,31 @@ read_materials(pugi::xml_node* node)
|
|||
}
|
||||
}
|
||||
|
||||
extern "C" void read_ce_cross_sections_c()
|
||||
{
|
||||
for (auto& mat : model::materials) {
|
||||
// Generate material bremsstrahlung data for electrons and positrons
|
||||
if (settings::photon_transport && settings::electron_treatment == ELECTRON_TTB) {
|
||||
mat->init_bremsstrahlung();
|
||||
}
|
||||
}
|
||||
|
||||
if (settings::photon_transport && settings::electron_treatment == ELECTRON_TTB) {
|
||||
// Determine if minimum/maximum energy for bremsstrahlung is greater/less
|
||||
// than the current minimum/maximum
|
||||
if (data::ttb_e_grid.size() >= 1) {
|
||||
// TODO: off-by-one
|
||||
int photon = static_cast<int>(ParticleType::photon) - 1;
|
||||
int n_e = data::ttb_e_grid.size();
|
||||
data::energy_min[photon] = std::max(data::energy_min[photon], data::ttb_e_grid(1));
|
||||
data::energy_max[photon] = std::min(data::energy_max[photon], data::ttb_e_grid(n_e - 1));
|
||||
}
|
||||
|
||||
// Take logarithm of energies since they are log-log interpolated
|
||||
data::ttb_e_grid = xt::log(data::ttb_e_grid);
|
||||
}
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// C API
|
||||
//==============================================================================
|
||||
|
|
@ -133,6 +358,11 @@ extern "C" {
|
|||
mat->fissionable = fissionable;
|
||||
}
|
||||
|
||||
void material_init_bremsstrahlung(Material* mat)
|
||||
{
|
||||
mat->init_bremsstrahlung();
|
||||
}
|
||||
|
||||
void extend_materials_c(int32_t n)
|
||||
{
|
||||
model::materials.reserve(model::materials.size() + n);
|
||||
|
|
|
|||
|
|
@ -5,7 +5,6 @@ module material_header
|
|||
use constants
|
||||
use dict_header, only: DictIntInt
|
||||
use error
|
||||
use math, only: spline, spline_integrate
|
||||
use nuclide_header
|
||||
use particle_header, only: Particle
|
||||
use photon_header
|
||||
|
|
@ -17,7 +16,6 @@ module material_header
|
|||
implicit none
|
||||
|
||||
private
|
||||
public :: bremsstrahlung_init
|
||||
public :: free_memory_material
|
||||
public :: openmc_extend_materials
|
||||
public :: openmc_get_material_index
|
||||
|
|
@ -487,6 +485,14 @@ contains
|
|||
integer :: i_element ! index into elements array
|
||||
real(8) :: atom_density ! atom density of a nuclide
|
||||
|
||||
interface
|
||||
subroutine photon_calculate_xs(i_element, E) bind(C)
|
||||
import C_INT, C_DOUBLE
|
||||
integer(C_INT), value :: i_element
|
||||
real(C_DOUBLE), value :: E
|
||||
end subroutine
|
||||
end interface
|
||||
|
||||
material_xs % coherent = ZERO
|
||||
material_xs % incoherent = ZERO
|
||||
material_xs % photoelectric = ZERO
|
||||
|
|
@ -502,8 +508,7 @@ contains
|
|||
|
||||
! Calculate microscopic cross section for this nuclide
|
||||
if (p % E /= micro_photon_xs(i_element) % last_E) then
|
||||
call elements(i_element) % calculate_xs(&
|
||||
p % E, micro_photon_xs(i_element))
|
||||
call photon_calculate_xs(i_element, p % E)
|
||||
end if
|
||||
|
||||
! ========================================================================
|
||||
|
|
@ -830,217 +835,6 @@ contains
|
|||
|
||||
end function openmc_material_set_densities
|
||||
|
||||
subroutine bremsstrahlung_init(this, i_material, particle)
|
||||
class(BremsstrahlungData), intent(inout) :: this
|
||||
integer, intent(in) :: i_material
|
||||
integer, intent(in) :: particle
|
||||
|
||||
integer :: i, j
|
||||
integer :: i_k
|
||||
integer :: n, n_e, n_k
|
||||
real(8) :: c
|
||||
real(8) :: k, k_l, k_r
|
||||
real(8) :: e, e_l, e_r
|
||||
real(8) :: w, w_l, w_r
|
||||
real(8) :: x, x_l, x_r
|
||||
real(8) :: t
|
||||
real(8) :: r
|
||||
real(8) :: awr
|
||||
real(8) :: beta
|
||||
real(8) :: Z_eq_sq
|
||||
real(8) :: atom_density
|
||||
real(8) :: mass_density
|
||||
real(8) :: sum_density
|
||||
real(8), allocatable :: stopping_power_collision(:)
|
||||
real(8), allocatable :: stopping_power_radiative(:)
|
||||
real(8), allocatable :: stopping_power(:)
|
||||
real(8), allocatable :: dcs(:,:)
|
||||
real(8), allocatable :: f(:)
|
||||
real(8), allocatable :: z(:)
|
||||
logical :: positron_
|
||||
type(Material), pointer :: mat
|
||||
type(PhotonInteraction), pointer :: elm
|
||||
|
||||
! Get pointer to this material
|
||||
mat => materials(i_material)
|
||||
|
||||
! Determine whether we are generating electron or positron data
|
||||
positron_ = (particle == POSITRON)
|
||||
|
||||
! Get the size of the energy grids
|
||||
n_k = size(ttb_k_grid)
|
||||
n_e = size(ttb_e_grid)
|
||||
|
||||
! Allocate arrays for TTB data
|
||||
allocate(this % pdf(n_e, n_e), source=ZERO)
|
||||
allocate(this % cdf(n_e, n_e), source=ZERO)
|
||||
allocate(this % yield(n_e))
|
||||
|
||||
! Allocate temporary arrays
|
||||
allocate(stopping_power_collision(n_e), source=ZERO)
|
||||
allocate(stopping_power_radiative(n_e), source=ZERO)
|
||||
allocate(stopping_power(n_e))
|
||||
allocate(dcs(n_k, n_e), source=ZERO)
|
||||
allocate(f(n_e))
|
||||
allocate(z(n_e))
|
||||
|
||||
Z_eq_sq = ZERO
|
||||
sum_density = ZERO
|
||||
|
||||
! Calculate the molecular DCS and the molecular total stopping power using
|
||||
! Bragg's additivity rule.
|
||||
! TODO: The collision stopping power cannot be accurately calculated using
|
||||
! Bragg's additivity rule since the mean excitation energies and the
|
||||
! density effect corrections cannot simply be summed together. Bragg's
|
||||
! additivity rule fails especially when a higher-density compound is
|
||||
! composed of elements that are in lower-density form at normal temperature
|
||||
! and pressure (at which the NIST stopping powers are given). It will be
|
||||
! used to approximate the collision stopping powers for now, but should be
|
||||
! fixed in the future.
|
||||
do i = 1, mat % n_nuclides
|
||||
! Get pointer to current element
|
||||
elm => elements(mat % element(i))
|
||||
|
||||
awr = nuclides(mat % nuclide(i)) % awr
|
||||
|
||||
! Get atomic density and mass density of nuclide given atom percent
|
||||
if (mat % atom_density(1) > ZERO) then
|
||||
atom_density = mat % atom_density(i)
|
||||
mass_density = mat % atom_density(i) * awr
|
||||
! Given weight percent
|
||||
else
|
||||
atom_density = -mat % atom_density(i) / awr
|
||||
mass_density = -mat % atom_density(i)
|
||||
end if
|
||||
|
||||
! Calculate the "equivalent" atomic number Zeq of the material
|
||||
Z_eq_sq = Z_eq_sq + atom_density * elm % Z**2
|
||||
sum_density = sum_density + atom_density
|
||||
|
||||
! Accumulate material DCS
|
||||
dcs = dcs + atom_density * elm % Z**2 * elm % dcs
|
||||
|
||||
! Accumulate material collision stopping power
|
||||
stopping_power_collision = stopping_power_collision + mass_density &
|
||||
* MASS_NEUTRON / N_AVOGADRO * elm % stopping_power_collision
|
||||
|
||||
! Accumulate material radiative stopping power
|
||||
stopping_power_radiative = stopping_power_radiative + mass_density &
|
||||
* MASS_NEUTRON / N_AVOGADRO * elm % stopping_power_radiative
|
||||
end do
|
||||
Z_eq_sq = Z_eq_sq / sum_density
|
||||
|
||||
! Calculate the positron DCS and radiative stopping power. These are
|
||||
! obtained by multiplying the electron DCS and radiative stopping powers by
|
||||
! a factor r, which is a numerical approximation of the ratio of the
|
||||
! radiative stopping powers for positrons and electrons. Source: F. Salvat,
|
||||
! J. M. Fernández-Varea, and J. Sempau, "PENELOPE-2011: A Code System for
|
||||
! Monte Carlo Simulation of Electron and Photon Transport," OECD-NEA,
|
||||
! Issy-les-Moulineaux, France (2011).
|
||||
if (positron_) then
|
||||
do i = 1, n_e
|
||||
t = log(ONE + 1.0e6_8*ttb_e_grid(i)/(Z_eq_sq*MASS_ELECTRON_EV))
|
||||
r = ONE - exp(-1.2359e-1_8*t + 6.1274e-2_8*t**2 - 3.1516e-2_8*t**3 + &
|
||||
7.7446e-3_8*t**4 - 1.0595e-3_8*t**5 + 7.0568e-5_8*t**6 - &
|
||||
1.808e-6_8*t**7)
|
||||
stopping_power_radiative(i) = r*stopping_power_radiative(i)
|
||||
dcs(:,i) = r*dcs(:,i)
|
||||
end do
|
||||
end if
|
||||
|
||||
! Total material stopping power
|
||||
stopping_power = stopping_power_collision + stopping_power_radiative
|
||||
|
||||
! Loop over photon energies
|
||||
do i = 1, n_e - 1
|
||||
w = ttb_e_grid(i)
|
||||
|
||||
! Loop over incident particle energies
|
||||
do j = i, n_e
|
||||
e = ttb_e_grid(j)
|
||||
|
||||
! Reduced photon energy
|
||||
k = w / e
|
||||
|
||||
! Find the lower bounding index of the reduced photon energy
|
||||
i_k = binary_search(ttb_k_grid, n_k, k)
|
||||
|
||||
! Get the interpolation bounds
|
||||
k_l = ttb_k_grid(i_k)
|
||||
k_r = ttb_k_grid(i_k+1)
|
||||
x_l = dcs(i_k, j)
|
||||
x_r = dcs(i_k+1, j)
|
||||
|
||||
! Find the value of the DCS using linear interpolation in reduced
|
||||
! photon energy k
|
||||
x = x_l + (k - k_l) * (x_r - x_l) / (k_r - k_l)
|
||||
|
||||
! Ratio of the velocity of the charged particle to the speed of light
|
||||
beta = sqrt(e*(e + TWO*MASS_ELECTRON_EV)) / (e + MASS_ELECTRON_EV)
|
||||
|
||||
! Compute the integrand of the PDF
|
||||
f(j) = x / (beta**2 * stopping_power(j) * w)
|
||||
end do
|
||||
|
||||
! Number of points to integrate
|
||||
n = n_e - i + 1
|
||||
|
||||
! Integrate the PDF using cubic spline integration over the incident
|
||||
! particle energy
|
||||
if (n > 2) then
|
||||
call spline(n, ttb_e_grid(i:), f(i:), z(i:))
|
||||
|
||||
c = ZERO
|
||||
do j = i, n_e - 1
|
||||
c = c + spline_integrate(n, ttb_e_grid(i:), f(i:), z(i:), &
|
||||
ttb_e_grid(j), ttb_e_grid(j+1))
|
||||
this % pdf(i,j+1) = c
|
||||
end do
|
||||
|
||||
! Integrate the last two points using trapezoidal rule in log-log space
|
||||
else
|
||||
e_l = log(ttb_e_grid(i))
|
||||
e_r = log(ttb_e_grid(i+1))
|
||||
x_l = log(f(i))
|
||||
x_r = log(f(i+1))
|
||||
|
||||
this % pdf(i,i+1) = HALF * (e_r - e_l) * (exp(e_l + x_l) + exp(e_r + x_r))
|
||||
end if
|
||||
end do
|
||||
|
||||
! Loop over incident particle energies
|
||||
do j = 2, n_e
|
||||
! Set last element of PDF to small non-zero value to enable log-log
|
||||
! interpolation
|
||||
this % pdf(j,j) = exp(-500.0_8)
|
||||
|
||||
! Loop over photon energies
|
||||
c = ZERO
|
||||
do i = 1, j - 1
|
||||
! Integrate the CDF from the PDF using the trapezoidal rule in log-log
|
||||
! space
|
||||
w_l = log(ttb_e_grid(i))
|
||||
w_r = log(ttb_e_grid(i+1))
|
||||
x_l = log(this % pdf(i,j))
|
||||
x_r = log(this % pdf(i+1,j))
|
||||
|
||||
c = c + HALF * (w_r - w_l) * (exp(w_l + x_l) + exp(w_r + x_r))
|
||||
this % cdf(i+1,j) = c
|
||||
end do
|
||||
|
||||
! Set photon number yield
|
||||
this % yield(j) = c
|
||||
end do
|
||||
|
||||
! Use logarithm of number yield since it is log-log interpolated
|
||||
where (this % yield > ZERO)
|
||||
this % yield = log(this % yield)
|
||||
elsewhere
|
||||
this % yield = -500.0_8
|
||||
end where
|
||||
|
||||
end subroutine bremsstrahlung_init
|
||||
|
||||
!===============================================================================
|
||||
! Fortran compatibility
|
||||
!===============================================================================
|
||||
|
|
@ -1058,4 +852,10 @@ contains
|
|||
end associate
|
||||
end function
|
||||
|
||||
function material_element(i_material) result(ptr) bind(C)
|
||||
integer(C_INT), value :: i_material
|
||||
type(C_PTR) :: ptr
|
||||
ptr = C_LOC(materials(i_material) % element(1))
|
||||
end function
|
||||
|
||||
end module material_header
|
||||
|
|
|
|||
40
src/math.F90
40
src/math.F90
|
|
@ -2,18 +2,12 @@ module math
|
|||
|
||||
use, intrinsic :: ISO_C_BINDING
|
||||
|
||||
use constants
|
||||
use random_lcg, only: prn
|
||||
|
||||
implicit none
|
||||
private
|
||||
public :: t_percentile
|
||||
public :: calc_pn
|
||||
public :: calc_rn
|
||||
public :: rotate_angle
|
||||
public :: spline
|
||||
public :: spline_interpolate
|
||||
public :: spline_integrate
|
||||
|
||||
interface
|
||||
|
||||
|
|
@ -49,40 +43,6 @@ module math
|
|||
real(C_DOUBLE), value, intent(in) :: mu
|
||||
real(C_DOUBLE), optional, intent(in) :: phi
|
||||
end subroutine rotate_angle_c_intfc
|
||||
|
||||
subroutine spline(n, x, y, z) bind(C, name='spline_c')
|
||||
use ISO_C_BINDING
|
||||
implicit none
|
||||
integer(C_INT), value, intent(in) :: n
|
||||
real(C_DOUBLE), intent(in) :: x(n)
|
||||
real(C_DOUBLE), intent(in) :: y(n)
|
||||
real(C_DOUBLE), intent(in) :: z(n)
|
||||
end subroutine spline
|
||||
|
||||
function spline_interpolate(n, x, y, z, xint) &
|
||||
bind(C, name='spline_interpolate_c') result(yint)
|
||||
use ISO_C_BINDING
|
||||
implicit none
|
||||
integer(C_INT), value, intent(in) :: n
|
||||
real(C_DOUBLE), intent(in) :: x(n)
|
||||
real(C_DOUBLE), intent(in) :: y(n)
|
||||
real(C_DOUBLE), intent(in) :: z(n)
|
||||
real(C_DOUBLE), value, intent(in) :: xint
|
||||
real(C_DOUBLE) :: yint
|
||||
end function spline_interpolate
|
||||
|
||||
function spline_integrate(n, x, y, z, xa, xb) &
|
||||
bind(C, name='spline_integrate_c') result(s)
|
||||
use ISO_C_BINDING
|
||||
implicit none
|
||||
integer(C_INT), value, intent(in) :: n
|
||||
real(C_DOUBLE), intent(in) :: x(n)
|
||||
real(C_DOUBLE), intent(in) :: y(n)
|
||||
real(C_DOUBLE), intent(in) :: z(n)
|
||||
real(C_DOUBLE), value, intent(in) :: xa
|
||||
real(C_DOUBLE), value, intent(in) :: xb
|
||||
real(C_DOUBLE) :: s
|
||||
end function spline_integrate
|
||||
end interface
|
||||
|
||||
contains
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
@ -681,7 +681,7 @@ double normal_variate(double mean, double standard_deviation) {
|
|||
while ( true ) {
|
||||
double v1 = 2 * prn() - 1.;
|
||||
double v2 = 2 * prn() - 1.;
|
||||
|
||||
|
||||
double r = std::pow(v1, 2) + std::pow(v2, 2);
|
||||
double r2 = std::pow(r, 2);
|
||||
if (r2 < 1) {
|
||||
|
|
@ -697,7 +697,7 @@ double muir_spectrum(double e0, double m_rat, double kt) {
|
|||
// 8 in https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-05411-MS
|
||||
double sigma = std::sqrt(2.*e0*kt/m_rat);
|
||||
return normal_variate(e0, sigma);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
double watt_spectrum(double a, double b) {
|
||||
|
|
@ -753,7 +753,7 @@ void broaden_wmp_polynomials(double E, double dopp, int n, double factors[])
|
|||
}
|
||||
|
||||
|
||||
void spline_c(int n, const double x[], const double y[], double z[])
|
||||
void spline(int n, const double x[], const double y[], double z[])
|
||||
{
|
||||
double c_new[n-1];
|
||||
|
||||
|
|
@ -780,8 +780,8 @@ void spline_c(int n, const double x[], const double y[], double z[])
|
|||
}
|
||||
|
||||
|
||||
double spline_interpolate_c(int n, const double x[], const double y[],
|
||||
const double z[], double xint)
|
||||
double spline_interpolate(int n, const double x[], const double y[],
|
||||
const double z[], double xint)
|
||||
{
|
||||
// Find the lower bounding index in x of xint
|
||||
int i = n - 1;
|
||||
|
|
@ -801,8 +801,8 @@ double spline_interpolate_c(int n, const double x[], const double y[],
|
|||
}
|
||||
|
||||
|
||||
double spline_integrate_c(int n, const double x[], const double y[],
|
||||
const double z[], double xa, double xb)
|
||||
double spline_integrate(int n, const double x[], const double y[],
|
||||
const double z[], double xa, double xb)
|
||||
{
|
||||
// Find the lower bounding index in x of the lower limit of integration.
|
||||
int ia = n - 1;
|
||||
|
|
|
|||
|
|
@ -5,6 +5,7 @@
|
|||
#include "openmc/error.h"
|
||||
#include "openmc/hdf5_interface.h"
|
||||
#include "openmc/message_passing.h"
|
||||
#include "openmc/photon.h"
|
||||
#include "openmc/random_lcg.h"
|
||||
#include "openmc/search.h"
|
||||
#include "openmc/settings.h"
|
||||
|
|
@ -935,12 +936,14 @@ extern "C" bool multipole_in_range(Nuclide* nuc, double E)
|
|||
extern "C" void nuclides_clear() { data::nuclides.clear(); }
|
||||
|
||||
extern "C" NuclideMicroXS* micro_xs_ptr();
|
||||
extern "C" ElementMicroXS* micro_photon_xs_ptr();
|
||||
|
||||
void set_micro_xs()
|
||||
{
|
||||
#pragma omp parallel
|
||||
{
|
||||
simulation::micro_xs = micro_xs_ptr();
|
||||
simulation::micro_photon_xs = micro_photon_xs_ptr();
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -10,6 +10,7 @@ module particle_restart
|
|||
use nuclide_header, only: micro_xs, n_nuclides
|
||||
use output, only: print_particle
|
||||
use particle_header
|
||||
use photon_header, only: micro_photon_xs, n_elements
|
||||
use random_lcg, only: set_particle_seed
|
||||
use settings
|
||||
use simulation_header
|
||||
|
|
@ -45,7 +46,8 @@ contains
|
|||
|
||||
!$omp parallel
|
||||
allocate(micro_xs(n_nuclides))
|
||||
!$omp end parallel
|
||||
allocate(micro_photon_xs(n_elements))
|
||||
!$omp end parallel
|
||||
call set_micro_xs()
|
||||
|
||||
! Initialize the particle to be tracked
|
||||
|
|
|
|||
775
src/photon.cpp
775
src/photon.cpp
|
|
@ -1,6 +1,21 @@
|
|||
#include "openmc/photon.h"
|
||||
|
||||
#include "openmc/bremsstrahlung.h"
|
||||
#include "openmc/constants.h"
|
||||
#include "openmc/hdf5_interface.h"
|
||||
#include "openmc/nuclide.h"
|
||||
#include "openmc/particle.h"
|
||||
#include "openmc/random_lcg.h"
|
||||
#include "openmc/search.h"
|
||||
#include "openmc/settings.h"
|
||||
|
||||
#include "xtensor/xbuilder.hpp"
|
||||
#include "xtensor/xoperation.hpp"
|
||||
#include "xtensor/xview.hpp"
|
||||
|
||||
#include <array>
|
||||
#include <cmath>
|
||||
#include <tuple> // for tie
|
||||
|
||||
namespace openmc {
|
||||
|
||||
|
|
@ -10,29 +25,783 @@ namespace openmc {
|
|||
|
||||
namespace data {
|
||||
|
||||
xt::xtensor<double, 1> compton_profile_pz;
|
||||
|
||||
std::vector<PhotonInteraction> elements;
|
||||
|
||||
} // namespace data
|
||||
|
||||
namespace simulation {
|
||||
ElementMicroXS* micro_photon_xs;
|
||||
} // namespace simulation
|
||||
|
||||
//==============================================================================
|
||||
// PhotonInteraction implementation
|
||||
//==============================================================================
|
||||
|
||||
PhotonInteraction::PhotonInteraction(hid_t group)
|
||||
PhotonInteraction::PhotonInteraction(hid_t group, int i_element)
|
||||
: i_element_{i_element}
|
||||
{
|
||||
// Get name of nuclide from group, removing leading '/'
|
||||
name_ = object_name(group).substr(1);
|
||||
|
||||
// Get atomic number
|
||||
read_attribute(group, "Z", Z_);
|
||||
|
||||
// Determine number of energies and read energy grid
|
||||
read_dataset(group, "energy", energy_);
|
||||
|
||||
// Read coherent scattering
|
||||
hid_t rgroup = open_group(group, "coherent");
|
||||
read_dataset(rgroup, "xs", coherent_);
|
||||
|
||||
hid_t dset = open_dataset(rgroup, "integrated_scattering_factor");
|
||||
coherent_int_form_factor_ = Tabulated1D{dset};
|
||||
close_dataset(dset);
|
||||
|
||||
if (object_exists(group, "anomalous_real")) {
|
||||
dset = open_dataset(rgroup, "anomalous_real");
|
||||
coherent_anomalous_real_ = Tabulated1D{dset};
|
||||
close_dataset(dset);
|
||||
}
|
||||
|
||||
if (object_exists(group, "anomalous_imag")) {
|
||||
dset = open_dataset(rgroup, "anomalous_imag");
|
||||
coherent_anomalous_imag_ = Tabulated1D{dset};
|
||||
close_dataset(dset);
|
||||
}
|
||||
close_group(rgroup);
|
||||
|
||||
// Read incoherent scattering
|
||||
rgroup = open_group(group, "incoherent");
|
||||
read_dataset(rgroup, "xs", incoherent_);
|
||||
dset = open_dataset(rgroup, "scattering_factor");
|
||||
incoherent_form_factor_ = Tabulated1D{dset};
|
||||
close_dataset(dset);
|
||||
close_group(rgroup);
|
||||
|
||||
// Read pair production
|
||||
rgroup = open_group(group, "pair_production_electron");
|
||||
read_dataset(rgroup, "xs", pair_production_electron_);
|
||||
close_group(rgroup);
|
||||
|
||||
// Read pair production
|
||||
if (object_exists(group, "pair_production_nuclear")) {
|
||||
rgroup = open_group(group, "pair_production_nuclear");
|
||||
read_dataset(rgroup, "xs", pair_production_nuclear_);
|
||||
close_group(rgroup);
|
||||
} else {
|
||||
pair_production_nuclear_ = xt::zeros_like(energy_);
|
||||
}
|
||||
|
||||
// Read photoelectric
|
||||
rgroup = open_group(group, "photoelectric");
|
||||
read_dataset(rgroup, "xs", photoelectric_total_);
|
||||
close_group(rgroup);
|
||||
|
||||
// Read subshell photoionization cross section and atomic relaxation data
|
||||
rgroup = open_group(group, "subshells");
|
||||
std::vector<std::string> designators;
|
||||
read_attribute(rgroup, "designators", designators);
|
||||
auto n_shell = designators.size();
|
||||
for (int i = 0; i < n_shell; ++i) {
|
||||
const auto& designator {designators[i]};
|
||||
|
||||
// TODO: Move to ElectronSubshell constructor
|
||||
|
||||
// Add empty shell
|
||||
shells_.emplace_back();
|
||||
|
||||
// Create mapping from designator to index
|
||||
int j = 1;
|
||||
for (const auto& subshell : SUBSHELLS) {
|
||||
if (designator == subshell) {
|
||||
shell_map_[j] = i;
|
||||
shells_[i].index_subshell = j;
|
||||
break;
|
||||
}
|
||||
++j;
|
||||
}
|
||||
|
||||
// Read binding energy and number of electrons
|
||||
auto& shell {shells_.back()};
|
||||
hid_t tgroup = open_group(rgroup, designator.c_str());
|
||||
read_attribute(tgroup, "binding_energy", shell.binding_energy);
|
||||
read_attribute(tgroup, "num_electrons", shell.n_electrons);
|
||||
|
||||
// Read subshell cross section
|
||||
dset = open_dataset(tgroup, "xs");
|
||||
read_attribute(dset, "threshold_idx", j);
|
||||
shell.threshold = j;
|
||||
close_dataset(dset);
|
||||
read_dataset(tgroup, "xs", shell.cross_section);
|
||||
|
||||
auto& xs = shell.cross_section;
|
||||
xs = xt::where(xs > 0.0, xt::log(xs), -500.0);
|
||||
|
||||
if (object_exists(tgroup, "transitions")) {
|
||||
// Determine dimensions of transitions
|
||||
dset = open_dataset(tgroup, "transitions");
|
||||
auto dims = object_shape(dset);
|
||||
close_dataset(dset);
|
||||
|
||||
int n_transition = dims[0];
|
||||
shell.n_transitions = n_transition;
|
||||
if (n_transition > 0) {
|
||||
xt::xtensor<double, 2> matrix;
|
||||
read_dataset(tgroup, "transitions", matrix);
|
||||
|
||||
shell.transition_subshells = xt::view(matrix, xt::all(), xt::range(0, 2));
|
||||
shell.transition_energy = xt::view(matrix, xt::all(), 2);
|
||||
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);
|
||||
}
|
||||
close_group(rgroup);
|
||||
|
||||
// Determine number of electron shells
|
||||
rgroup = open_group(group, "compton_profiles");
|
||||
|
||||
// Read electron shell PDF and binding energies
|
||||
read_dataset(rgroup, "num_electrons", electron_pdf_);
|
||||
electron_pdf_ /= xt::sum(electron_pdf_);
|
||||
read_dataset(rgroup, "binding_energy", binding_energy_);
|
||||
|
||||
// Read Compton profiles
|
||||
read_dataset(rgroup, "J", profile_pdf_);
|
||||
|
||||
// Get Compton profile momentum grid. By deafult, an xtensor has a size of 1.
|
||||
// TODO: Update version of xtensor and change to 0
|
||||
if (data::compton_profile_pz.size() == 1) {
|
||||
read_dataset(rgroup, "pz", data::compton_profile_pz);
|
||||
}
|
||||
close_group(rgroup);
|
||||
|
||||
// Create Compton profile CDF
|
||||
auto n_profile = data::compton_profile_pz.size();
|
||||
profile_cdf_ = xt::empty<double>({n_shell, n_profile});
|
||||
for (int i = 0; i < n_shell; ++i) {
|
||||
double c = 0.0;
|
||||
profile_cdf_(i,0) = 0.0;
|
||||
for (int j = 0; j < n_profile - 1; ++j) {
|
||||
c += 0.5*(data::compton_profile_pz(j+1) - data::compton_profile_pz(j)) *
|
||||
(profile_pdf_(i,j) + profile_pdf_(i,j+1));
|
||||
profile_cdf_(i,j+1) = c;
|
||||
}
|
||||
}
|
||||
|
||||
// Calculate total pair production
|
||||
pair_production_total_ = pair_production_nuclear_ + pair_production_electron_;
|
||||
|
||||
if (settings::electron_treatment == ELECTRON_TTB) {
|
||||
// Read bremsstrahlung scaled DCS
|
||||
rgroup = open_group(group, "bremsstrahlung");
|
||||
read_dataset(rgroup, "dcs", dcs_);
|
||||
auto n_e = dcs_.shape()[0];
|
||||
auto n_k = dcs_.shape()[1];
|
||||
|
||||
// Get energy grids used for bremsstrahlung DCS and for stopping powers
|
||||
xt::xtensor<double, 1> electron_energy;
|
||||
read_dataset(rgroup, "electron_energy", electron_energy);
|
||||
if (data::ttb_k_grid.size() == 1) {
|
||||
read_dataset(rgroup, "photon_energy", data::ttb_k_grid);
|
||||
}
|
||||
close_group(rgroup);
|
||||
|
||||
// Read stopping power data
|
||||
if (Z_ < 99) {
|
||||
rgroup = open_group(group, "stopping_powers");
|
||||
read_dataset(rgroup, "s_collision", stopping_power_collision_);
|
||||
read_dataset(rgroup, "s_radiative", stopping_power_radiative_);
|
||||
read_attribute(rgroup, "I", I_);
|
||||
close_group(rgroup);
|
||||
}
|
||||
|
||||
// Truncate the bremsstrahlung data at the cutoff energy
|
||||
int photon = static_cast<int>(ParticleType::photon) - 1;
|
||||
const auto& E {electron_energy};
|
||||
double cutoff = settings::energy_cutoff[photon];
|
||||
if (cutoff > E(0)) {
|
||||
size_t i_grid = lower_bound_index(E.cbegin(), E.cend(),
|
||||
settings::energy_cutoff[photon]);
|
||||
|
||||
// calculate interpolation factor
|
||||
double f = (std::log(cutoff) - std::log(E(i_grid))) /
|
||||
(std::log(E(i_grid+1)) - std::log(E(i_grid)));
|
||||
|
||||
// Interpolate collision stopping power at the cutoff energy and truncate
|
||||
auto& s_col {stopping_power_collision_};
|
||||
double y = std::exp(std::log(s_col(i_grid)) + f*(std::log(s_col(i_grid+1)) -
|
||||
std::log(s_col(i_grid))));
|
||||
xt::xtensor<double, 1> frst {y};
|
||||
stopping_power_collision_ = xt::concatenate(xt::xtuple(
|
||||
frst, xt::view(s_col, xt::range(i_grid+1, n_e))));
|
||||
|
||||
// Interpolate radiative stopping power at the cutoff energy and truncate
|
||||
auto& s_rad {stopping_power_radiative_};
|
||||
y = std::exp(std::log(s_rad(i_grid)) + f*(std::log(s_rad(i_grid+1)) -
|
||||
std::log(s_rad(i_grid))));
|
||||
frst(0) = y;
|
||||
stopping_power_radiative_ = xt::concatenate(xt::xtuple(
|
||||
frst, xt::view(s_rad, xt::range(i_grid+1, n_e))));
|
||||
|
||||
// Interpolate bremsstrahlung DCS at the cutoff energy and truncate
|
||||
xt::xtensor<double, 2> dcs({n_e - i_grid, n_k});
|
||||
for (int i = 0; i < n_k; ++i) {
|
||||
y = std::exp(std::log(dcs_(i_grid,i)) +
|
||||
f*(std::log(dcs_(i_grid+1,i)) - std::log(dcs_(i_grid,i))));
|
||||
auto col_i = xt::view(dcs, xt::all(), i);
|
||||
col_i(0) = y;
|
||||
for (int j = i_grid + 1; j < n_e; ++j) {
|
||||
col_i(j - i_grid) = dcs_(j, i);
|
||||
}
|
||||
}
|
||||
dcs_ = dcs;
|
||||
|
||||
frst(0) = cutoff;
|
||||
electron_energy = xt::concatenate(xt::xtuple(
|
||||
frst, xt::view(electron_energy, xt::range(i_grid+1, n_e))));
|
||||
}
|
||||
|
||||
// Set incident particle energy grid
|
||||
// TODO: Change to zero when xtensor is updated
|
||||
if (data::ttb_e_grid.size() == 1) {
|
||||
data::ttb_e_grid = electron_energy;
|
||||
}
|
||||
}
|
||||
|
||||
// Take logarithm of energies and cross sections since they are log-log
|
||||
// interpolated
|
||||
energy_ = xt::log(energy_);
|
||||
coherent_ = xt::where(coherent_ > 0.0, xt::log(coherent_), -500.0);
|
||||
incoherent_ = xt::where(incoherent_ > 0.0, xt::log(incoherent_), -500.0);
|
||||
photoelectric_total_ = xt::where(photoelectric_total_ > 0.0,
|
||||
xt::log(photoelectric_total_), -500.0);
|
||||
pair_production_total_ = xt::where(pair_production_total_ > 0.0,
|
||||
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;
|
||||
}
|
||||
|
||||
void PhotonInteraction::calculate_xs(double E) const
|
||||
{
|
||||
// Perform binary search on the element energy grid in order to determine
|
||||
// which points to interpolate between
|
||||
int n_grid = energy_.size();
|
||||
double log_E = std::log(E);
|
||||
int i_grid;
|
||||
if (log_E <= energy_[0]) {
|
||||
i_grid = 0;
|
||||
} else if (log_E > energy_(n_grid - 1)) {
|
||||
i_grid = n_grid - 2;
|
||||
} else {
|
||||
// We use upper_bound_index here because sometimes photons are created with
|
||||
// energies that exactly match a grid point
|
||||
i_grid = upper_bound_index(energy_.cbegin(), energy_.cend(), log_E);
|
||||
}
|
||||
|
||||
// check for case where two energy points are the same
|
||||
if (energy_(i_grid) == energy_(i_grid+1)) ++i_grid;
|
||||
|
||||
// calculate interpolation factor
|
||||
double f = (log_E - energy_(i_grid)) / (energy_(i_grid+1) - energy_(i_grid));
|
||||
|
||||
auto& xs {simulation::micro_photon_xs[i_element_]};
|
||||
xs.index_grid = i_grid;
|
||||
xs.interp_factor = f;
|
||||
|
||||
// Calculate microscopic coherent cross section
|
||||
xs.coherent = std::exp(coherent_(i_grid) +
|
||||
f*(coherent_(i_grid+1) - coherent_(i_grid)));
|
||||
|
||||
// Calculate microscopic incoherent cross section
|
||||
xs.incoherent = std::exp(incoherent_(i_grid) +
|
||||
f*(incoherent_(i_grid+1) - incoherent_(i_grid)));
|
||||
|
||||
// Calculate microscopic photoelectric cross section
|
||||
xs.photoelectric = 0.0;
|
||||
for (const auto& shell : shells_) {
|
||||
// Check threshold of reaction
|
||||
int i_start = shell.threshold;
|
||||
if (i_grid < i_start) continue;
|
||||
|
||||
// Evaluation subshell photoionization cross section
|
||||
xs.photoelectric +=
|
||||
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)));
|
||||
}
|
||||
|
||||
// Calculate microscopic pair production cross section
|
||||
xs.pair_production = std::exp(
|
||||
pair_production_total_(i_grid) + f*(
|
||||
pair_production_total_(i_grid+1) -
|
||||
pair_production_total_(i_grid)));
|
||||
|
||||
// Calculate microscopic total cross section
|
||||
xs.total = xs.coherent + xs.incoherent + xs.photoelectric + xs.pair_production;
|
||||
xs.last_E = E;
|
||||
}
|
||||
|
||||
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};
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// Fortran compatibility
|
||||
//==============================================================================
|
||||
|
||||
extern "C" void photon_from_hdf5_c(hid_t group)
|
||||
extern "C" void photon_from_hdf5(hid_t group)
|
||||
{
|
||||
data::elements.emplace_back(group);
|
||||
data::elements.emplace_back(group, data::elements.size());
|
||||
|
||||
// Determine if minimum/maximum energy for this element is greater/less than
|
||||
// the previous
|
||||
const auto& element {data::elements.back()};
|
||||
if (element.energy_.size() >= 1) {
|
||||
// TODO: off-by-one
|
||||
int photon = static_cast<int>(ParticleType::photon) - 1;
|
||||
int n = element.energy_.size();
|
||||
data::energy_min[photon] = std::max(data::energy_min[photon],
|
||||
std::exp(element.energy_(1)));
|
||||
data::energy_max[photon] = std::min(data::energy_max[photon],
|
||||
std::exp(element.energy_(n - 1)));
|
||||
}
|
||||
}
|
||||
|
||||
extern "C" void photon_calculate_xs(int i_element, double E)
|
||||
{
|
||||
data::elements[i_element - 1].calculate_xs(E);
|
||||
}
|
||||
|
||||
extern "C" void free_memory_photon_c()
|
||||
{
|
||||
data::elements.clear();
|
||||
data::compton_profile_pz.resize({0});
|
||||
data::ttb_e_grid.resize({0});
|
||||
data::ttb_k_grid.resize({0});
|
||||
}
|
||||
|
||||
} // namespace openmc
|
||||
|
|
|
|||
|
|
@ -1,520 +1,57 @@
|
|||
module photon_header
|
||||
|
||||
use algorithm, only: binary_search
|
||||
use, intrinsic :: ISO_C_BINDING
|
||||
|
||||
use constants
|
||||
use dict_header, only: DictIntInt, DictCharInt
|
||||
use endf_header, only: Tabulated1D
|
||||
use hdf5_interface
|
||||
use nuclide_header, only: nuclides
|
||||
use settings
|
||||
use dict_header, only: DictCharInt
|
||||
|
||||
real(8), allocatable :: compton_profile_pz(:)
|
||||
real(8), allocatable :: ttb_e_grid(:) ! energy T of incident electron
|
||||
real(8), allocatable :: ttb_k_grid(:) ! reduced energy W/T of emitted photon
|
||||
|
||||
type ElectronSubshell
|
||||
integer :: index_subshell ! index in SUBSHELLS
|
||||
integer :: threshold
|
||||
real(8) :: n_electrons
|
||||
real(8) :: binding_energy
|
||||
real(8), allocatable :: cross_section(:)
|
||||
|
||||
! Transition data
|
||||
integer :: n_transitions
|
||||
integer, allocatable :: transition_subshells(:,:)
|
||||
real(8), allocatable :: transition_energy(:)
|
||||
real(8), allocatable :: transition_probability(:)
|
||||
end type ElectronSubshell
|
||||
|
||||
type PhotonInteraction
|
||||
character(3) :: name ! atomic symbol, e.g. 'Zr'
|
||||
integer :: Z ! atomic number
|
||||
|
||||
! Microscopic cross sections
|
||||
real(8), allocatable :: energy(:)
|
||||
real(8), allocatable :: coherent(:)
|
||||
real(8), allocatable :: incoherent(:)
|
||||
real(8), allocatable :: photoelectric_total(:)
|
||||
real(8), allocatable :: pair_production_total(:)
|
||||
real(8), allocatable :: pair_production_electron(:)
|
||||
real(8), allocatable :: pair_production_nuclear(:)
|
||||
|
||||
! Form factors
|
||||
type(Tabulated1D) :: incoherent_form_factor
|
||||
type(Tabulated1D) :: coherent_int_form_factor
|
||||
type(Tabulated1D) :: coherent_anomalous_real
|
||||
type(Tabulated1D) :: coherent_anomalous_imag
|
||||
|
||||
! Photoionization and atomic relaxation data
|
||||
type(DictIntInt) :: shell_dict ! Given a shell designator, e.g. 3, this
|
||||
! dictionary gives an index in shells(:)
|
||||
type(ElectronSubshell), allocatable :: shells(:)
|
||||
|
||||
! Compton profile data
|
||||
real(8), allocatable :: profile_pdf(:,:)
|
||||
real(8), allocatable :: profile_cdf(:,:)
|
||||
real(8), allocatable :: binding_energy(:)
|
||||
real(8), allocatable :: electron_pdf(:)
|
||||
|
||||
! Stopping power data
|
||||
real(8) :: I ! mean excitation energy
|
||||
real(8), allocatable :: stopping_power_collision(:)
|
||||
real(8), allocatable :: stopping_power_radiative(:)
|
||||
|
||||
! Bremsstrahlung scaled DCS
|
||||
real(8), allocatable :: dcs(:,:)
|
||||
|
||||
contains
|
||||
procedure :: from_hdf5 => photon_from_hdf5
|
||||
procedure :: calculate_xs => photon_calculate_xs
|
||||
end type PhotonInteraction
|
||||
|
||||
type BremsstrahlungData
|
||||
real(8), allocatable :: pdf(:,:) ! Bremsstrahlung energy PDF
|
||||
real(8), allocatable :: cdf(:,:) ! Bremsstrahlung energy CDF
|
||||
real(8), allocatable :: yield(:) ! Photon number yield
|
||||
end type BremsstrahlungData
|
||||
|
||||
type Bremsstrahlung
|
||||
type(BremsstrahlungData) :: electron
|
||||
type(BremsstrahlungData) :: positron
|
||||
end type Bremsstrahlung
|
||||
|
||||
type(PhotonInteraction), allocatable, target :: elements(:) ! Photon cross sections
|
||||
integer :: n_elements ! Number of photon cross section tables
|
||||
|
||||
type(DictCharInt) :: element_dict
|
||||
|
||||
type(Bremsstrahlung), allocatable, target :: ttb(:) ! Bremsstrahlung data
|
||||
|
||||
!===============================================================================
|
||||
! ELEMENTMICROXS contains cached microscopic photon cross sections for a
|
||||
! particular element at the current energy
|
||||
!===============================================================================
|
||||
|
||||
type ElementMicroXS
|
||||
integer :: index_grid ! index on element energy grid
|
||||
real(8) :: last_E = ZERO ! last evaluated energy
|
||||
real(8) :: interp_factor ! interpolation factor on energy grid
|
||||
real(8) :: total ! microscropic total photon xs
|
||||
real(8) :: coherent ! microscopic coherent xs
|
||||
real(8) :: incoherent ! microscopic incoherent xs
|
||||
real(8) :: photoelectric ! microscopic photoelectric xs
|
||||
real(8) :: pair_production ! microscopic pair production xs
|
||||
type, bind(C) :: ElementMicroXS
|
||||
integer(C_INT) :: index_grid ! index on element energy grid
|
||||
real(C_DOUBLE) :: last_E = ZERO ! last evaluated energy
|
||||
real(C_DOUBLE) :: interp_factor ! interpolation factor on energy grid
|
||||
real(C_DOUBLE) :: total ! microscropic total photon xs
|
||||
real(C_DOUBLE) :: coherent ! microscopic coherent xs
|
||||
real(C_DOUBLE) :: incoherent ! microscopic incoherent xs
|
||||
real(C_DOUBLE) :: photoelectric ! microscopic photoelectric xs
|
||||
real(C_DOUBLE) :: pair_production ! microscopic pair production xs
|
||||
end type ElementMicroXS
|
||||
|
||||
type(ElementMicroXS), allocatable :: micro_photon_xs(:) ! Cache for each element
|
||||
type(ElementMicroXS), allocatable, target :: micro_photon_xs(:) ! Cache for each element
|
||||
!$omp threadprivate(micro_photon_xs)
|
||||
|
||||
contains
|
||||
|
||||
subroutine photon_from_hdf5(this, group_id)
|
||||
class(PhotonInteraction), intent(inout) :: this
|
||||
integer(HID_T), intent(in) :: group_id
|
||||
|
||||
integer :: i, j
|
||||
integer(HID_T) :: rgroup, tgroup
|
||||
integer(HID_T) :: dset_id
|
||||
integer(HSIZE_T) :: dims(1), dims2(2)
|
||||
integer :: n_energy
|
||||
integer :: n_shell
|
||||
integer :: n_profile
|
||||
integer :: n_transition
|
||||
integer :: n_k
|
||||
integer :: n_e
|
||||
character(3), allocatable :: designators(:)
|
||||
real(8) :: c
|
||||
real(8) :: f
|
||||
real(8) :: y
|
||||
real(8), allocatable :: electron_energy(:)
|
||||
real(8), allocatable :: matrix(:,:)
|
||||
real(8), allocatable :: dcs(:,:)
|
||||
|
||||
interface
|
||||
subroutine photon_from_hdf5_c(group) bind(C)
|
||||
import HID_T
|
||||
integer(HID_T), value :: group
|
||||
end subroutine
|
||||
end interface
|
||||
|
||||
! Read element data on C++ side
|
||||
call photon_from_hdf5_c(group_id)
|
||||
|
||||
! Get name of nuclide from group
|
||||
this % name = get_name(group_id)
|
||||
|
||||
! Get rid of leading '/'
|
||||
this % name = trim(this % name(2:))
|
||||
|
||||
! Get atomic number
|
||||
call read_attribute(this % Z, group_id, 'Z')
|
||||
|
||||
! Determine number of energies and read energy grid
|
||||
dset_id = open_dataset(group_id, 'energy')
|
||||
call get_shape(dset_id, dims)
|
||||
n_energy = int(dims(1), 4)
|
||||
allocate(this % energy(dims(1)))
|
||||
call read_dataset(this % energy, dset_id)
|
||||
call close_dataset(dset_id)
|
||||
|
||||
! Allocate arrays
|
||||
allocate(this % coherent(n_energy))
|
||||
allocate(this % incoherent(n_energy))
|
||||
allocate(this % pair_production_total(n_energy))
|
||||
allocate(this % pair_production_nuclear(n_energy))
|
||||
allocate(this % pair_production_electron(n_energy))
|
||||
allocate(this % photoelectric_total(n_energy))
|
||||
|
||||
! Read coherent scattering
|
||||
rgroup = open_group(group_id, 'coherent')
|
||||
call read_dataset(this % coherent, rgroup, 'xs')
|
||||
|
||||
dset_id = open_dataset(rgroup, 'integrated_scattering_factor')
|
||||
call this % coherent_int_form_factor % from_hdf5(dset_id)
|
||||
call close_dataset(dset_id)
|
||||
|
||||
if (object_exists(group_id, 'anomalous_real')) then
|
||||
dset_id = open_dataset(rgroup, 'anomalous_real')
|
||||
call this % coherent_anomalous_real % from_hdf5(dset_id)
|
||||
call close_dataset(dset_id)
|
||||
end if
|
||||
|
||||
if (object_exists(group_id, 'anomalous_imag')) then
|
||||
dset_id = open_dataset(rgroup, 'anomalous_imag')
|
||||
call this % coherent_anomalous_imag % from_hdf5(dset_id)
|
||||
call close_dataset(dset_id)
|
||||
call close_group(rgroup)
|
||||
end if
|
||||
|
||||
! Read incoherent scattering
|
||||
rgroup = open_group(group_id, 'incoherent')
|
||||
call read_dataset(this % incoherent, rgroup, 'xs')
|
||||
dset_id = open_dataset(rgroup, 'scattering_factor')
|
||||
call this % incoherent_form_factor % from_hdf5(dset_id)
|
||||
call close_dataset(dset_id)
|
||||
call close_group(rgroup)
|
||||
|
||||
! Read pair production
|
||||
rgroup = open_group(group_id, 'pair_production_electron')
|
||||
call read_dataset(this % pair_production_electron, rgroup, 'xs')
|
||||
call close_group(rgroup)
|
||||
|
||||
! Read pair production
|
||||
if (object_exists(group_id, 'pair_production_nuclear')) then
|
||||
rgroup = open_group(group_id, 'pair_production_nuclear')
|
||||
call read_dataset(this % pair_production_nuclear, rgroup, 'xs')
|
||||
call close_group(rgroup)
|
||||
else
|
||||
this % pair_production_nuclear(:) = ZERO
|
||||
end if
|
||||
|
||||
! Read photoelectric
|
||||
rgroup = open_group(group_id, 'photoelectric')
|
||||
call read_dataset(this % photoelectric_total, rgroup, 'xs')
|
||||
call close_group(rgroup)
|
||||
|
||||
! Read subshell photoionization cross section and atomic relaxation data
|
||||
rgroup = open_group(group_id, 'subshells')
|
||||
call read_attribute(designators, rgroup, 'designators')
|
||||
n_shell = size(designators)
|
||||
allocate(this % shells(n_shell))
|
||||
do i = 1, n_shell
|
||||
! Create mapping from designator to index
|
||||
do j = 1, size(SUBSHELLS)
|
||||
if (designators(i) == SUBSHELLS(j)) then
|
||||
call this % shell_dict % set(j, i)
|
||||
this % shells(i) % index_subshell = j
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
|
||||
! Read binding energy and number of electrons
|
||||
tgroup = open_group(rgroup, trim(designators(i)))
|
||||
call read_attribute(this % shells(i) % binding_energy, tgroup, &
|
||||
'binding_energy')
|
||||
call read_attribute(this % shells(i) % n_electrons, tgroup, &
|
||||
'num_electrons')
|
||||
|
||||
! Read subshell cross section
|
||||
dset_id = open_dataset(tgroup, 'xs')
|
||||
call read_attribute(j, dset_id, 'threshold_idx')
|
||||
this % shells(i) % threshold = j
|
||||
allocate(this % shells(i) % cross_section(n_energy - j))
|
||||
call read_dataset(this % shells(i) % cross_section, dset_id)
|
||||
call close_dataset(dset_id)
|
||||
where (this % shells(i) % cross_section > ZERO)
|
||||
this % shells(i) % cross_section = log(this % shells(i) % cross_section)
|
||||
elsewhere
|
||||
this % shells(i) % cross_section = -500.0_8
|
||||
end where
|
||||
|
||||
if (object_exists(tgroup, 'transitions')) then
|
||||
dset_id = open_dataset(tgroup, 'transitions')
|
||||
call get_shape(dset_id, dims2)
|
||||
n_transition = int(dims2(2), 4)
|
||||
this % shells(i) % n_transitions = n_transition
|
||||
if (n_transition > 0) then
|
||||
allocate(this % shells(i) % transition_subshells(2, n_transition))
|
||||
allocate(this % shells(i) % transition_energy(n_transition))
|
||||
allocate(this % shells(i) % transition_probability(n_transition))
|
||||
|
||||
allocate(matrix(dims2(1), dims2(2)))
|
||||
call read_dataset(matrix, dset_id)
|
||||
|
||||
this % shells(i) % transition_subshells(:,:) = int(matrix(1:2, :), 4)
|
||||
this % shells(i) % transition_energy(:) = matrix(3, :)
|
||||
this % shells(i) % transition_probability(:) = matrix(4, :) &
|
||||
/ sum(matrix(4, :))
|
||||
deallocate(matrix)
|
||||
end if
|
||||
call close_dataset(dset_id)
|
||||
else
|
||||
this % shells(i) % n_transitions = 0
|
||||
end if
|
||||
call close_group(tgroup)
|
||||
end do
|
||||
call close_group(rgroup)
|
||||
deallocate(designators)
|
||||
|
||||
! Determine number of electron shells
|
||||
rgroup = open_group(group_id, 'compton_profiles')
|
||||
|
||||
! Determine number of shells
|
||||
dset_id = open_dataset(rgroup, 'num_electrons')
|
||||
call get_shape(dset_id, dims)
|
||||
n_shell = int(dims(1), 4)
|
||||
|
||||
! Read electron shell PDF and binding energies
|
||||
allocate(this % electron_pdf(n_shell), this % binding_energy(n_shell))
|
||||
call read_dataset(this % electron_pdf, dset_id)
|
||||
call close_dataset(dset_id)
|
||||
call read_dataset(this % binding_energy, rgroup, 'binding_energy')
|
||||
this % electron_pdf(:) = this % electron_pdf / sum(this % electron_pdf)
|
||||
|
||||
! Read Compton profiles
|
||||
dset_id = open_dataset(rgroup, 'J')
|
||||
call get_shape(dset_id, dims2)
|
||||
n_profile = int(dims2(1), 4)
|
||||
allocate(this % profile_pdf(n_profile, n_shell))
|
||||
call read_dataset(this % profile_pdf, dset_id)
|
||||
call close_dataset(dset_id)
|
||||
|
||||
! Get Compton profile momentum grid
|
||||
if (.not. allocated(compton_profile_pz)) then
|
||||
allocate(compton_profile_pz(n_profile))
|
||||
call read_dataset(compton_profile_pz, rgroup, 'pz')
|
||||
end if
|
||||
call close_group(rgroup)
|
||||
|
||||
! Create Compton profile CDF
|
||||
allocate(this % profile_cdf(n_profile, n_shell))
|
||||
do i = 1, n_shell
|
||||
c = ZERO
|
||||
this % profile_cdf(1,i) = ZERO
|
||||
do j = 1, n_profile - 1
|
||||
c = c + HALF*(compton_profile_pz(j+1) - compton_profile_pz(j)) * &
|
||||
(this%profile_pdf(j,i) + this%profile_pdf(j+1,i))
|
||||
this % profile_cdf(j+1,i) = c
|
||||
end do
|
||||
end do
|
||||
|
||||
! Calculate total pair production
|
||||
this % pair_production_total(:) = this % pair_production_nuclear + &
|
||||
this % pair_production_electron
|
||||
|
||||
if (electron_treatment == ELECTRON_TTB) then
|
||||
! Read bremsstrahlung scaled DCS
|
||||
rgroup = open_group(group_id, 'bremsstrahlung')
|
||||
dset_id = open_dataset(rgroup, 'dcs')
|
||||
call get_shape(dset_id, dims2)
|
||||
n_k = int(dims2(1), 4)
|
||||
n_e = int(dims2(2), 4)
|
||||
allocate(this % dcs(n_k, n_e))
|
||||
call read_dataset(this % dcs, dset_id)
|
||||
call close_dataset(dset_id)
|
||||
|
||||
! Get energy grids used for bremsstrahlung DCS and for stopping powers
|
||||
allocate(electron_energy(n_e))
|
||||
call read_dataset(electron_energy, rgroup, 'electron_energy')
|
||||
if (.not. allocated(ttb_k_grid)) then
|
||||
allocate(ttb_k_grid(n_k))
|
||||
call read_dataset(ttb_k_grid, rgroup, 'photon_energy')
|
||||
end if
|
||||
call close_group(rgroup)
|
||||
|
||||
! Read stopping power data
|
||||
if (this % Z < 99) then
|
||||
rgroup = open_group(group_id, 'stopping_powers')
|
||||
allocate(this % stopping_power_collision(n_e))
|
||||
allocate(this % stopping_power_radiative(n_e))
|
||||
call read_dataset(this % stopping_power_collision, rgroup, 's_collision')
|
||||
call read_dataset(this % stopping_power_radiative, rgroup, 's_radiative')
|
||||
call read_attribute(this % I, rgroup, 'I')
|
||||
call close_group(rgroup)
|
||||
end if
|
||||
|
||||
! Truncate the bremsstrahlung data at the cutoff energy
|
||||
if (energy_cutoff(PHOTON) > electron_energy(1)) then
|
||||
i_grid = binary_search(electron_energy, n_e, energy_cutoff(PHOTON))
|
||||
|
||||
! calculate interpolation factor
|
||||
f = (log(energy_cutoff(PHOTON)) - log(electron_energy(i_grid))) / &
|
||||
(log(electron_energy(i_grid+1)) - log(electron_energy(i_grid)))
|
||||
|
||||
! Interpolate collision stopping power at the cutoff energy and
|
||||
! truncate
|
||||
y = exp(log(this % stopping_power_collision(i_grid)) + &
|
||||
f*(log(this % stopping_power_collision(i_grid+1)) - &
|
||||
log(this % stopping_power_collision(i_grid))))
|
||||
this % stopping_power_collision = &
|
||||
[y, this % stopping_power_collision(i_grid+1:n_e)]
|
||||
|
||||
! Interpolate radiative stopping power at the cutoff energy and
|
||||
! truncate
|
||||
y = exp(log(this % stopping_power_radiative(i_grid)) + &
|
||||
f*(log(this % stopping_power_radiative(i_grid+1)) - &
|
||||
log(this % stopping_power_radiative(i_grid))))
|
||||
this % stopping_power_radiative = &
|
||||
[y, this % stopping_power_radiative(i_grid+1:n_e)]
|
||||
|
||||
! Interpolate bremsstrahlung DCS at the cutoff energy and truncate
|
||||
allocate(dcs(n_k, n_e-i_grid+1))
|
||||
do i = 1, n_k
|
||||
y = exp(log(this % dcs(i,i_grid)) + &
|
||||
f*(log(this % dcs(i,i_grid+1)) - log(this % dcs(i,i_grid))))
|
||||
dcs(i,:) = [y, this % dcs(i,i_grid+1:n_e)]
|
||||
end do
|
||||
call move_alloc(dcs, this % dcs)
|
||||
|
||||
electron_energy = [energy_cutoff(PHOTON), electron_energy(i_grid+1:n_e)]
|
||||
end if
|
||||
|
||||
! Set incident particle energy grid
|
||||
if (.not. allocated(ttb_e_grid)) then
|
||||
call move_alloc(electron_energy, ttb_e_grid)
|
||||
end if
|
||||
end if
|
||||
|
||||
! Take logarithm of energies and cross sections since they are log-log
|
||||
! interpolated
|
||||
this % energy = log(this % energy)
|
||||
|
||||
where (this % coherent > ZERO)
|
||||
this % coherent = log(this % coherent)
|
||||
elsewhere
|
||||
this % coherent = -500.0_8
|
||||
end where
|
||||
|
||||
where (this % incoherent > ZERO)
|
||||
this % incoherent = log(this % incoherent)
|
||||
elsewhere
|
||||
this % incoherent = -500.0_8
|
||||
end where
|
||||
|
||||
where (this % photoelectric_total > ZERO)
|
||||
this % photoelectric_total = log(this % photoelectric_total)
|
||||
elsewhere
|
||||
this % photoelectric_total = -500.0_8
|
||||
end where
|
||||
|
||||
where (this % pair_production_total > ZERO)
|
||||
this % pair_production_total = log(this % pair_production_total)
|
||||
elsewhere
|
||||
this % pair_production_total = -500.0_8
|
||||
end where
|
||||
|
||||
end subroutine photon_from_hdf5
|
||||
|
||||
!===============================================================================
|
||||
! CALCULATE_ELEMENT_XS determines microscopic photon cross sections for an
|
||||
! element of a given index in the elements array at the energy of the given
|
||||
! particle
|
||||
!===============================================================================
|
||||
|
||||
subroutine photon_calculate_xs(this, E, xs)
|
||||
class(PhotonInteraction), intent(in) :: this ! index into elements array
|
||||
real(8), intent(in) :: E ! energy
|
||||
type(ElementMicroXS), intent(inout) :: xs
|
||||
|
||||
integer :: i_grid ! index on element energy grid
|
||||
integer :: i_shell ! index in subshells
|
||||
integer :: i_start ! threshold index
|
||||
integer :: n_grid ! number of grid points
|
||||
real(8) :: f ! interp factor on element energy grid
|
||||
real(8) :: log_E ! logarithm of the energy
|
||||
|
||||
! Perform binary search on the element energy grid in order to determine
|
||||
! which points to interpolate between
|
||||
n_grid = size(this % energy)
|
||||
log_E = log(E)
|
||||
if (log_E <= this % energy(1)) then
|
||||
i_grid = 1
|
||||
elseif (log_E > this % energy(n_grid)) then
|
||||
i_grid = n_grid - 1
|
||||
else
|
||||
i_grid = binary_search(this % energy, n_grid, log_E)
|
||||
end if
|
||||
|
||||
! check for case where two energy points are the same
|
||||
if (this % energy(i_grid) == this % energy(i_grid+1)) i_grid = i_grid + 1
|
||||
|
||||
! calculate interpolation factor
|
||||
f = (log_E - this % energy(i_grid)) / &
|
||||
(this % energy(i_grid+1) - this % energy(i_grid))
|
||||
|
||||
xs % index_grid = i_grid
|
||||
xs % interp_factor = f
|
||||
|
||||
! Calculate microscopic coherent cross section
|
||||
xs % coherent = exp(this % coherent(i_grid) + f * &
|
||||
(this % coherent(i_grid+1) - this % coherent(i_grid)))
|
||||
|
||||
! Calculate microscopic incoherent cross section
|
||||
xs % incoherent = exp(this % incoherent(i_grid) + &
|
||||
f*(this % incoherent(i_grid+1) - this % incoherent(i_grid)))
|
||||
|
||||
! Calculate microscopic photoelectric cross section
|
||||
xs % photoelectric = ZERO
|
||||
do i_shell = 1, size(this % shells)
|
||||
! Check threshold of reaction
|
||||
i_start = this % shells(i_shell) % threshold
|
||||
if (i_grid <= i_start) cycle
|
||||
|
||||
! Evaluation subshell photoionization cross section
|
||||
xs % photoelectric = xs % photoelectric + &
|
||||
exp(this % shells(i_shell) % cross_section(i_grid-i_start) + &
|
||||
f*(this % shells(i_shell) % cross_section(i_grid+1-i_start) - &
|
||||
this % shells(i_shell) % cross_section(i_grid-i_start)))
|
||||
end do
|
||||
|
||||
! Calculate microscopic pair production cross section
|
||||
xs % pair_production = exp(&
|
||||
this % pair_production_total(i_grid) + f*(&
|
||||
this % pair_production_total(i_grid+1) - &
|
||||
this % pair_production_total(i_grid)))
|
||||
|
||||
! Calculate microscopic total cross section
|
||||
xs % total = xs % coherent + xs % incoherent + xs % photoelectric + &
|
||||
xs % pair_production
|
||||
|
||||
xs % last_E = E
|
||||
|
||||
end subroutine photon_calculate_xs
|
||||
|
||||
!===============================================================================
|
||||
! FREE_MEMORY_PHOTON deallocates/resets global variables in this module
|
||||
!===============================================================================
|
||||
|
||||
subroutine free_memory_photon()
|
||||
interface
|
||||
subroutine free_memory_photon_c() bind(C)
|
||||
end subroutine
|
||||
end interface
|
||||
|
||||
! Deallocate photon cross section data
|
||||
if (allocated(elements)) deallocate(elements)
|
||||
if (allocated(compton_profile_pz)) deallocate(compton_profile_pz)
|
||||
n_elements = 0
|
||||
call element_dict % clear()
|
||||
|
||||
! Clear TTB-related arrays
|
||||
if (allocated(ttb_e_grid)) deallocate(ttb_e_grid)
|
||||
if (allocated(ttb)) deallocate(ttb)
|
||||
end subroutine free_memory_photon
|
||||
|
||||
function micro_photon_xs_ptr() result(ptr) bind(C)
|
||||
type(C_PTR) :: ptr
|
||||
if (size(micro_photon_xs) > 0) then
|
||||
ptr = C_LOC(micro_photon_xs(1))
|
||||
else
|
||||
ptr = C_NULL_PTR
|
||||
end if
|
||||
end function
|
||||
|
||||
end module photon_header
|
||||
|
|
|
|||
|
|
@ -1,643 +0,0 @@
|
|||
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 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
|
||||
!===============================================================================
|
||||
|
||||
subroutine thick_target_bremsstrahlung(p, E_lost) bind(C)
|
||||
type(Particle), intent(inout) :: p
|
||||
real(C_DOUBLE), intent(out) :: E_lost
|
||||
|
||||
integer :: i, j
|
||||
integer :: i_e, i_w
|
||||
integer :: n
|
||||
integer :: n_e
|
||||
real(8) :: a
|
||||
real(8) :: f
|
||||
real(8) :: e, e_l, e_r
|
||||
real(8) :: y, y_l, y_r
|
||||
real(8) :: w, w_l, w_r
|
||||
real(8) :: p_l, p_r
|
||||
real(8) :: c, c_l, c_max
|
||||
type(BremsstrahlungData), pointer :: mat
|
||||
|
||||
if (p % material == MATERIAL_VOID) return
|
||||
|
||||
if (p % E < energy_cutoff(PHOTON)) return
|
||||
|
||||
! Get bremsstrahlung data for this material and particle type
|
||||
if (p % type == POSITRON) then
|
||||
mat => ttb(p % material) % positron
|
||||
else
|
||||
mat => ttb(p % material) % electron
|
||||
end if
|
||||
|
||||
e = log(p % E)
|
||||
n_e = size(ttb_e_grid)
|
||||
|
||||
! Find the lower bounding index of the incident electron energy
|
||||
j = binary_search(ttb_e_grid, n_e, e)
|
||||
if (j == n_e) j = j - 1
|
||||
|
||||
! Get the interpolation bounds
|
||||
e_l = ttb_e_grid(j)
|
||||
e_r = ttb_e_grid(j+1)
|
||||
y_l = mat % yield(j)
|
||||
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
|
||||
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
|
||||
y = exp(y_l + (y_r - y_l)*f)
|
||||
|
||||
! Sample number of secondary bremsstrahlung photons
|
||||
n = int(y + prn())
|
||||
|
||||
E_lost = ZERO
|
||||
if (n == 0) return
|
||||
|
||||
! Sample index of the tabulated PDF in the energy grid, j or j+1
|
||||
if (prn() <= f .or. j == 1) then
|
||||
i_e = j + 1
|
||||
|
||||
! Interpolate the maximum value of the CDF at the incoming particle
|
||||
! energy on a log-log scale
|
||||
p_l = mat % pdf(i_e-1, i_e)
|
||||
p_r = mat % pdf(i_e, i_e)
|
||||
c_l = mat % cdf(i_e-1, i_e)
|
||||
a = log(p_r/p_l)/(e_r - e_l) + ONE
|
||||
c_max = c_l + exp(e_l)*p_l/a*(exp(a*(e - e_l)) - ONE)
|
||||
else
|
||||
i_e = j
|
||||
|
||||
! Maximum value of the CDF
|
||||
c_max = mat % cdf(i_e, i_e)
|
||||
end if
|
||||
|
||||
! Sample the energies of the emitted photons
|
||||
do i = 1, n
|
||||
! Generate a random number r and determine the index i for which
|
||||
! cdf(i) <= r*cdf,max <= cdf(i+1)
|
||||
c = prn()*c_max
|
||||
i_w = binary_search(mat % cdf(:i_e,i_e), i_e, c)
|
||||
|
||||
! Sample the photon energy
|
||||
w_l = ttb_e_grid(i_w)
|
||||
w_r = ttb_e_grid(i_w+1)
|
||||
p_l = mat % pdf(i_w, i_e)
|
||||
p_r = mat % pdf(i_w+1, i_e)
|
||||
c_l = mat % cdf(i_w, i_e)
|
||||
a = log(p_r/p_l)/(w_r - w_l) + ONE
|
||||
w = exp(w_l)*(a*(c - c_l)/(exp(w_l)*p_l) + ONE)**(ONE/a)
|
||||
|
||||
if (w > energy_cutoff(PHOTON)) then
|
||||
! Create secondary photon
|
||||
call particle_create_secondary(p, p % coord(1) % uvw, w, PHOTON, &
|
||||
run_ce=.true._C_BOOL)
|
||||
E_lost = E_lost + w
|
||||
end if
|
||||
end do
|
||||
|
||||
end subroutine thick_target_bremsstrahlung
|
||||
|
||||
end module photon_physics
|
||||
251
src/physics.F90
251
src/physics.F90
|
|
@ -1,251 +0,0 @@
|
|||
module physics
|
||||
|
||||
use constants
|
||||
use error, only: fatal_error, warning, write_message
|
||||
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
|
||||
|
||||
interface
|
||||
subroutine collision(p) bind(C)
|
||||
import Particle
|
||||
type(Particle), intent(inout) :: p
|
||||
end subroutine
|
||||
end interface
|
||||
|
||||
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)
|
||||
type(Particle), intent(in) :: p
|
||||
integer :: i_element
|
||||
|
||||
integer :: i
|
||||
real(8) :: prob
|
||||
real(8) :: cutoff
|
||||
real(8) :: atom_density ! atom density of nuclide in atom/b-cm
|
||||
real(8) :: sigma ! microscopic total xs for nuclide
|
||||
|
||||
associate (mat => materials(p % material))
|
||||
! Sample cumulative distribution function
|
||||
cutoff = prn() * material_xs % total
|
||||
|
||||
i = 0
|
||||
prob = ZERO
|
||||
do while (prob < cutoff)
|
||||
i = i + 1
|
||||
|
||||
! Check to make sure that a nuclide was sampled
|
||||
if (i > mat % n_nuclides) then
|
||||
call particle_write_restart(p)
|
||||
call fatal_error("Did not sample any element during collision.")
|
||||
end if
|
||||
|
||||
! Find atom density
|
||||
i_element = mat % element(i)
|
||||
atom_density = mat % atom_density(i)
|
||||
|
||||
! 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
|
||||
|
||||
end function sample_element
|
||||
|
||||
end module physics
|
||||
328
src/physics.cpp
328
src/physics.cpp
|
|
@ -1,6 +1,7 @@
|
|||
#include "openmc/physics.h"
|
||||
|
||||
#include "openmc/bank.h"
|
||||
#include "openmc/bremsstrahlung.h"
|
||||
#include "openmc/constants.h"
|
||||
#include "openmc/eigenvalue.h"
|
||||
#include "openmc/error.h"
|
||||
|
|
@ -214,154 +215,167 @@ 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);
|
||||
// TODO: off-by-one
|
||||
p->event_nuclide = i_element + 1;
|
||||
const auto& micro {simulation::micro_photon_xs[i_element]};
|
||||
const auto& element {data::elements[i_element]};
|
||||
|
||||
// // 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
|
||||
int i_grid = micro.index_grid;
|
||||
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)
|
||||
{
|
||||
|
|
@ -369,7 +383,7 @@ void sample_electron_reaction(Particle* p)
|
|||
|
||||
if (settings::electron_treatment == ELECTRON_TTB) {
|
||||
double E_lost;
|
||||
thick_target_bremsstrahlung(p, &E_lost);
|
||||
thick_target_bremsstrahlung(*p, &E_lost);
|
||||
}
|
||||
|
||||
p->E = 0.0;
|
||||
|
|
@ -382,7 +396,7 @@ void sample_positron_reaction(Particle* p)
|
|||
|
||||
if (settings::electron_treatment == ELECTRON_TTB) {
|
||||
double E_lost;
|
||||
thick_target_bremsstrahlung(p, &E_lost);
|
||||
thick_target_bremsstrahlung(*p, &E_lost);
|
||||
}
|
||||
|
||||
// Sample angle isotropically
|
||||
|
|
@ -465,37 +479,43 @@ void sample_nuclide(const Particle* p, int mt, int* i_nuclide, int* i_nuc_mat)
|
|||
}
|
||||
}
|
||||
|
||||
// TODO: Finish converting photon physics functions
|
||||
int sample_element(Particle* p)
|
||||
{
|
||||
// Sample cumulative distribution function
|
||||
double cutoff = prn() * simulation::material_xs.total;
|
||||
|
||||
// void sample_element(Particle* p)
|
||||
// {
|
||||
// associate (mat => materials(p->material))
|
||||
// // Sample cumulative distribution function
|
||||
// cutoff = prn() * simulation::material_xs.total
|
||||
// Get pointers to elements, densities
|
||||
int* nuclide;
|
||||
double* density;
|
||||
int n;
|
||||
openmc_material_get_densities(p->material, &nuclide, &density, &n);
|
||||
int* element = material_element(p->material);
|
||||
|
||||
// i = 0
|
||||
// prob = 0.0
|
||||
// do while (prob < cutoff)
|
||||
// i = i + 1
|
||||
int i = 0;
|
||||
double prob = 0.0;
|
||||
int i_element;
|
||||
while (prob < cutoff) {
|
||||
// Check to make sure that a nuclide was sampled
|
||||
if (i >= n) {
|
||||
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
|
||||
// TODO: off-by-one
|
||||
i_element = element[i] - 1;
|
||||
double atom_density = 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 += sigma;
|
||||
++i;
|
||||
}
|
||||
|
||||
// // Increment probability to compare to cutoff
|
||||
// prob = prob + sigma
|
||||
// end do
|
||||
// end associate
|
||||
// }
|
||||
return i_element;
|
||||
}
|
||||
|
||||
Reaction* sample_fission(int i_nuclide, double E)
|
||||
{
|
||||
|
|
|
|||
|
|
@ -1,20 +0,0 @@
|
|||
module physics_common
|
||||
|
||||
use, intrinsic :: ISO_C_BINDING
|
||||
|
||||
use particle_header, only: Particle
|
||||
|
||||
implicit none
|
||||
|
||||
!===============================================================================
|
||||
! RUSSIAN_ROULETTE FROM C
|
||||
!===============================================================================
|
||||
|
||||
interface
|
||||
subroutine russian_roulette(p) bind(C)
|
||||
import Particle
|
||||
type(Particle), intent(inout) :: p
|
||||
end subroutine russian_roulette
|
||||
end interface
|
||||
|
||||
end module physics_common
|
||||
|
|
@ -23,7 +23,6 @@ module settings
|
|||
integer(C_INT), bind(C) :: n_log_bins ! number of bins for logarithmic grid
|
||||
|
||||
logical(C_BOOL), bind(C) :: photon_transport
|
||||
integer(C_INT), bind(C) :: electron_treatment
|
||||
|
||||
! ============================================================================
|
||||
! MULTI-GROUP CROSS SECTION RELATED VARIABLES
|
||||
|
|
|
|||
|
|
@ -16,7 +16,6 @@ module tracking
|
|||
use mgxs_interface
|
||||
use nuclide_header
|
||||
use particle_header
|
||||
use physics, only: collision
|
||||
use random_lcg, only: prn, prn_set_stream
|
||||
use settings
|
||||
use simulation_header
|
||||
|
|
@ -33,11 +32,15 @@ module tracking
|
|||
implicit none
|
||||
|
||||
interface
|
||||
subroutine collision(p) bind(C)
|
||||
import Particle
|
||||
type(Particle), intent(inout) :: p
|
||||
end subroutine
|
||||
|
||||
subroutine collision_mg(p) bind(C)
|
||||
import Particle, C_DOUBLE
|
||||
type(Particle), intent(inout) :: p
|
||||
end subroutine collision_mg
|
||||
|
||||
end interface
|
||||
|
||||
contains
|
||||
|
|
|
|||
|
|
@ -1,9 +1,5 @@
|
|||
#!/bin/bash
|
||||
set -ex
|
||||
|
||||
# Allow tests that require GUI as described at:
|
||||
# https://docs.travis-ci.com/user/gui-and-headless-browsers/#Using-xvfb-to-Run-Tests-That-Require-a-GUI
|
||||
sh -e /etc/init.d/xvfb start
|
||||
|
||||
# Download NNDC HDF5 data, ENDF/B-VII.1 distribution, multipole library
|
||||
source tools/ci/download-xs.sh
|
||||
|
|
|
|||
|
|
@ -9,11 +9,10 @@ if [[ $DAGMC = 'y' ]]; then
|
|||
./tools/ci/travis-install-dagmc.sh
|
||||
fi
|
||||
|
||||
# Upgrade pip before doing anything else
|
||||
# Upgrade pip, pytest, numpy before doing anything else
|
||||
pip install --upgrade pip
|
||||
|
||||
# pytest installed by default -- make sure we get latest
|
||||
pip install --upgrade pytest
|
||||
pip install --upgrade numpy
|
||||
|
||||
# Install mpi4py for MPI configurations
|
||||
if [[ $MPI == 'y' ]]; then
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue