crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Remove bremsstrahlung code on Fortran side

Browse source
This commit is contained in:
Paul Romano 2019-01-15 06:51:18 -06:00
parent 2eacd5f72c
commit 9f8ebf527b
13 changed files with 81 additions and 636 deletions
Download patch file Download diff file Expand all files Collapse all files

1
CMakeLists.txt
View file

@ -331,7 +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

17
include/openmc/math_functions.h
View file

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

2
include/openmc/settings.h
View file

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

26
src/constants.F90
View file

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

59
src/input_xml.F90
View file

@ -2035,19 +2035,16 @@ contains
type(SetChar) :: element_already_read
interface
subroutine material_init_bremsstrahlung(ptr) bind(C)
import C_PTR
type(C_PTR), value :: ptr
subroutine photon_from_hdf5(group) bind(C)
import HID_T
integer(HID_T), value :: group
end subroutine
subroutine set_log_ttb_e_grid() bind(C)
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)
@ -2103,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
@ -2133,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 material_init_bremsstrahlung(materials(i) % ptr)
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)
call set_log_ttb_e_grid()
end if
call read_ce_cross_sections_c()
! Set up logarithmic grid for nuclides
do i = 1, size(nuclides)

31
src/material.cpp
View file

@ -1,5 +1,6 @@
#include "openmc/material.h"
#include <algorithm> // for min, max
#include <cmath>
#include <string>
#include <sstream>
@ -13,6 +14,7 @@
#include "openmc/nuclide.h"
#include "openmc/photon.h"
#include "openmc/search.h"
#include "openmc/settings.h"
#include "openmc/xml_interface.h"
namespace openmc {
@ -189,11 +191,11 @@ void Material::init_bremsstrahlung()
// Integrate the PDF using cubic spline integration over the incident
// particle energy
if (n > 2) {
spline_c(n, &data::ttb_e_grid(i), &f(i), &z(i));
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_c(n, &data::ttb_e_grid(i), &f(i), &z(i),
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;
@ -266,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
//==============================================================================

213
src/material_header.F90
View file

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

40
src/math.F90
View file

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

14
src/math_functions.cpp
View file

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

24
src/photon.cpp
View file

@ -3,6 +3,7 @@
#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"
@ -772,9 +773,22 @@ std::pair<double, double> klein_nishina(double alpha)
// 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.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)
@ -782,4 +796,12 @@ 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

173
src/photon_header.F90
View file

@ -2,53 +2,13 @@ module photon_header
use, intrinsic :: ISO_C_BINDING
use algorithm, only: binary_search
use constants
use dict_header, only: DictCharInt
use hdf5_interface
use settings
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 PhotonInteraction
character(3) :: name ! atomic symbol, e.g. 'Zr'
integer :: Z ! atomic number
! Microscopic cross sections
real(8), allocatable :: energy(:)
! 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
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
@ -70,142 +30,19 @@ module photon_header
contains
subroutine photon_from_hdf5(this, group_id)
class(PhotonInteraction), intent(inout) :: this
integer(HID_T), intent(in) :: group_id
integer :: i
integer(HID_T) :: rgroup
integer(HID_T) :: dset_id
integer(HSIZE_T) :: dims(1), dims2(2)
integer :: n_energy
integer :: n_k
integer :: n_e
real(8) :: f
real(8) :: y
real(8), allocatable :: electron_energy(:)
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)
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)
end subroutine photon_from_hdf5
!===============================================================================
! 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)

116
src/photon_physics.F90
View file

@ -1,116 +0,0 @@
module photon_physics
use algorithm, only: binary_search
use constants
use particle_header
use photon_header, only: BremsstrahlungData, ttb_e_grid, ttb
use random_lcg, only: prn
use settings
contains
!===============================================================================
! 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

1
src/settings.F90
View file

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