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

Merge branch 'photon-new-ttb' into photon-new

Browse source
This commit is contained in:
amandalund 2018-06-26 11:59:09 -05:00
commit 49e23fa605
8 changed files with 354 additions and 290 deletions
Download patch file Download diff file Expand all files Collapse all files

64
openmc/data/photon.py
View file

@ -94,9 +94,8 @@ _COMPTON_PROFILES = {}
# Stopping powers are read from a pre-generated HDF5 file when they are first
# needed. The dictionary stores an array of energy values at which the other
# quantities are tabulated with the key 'energy' and for each element has the
# mass density, the mean excitation energy, and arrays containing the collision
# stopping powers, radiative stopping powers, and the density effect parameter
# stored on the key 'Z'.
# mean excitation energy and arrays containing the collision stopping powers
# and radiative stopping powers stored on the key 'Z'.
_STOPPING_POWERS = {}
# Scaled bremsstrahlung DCSs are read from a data file provided by Selzter and
@ -106,23 +105,6 @@ _STOPPING_POWERS = {}
# for each element are in a 2D array with shape (n, k) stored on the key 'Z'.
_BREMSSTRAHLUNG = {}
# Reduced screening radii for Z = 1-99 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).
_REDUCED_SCREENING_RADIUS = [
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
]
class AtomicRelaxation(EqualityMixin):
"""Atomic relaxation data.
@ -384,17 +366,11 @@ class IncidentPhoton(EqualityMixin):
reactions : collections.OrderedDict
Contains the cross sections for each photon reaction. The keys are MT
values and the values are instances of :class:`PhotonReaction`.
reduced_screening_radius : float
Reduced screening radius :math:`R m_e c/\hbar`, where R is the screening
radius for an atom of atomic number Z under the assumption that the
Coulomb field of the nucleus is exponentially screened by atomic electrons.
:math:`\hbar/m_e c` is the Compton wavelength of the electron.
stopping_powers : dict
Dictionary of stopping power data with keys 'energy' (in eV), 'density'
(mass density in g/cm:sup:`3`), 'I' (mean excitation energy),
's_collision' (collision stopping power in eV cm:sup:`2`/g),
's_radiative' (radiative stopping power in eV cm:sup:`2`/g), and
'density_effect' (density effect parameter).
Dictionary of stopping power data with keys 'energy' (in eV), 'I' (mean
excitation energy), 's_collision' (collision stopping power in
eV cm:sup:`2`/g), and 's_radiative' (radiative stopping power in
eV cm:sup:`2`/g)
summed_reactions : collections.OrderedDict
Contains summed cross sections. The keys are MT values and the values
are instances of :class:`PhotonReaction`.
@ -439,14 +415,6 @@ class IncidentPhoton(EqualityMixin):
def name(self):
return ATOMIC_SYMBOL[self.atomic_number]
@property
def reduced_screening_radius(self):
if self.atomic_number < 100:
return _REDUCED_SCREENING_RADIUS[self.atomic_number - 1]
else:
raise IndexError('No reduced screening radius for '
'Z={}.'.format(self.atomic_number))
@atomic_number.setter
def atomic_number(self, atomic_number):
cv.check_type('atomic number', atomic_number, Integral)
@ -620,11 +588,9 @@ class IncidentPhoton(EqualityMixin):
_STOPPING_POWERS['energy'] = f['energy'].value*EV_PER_MEV
for i in range(1, 99):
group = f['{:03}'.format(i)]
_STOPPING_POWERS[i] = {'density': group.attrs['density'],
'I': group.attrs['I'],
_STOPPING_POWERS[i] = {'I': group.attrs['I'],
's_collision': group['s_collision'].value,
's_radiative': group['s_radiative'].value,
'density_effect': group['density_effect'].value}
's_radiative': group['s_radiative'].value}
# Units are in MeV cm^2/g; convert to eV cm^2/g
_STOPPING_POWERS[i]['s_collision'] *= EV_PER_MEV
@ -665,15 +631,15 @@ class IncidentPhoton(EqualityMixin):
# Get the scaled cross section values for each electron energy and
# reduced photon energy for this Z
logy = np.log(np.reshape(np.fromiter(brem[p:p+n*k], float, n*k), (n, k)))
y = np.reshape(np.fromiter(brem[p:p+n*k], float, n*k), (n, k))
p += k*n
for j in range(k):
# Cubic spline log-log interpolation
cs = CubicSpline(logx, logy[:,j])
# Cubic spline interpolation in log energy and linear DCS
cs = CubicSpline(logx, y[:,j])
# Get scaled DCS values (millibarns) on new energy grid
dcs[:,j] = np.exp(cs(log_energy))
dcs[:,j] = cs(log_energy)
_BREMSSTRAHLUNG[i] = {'dcs': dcs}
@ -784,10 +750,6 @@ class IncidentPhoton(EqualityMixin):
shell_group.attrs['designators'] = np.array(designators, dtype='S')
# Write reduced screening radius
if Z < 100:
group.attrs['reduced_screening_radius'] = self.reduced_screening_radius
# Write Compton profiles
if self.compton_profiles:
compton_group = group.create_group('compton_profiles')
@ -810,7 +772,7 @@ class IncidentPhoton(EqualityMixin):
s_group = group.create_group('stopping_powers')
for key, value in self.stopping_powers.items():
if key in ('density', 'I'):
if key == 'I':
s_group.attrs[key] = value
else:
s_group.create_dataset(key, data=value)

8
scripts/openmc-make-stopping-powers
View file

@ -12,7 +12,7 @@ from openmc.data import ATOMIC_SYMBOL
base_url = 'https://physics.nist.gov/cgi-bin/Star/e_table-t.pl'
energies = np.logspace(-3, 3, 200)
data = {'matno': '', 'Energies': '\n'.join(str(x) for x in energies)}
columns = {1: 's_collision', 2: 's_radiative', 4: 'density_effect'}
columns = {1: 's_collision', 2: 's_radiative'}
# ==============================================================================
# SCRAPE DATA FROM ESTAR SITE AND GENERATE STOPPING POWER HDF5 FILE
@ -41,12 +41,10 @@ with h5py.File('stopping_powers.h5', 'w') as f:
# Create group for this element
group = f.create_group('{:03}'.format(Z))
# Write the density and mean excitation energy
# Write the mean excitation energy
attributes = np.fromstring(r[3], sep=' ')
group.attrs['density'] = attributes[1]
group.attrs['I'] = attributes[2]
# Write collision and radiative stopping powers and density effect
# parameter
# Write collision and radiative stopping powers
for i in columns:
group.create_dataset(columns[i], data=values[i])

15
src/input_xml.F90
View file

@ -4370,9 +4370,10 @@ contains
end if
end do
! Generate material bremsstrahlung data
! Generate material bremsstrahlung data for electrons and positrons
if (photon_transport .and. electron_treatment == ELECTRON_TTB) then
call bremsstrahlung_init(ttb(i), i)
call bremsstrahlung_init(ttb(i) % electron, i, ELECTRON)
call bremsstrahlung_init(ttb(i) % positron, i, POSITRON)
end if
end do
@ -4385,9 +4386,17 @@ contains
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
! Take logarithm of electron energies since they are log-log interpolated
! 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)))
end if
! Take logarithm of energies since they are log-log interpolated
ttb_e_grid = log(ttb_e_grid)
end if

284
src/material_header.F90
View file

@ -699,165 +699,221 @@ contains
end function openmc_material_set_densities
subroutine bremsstrahlung_init(this, i_material)
class(Bremsstrahlung), intent(inout) :: this
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_e, n_k
real(8) :: e
integer :: n, n_e, n_k
real(8) :: c
real(8) :: k, k_l, k_r, k_c
real(8) :: x_l, x_r, x_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) :: density
real(8) :: density_gpcc
real(8) :: Z_eq_sq
real(8) :: beta
real(8) :: atom_sum
real(8) :: mass_sum
real(8), allocatable :: atom_fraction(:)
real(8), allocatable :: mass_fraction(:)
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 :: mfp_inv(:)
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)
this % i_material = i_material
! Allocate and initialize arrays
n_k = size(ttb_k_grid)
n_e = size(ttb_e_grid)
allocate(atom_fraction(mat % n_nuclides))
allocate(mass_fraction(mat % n_nuclides))
allocate(stopping_power(n_e))
allocate(mfp_inv(n_e))
allocate(this % yield(n_e))
allocate(this % dcs(n_k, n_e))
allocate(this % cdf(n_k, n_e))
allocate(z(n_e))
stopping_power(:) = ZERO
mfp_inv(:) = ZERO
this % dcs(:,:) = ZERO
this % cdf(:,:) = ZERO
! Calculate the "equivalent" atomic number Zeq, the atomic fraction and the
! mass fraction of each element, and the material density in atom/b-cm and
! in g/cm^3
Z_eq_sq = ZERO
do i = 1, mat % n_nuclides
awr = nuclides(mat % nuclide(i)) % awr
! Given atom percent
if (mat % atom_density(1) > ZERO) then
atom_fraction(i) = mat % atom_density(i)
mass_fraction(i) = mat % atom_density(i) * awr
! Given weight percent
else
atom_fraction(i) = -mat % atom_density(i) / awr
mass_fraction(i) = -mat % atom_density(i)
end if
Z_eq_sq = Z_eq_sq + atom_fraction(i) * nuclides(mat % nuclide(i)) % Z**2
end do
atom_sum = sum(atom_fraction)
mass_sum = sum(mass_fraction)
! Given material density in g/cm^3
if (mat % density < ZERO) then
density = -mat % density * (atom_sum / mass_sum) * N_AVOGADRO / MASS_NEUTRON
density_gpcc = -mat % density
! Given material density in atom/b-cm
! Determine whether we are generating electron or positron data
if (particle == POSITRON) then
positron_ = .true.
else
density = mat % density
density_gpcc = mat % density * (mass_sum / atom_sum) * MASS_NEUTRON / &
N_AVOGADRO
positron_ = .false.
end if
Z_eq_sq = Z_eq_sq / atom_sum
atom_fraction = atom_fraction / atom_sum
mass_fraction = mass_fraction / mass_sum
! 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. Note: 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.
! 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))
! TODO: for molecular DCS, atom_fraction should actually be the number of
! atoms in the molecule.
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
this % dcs = this % dcs + atom_fraction(i) * elm % Z**2 / Z_eq_sq * elm % dcs
dcs = dcs + atom_density * elm % Z**2 * elm % dcs
! Accumulate material total stopping power
stopping_power = stopping_power + mass_fraction(i) * density_gpcc * &
(elm % stopping_power_collision + elm % stopping_power_radiative)
! 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 inverse bremsstrahlung mean free path
do i = 1, n_e
e = ttb_e_grid(i)
if (e <= energy_cutoff(PHOTON)) cycle
! 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))
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
! Ratio of the velocity of the charged particle to the speed of light
beta = sqrt(e*(e + TWO*MASS_ELECTRON)) / (e + MASS_ELECTRON)
! Total material stopping power
stopping_power = stopping_power_collision + stopping_power_radiative
! Integration lower bound
k_c = energy_cutoff(PHOTON) / e
! Loop over photon energies
do i = 1, n_e - 1
w = ttb_e_grid(i)
! Find the upper bounding index of the reduced photon cutoff energy
i_k = binary_search(ttb_k_grid, n_k, k_c) + 1
! Loop over incident particle energies
do j = i, n_e
e = ttb_e_grid(j)
! Get the interpolation bounds
k_l = ttb_k_grid(i_k-1)
k_r = ttb_k_grid(i_k)
x_l = this % dcs(i_k-1, i)
x_r = this % dcs(i_k, i)
! Reduced photon energy
k = w / e
! Use linear interpolation in reduced photon energy k to find value of
! the DCS at the cutoff energy
x_c = (x_l * (k_r - k_c) + x_r * (k_c - k_l)) / (k_r - k_l)
! Find the lower bounding index of the reduced photon energy
i_k = binary_search(ttb_k_grid, n_k, k)
! Calculate the CDF using the trapezoidal rule in log-log space
c = HALF * (log(k_r) - log(k_c)) * (x_c + x_r)
this % cdf(i_k,i) = c
do j = i_k, n_k - 1
c = c + HALF * (log(ttb_k_grid(j+1)) - log(ttb_k_grid(j))) * &
(this % dcs(j,i) + this % dcs(j+1,i))
this % cdf(j+1,i) = c
! 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)) / (e + MASS_ELECTRON)
! Compute the integrand of the PDF
f(j) = (1.0e-3_8 * x) / (beta**2 * stopping_power(j) * w)
end do
! Calculate the inverse bremsstrahlung mean free path
mfp_inv(i) = c * density * Z_eq_sq / beta**2 * 1.0e-3_8
! 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(ttb_e_grid(i:), f(i:), z(i:), n)
c = ZERO
do j = i, n_e - 1
c = c + spline_integrate(ttb_e_grid(i:), f(i:), z(i:), n, &
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
! Calculate photon number yield
mfp_inv(:) = mfp_inv(:) / stopping_power(:)
call spline(ttb_e_grid, mfp_inv, z, n_e)
do i = 1, n_e
this % yield(i) = spline_integrate(ttb_e_grid, mfp_inv, z, n_e, &
energy_cutoff(PHOTON), ttb_e_grid(i))
! 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
deallocate(atom_fraction, mass_fraction, stopping_power, mfp_inv, z)
deallocate(stopping_power_collision, stopping_power_radiative, &
stopping_power, dcs, f, z)
end subroutine bremsstrahlung_init

2
src/math.F90
View file

@ -925,7 +925,7 @@ contains
integer :: i
integer :: ia, ib
real(8) :: h, r
real(8) :: a, b, c, d
real(8) :: b, c, d
! Find the lower bounding index in x of the lower limit of integration.
if (xa < x(1)) then

97
src/photon_header.F90
View file

@ -12,8 +12,8 @@ module photon_header
use settings
real(8), allocatable :: compton_profile_pz(:)
real(8), allocatable :: ttb_e_grid(:) ! incident electron energy grid
real(8), allocatable :: ttb_k_grid(:) ! reduced photon energy grid
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
@ -53,11 +53,6 @@ module photon_header
! dictionary gives an index in shells(:)
type(ElectronSubshell), allocatable :: shells(:)
! Pair production data
real(8) :: reduced_screening_radius
real(8) :: coulomb_correction
real(8) :: correction_factor_coeffs(4)
! Compton profile data
real(8), allocatable :: profile_pdf(:,:)
real(8), allocatable :: profile_cdf(:,:)
@ -65,7 +60,7 @@ module photon_header
real(8), allocatable :: electron_pdf(:)
! Stopping power data
real(8) :: density
real(8) :: I ! mean excitation energy
real(8), allocatable :: stopping_power_collision(:)
real(8), allocatable :: stopping_power_radiative(:)
@ -77,12 +72,15 @@ module photon_header
procedure :: calculate_xs => photon_calculate_xs
end type PhotonInteraction
type Bremsstrahlung
integer :: i_material ! Index in materials array
real(8), allocatable :: yield(:) ! Photon number yield
real(8), allocatable :: dcs(:,:) ! Bremsstrahlung scaled DCS
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
@ -90,7 +88,7 @@ module photon_header
type(DictCharInt) :: element_dict
type(Bremsstrahlung), allocatable, target :: ttb(:) ! Bremsstrahlung cross sections
type(Bremsstrahlung), allocatable, target :: ttb(:) ! Bremsstrahlung data
!===============================================================================
! ELEMENTMICROXS contains cached microscopic photon cross sections for a
@ -129,9 +127,12 @@ contains
integer :: n_k
integer :: n_e
character(3), allocatable :: designators(:)
real(8) :: a
real(8) :: c
real(8) :: f
real(8) :: y
real(8), allocatable :: electron_energy(:)
real(8), allocatable :: matrix(:,:)
real(8), allocatable :: dcs(:,:)
! Get name of nuclide from group
name_len = len(this % name)
@ -281,22 +282,6 @@ contains
call read_dataset(this % binding_energy, rgroup, 'binding_energy')
this % electron_pdf(:) = this % electron_pdf / sum(this % electron_pdf)
! Get reduced screening radius
call read_attribute(this % reduced_screening_radius, group_id, &
'reduced_screening_radius')
! Compute the high-energy Coulomb correction
a = this % Z / FINE_STRUCTURE
this % coulomb_correction = 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)
! Compute the coefficients of the correction factor
this % correction_factor_coeffs(1) = -0.1774_8 - 12.10_8*a + 11.18_8*a**2
this % correction_factor_coeffs(2) = 8.523_8 + 73.26_8*a - 44.41_8*a**2
this % correction_factor_coeffs(3) = -13.52_8 - 121.1_8*a + 96.41_8*a**2
this % correction_factor_coeffs(4) = 8.946_8 + 62.05_8*a - 63.41_8*a**2
! Read Compton profiles
dset_id = open_dataset(rgroup, 'J')
call get_shape(dset_id, dims2)
@ -340,10 +325,8 @@ contains
call close_dataset(dset_id)
! Get energy grids used for bremsstrahlung DCS and for stopping powers
if (.not. allocated(ttb_e_grid)) then
allocate(ttb_e_grid(n_e))
call read_dataset(ttb_e_grid, rgroup, 'electron_energy')
end if
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')
@ -357,9 +340,50 @@ contains
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 % density, rgroup, 'density')
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
@ -481,7 +505,6 @@ contains
! Clear TTB-related arrays
if (allocated(ttb_e_grid)) deallocate(ttb_e_grid)
if (allocated(ttb_k_grid)) deallocate(ttb_k_grid)
if (allocated(ttb)) deallocate(ttb)
end subroutine free_memory_photon

173
src/photon_physics.F90
View file

@ -3,8 +3,8 @@ module photon_physics
use algorithm, only: binary_search
use constants
use particle_header, only: Particle
use photon_header, only: PhotonInteraction, Bremsstrahlung, &
compton_profile_pz, ttb_e_grid, ttb_k_grid, ttb
use photon_header, only: PhotonInteraction, BremsstrahlungData, &
compton_profile_pz, ttb_e_grid, ttb
use random_lcg, only: prn
use settings
@ -402,24 +402,35 @@ contains
integer :: i
real(8) :: f
real(8) :: c
real(8) :: a
real(8) :: b
real(8) :: r
real(8) :: q
real(8) :: rn
real(8) :: beta
real(8) :: mu
real(8) :: phi
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) :: c(4)
! 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
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
@ -427,23 +438,37 @@ contains
! 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.
r = elm % reduced_screening_radius
! 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.
a = sqrt(TWO/alpha)
c = elm % correction_factor_coeffs
f = c(1)*a + c(2)*a**2 + c(3)*a**3 + c(4)*a**4
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/alpha
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) - FOUR*elm % coulomb_correction + f
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
@ -478,7 +503,7 @@ contains
end if
! Calculate phi_i(e) and deliver e if rn <= U_i(e)
b = r/(TWO*alpha*e*(ONE - 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))
@ -520,32 +545,33 @@ contains
real(8), intent(inout) :: E_lost
integer :: i, j
integer :: i_e, i_k
integer :: i_e, i_w
integer :: n
integer :: n_e, n_k
real(8) :: c_max
integer :: n_e
real(8) :: a
real(8) :: f
real(8) :: w
real(8) :: r
real(8) :: e, e_l, e_r
real(8) :: y, y_l, y_r
real(8) :: k, k_l, k_r, k_c
real(8) :: x, x_l, x_r
type(Bremsstrahlung), pointer :: mat
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 % E < energy_cutoff(PHOTON)) return
! Get bremsstrahlung data for this material
mat => ttb(p % material)
k_c = energy_cutoff(PHOTON) / p % E
! 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)
n_k = size(ttb_k_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)
@ -556,65 +582,56 @@ contains
! 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)
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)
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
! Sample index of the tabulated PDF in the energy grid, j or j+1
if (prn() > f) then
i_e = j
else
i_e = j + 1
! 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)
! TODO: interpolate maximum value of the CDF
! 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
if (w > energy_cutoff(PHOTON)) then
! Create secondary photon
call p % create_secondary(p % coord(1) % uvw, w, PHOTON, run_ce=.true.)
E_lost = E_lost + w
end if
! Maximum value of the CDF
c_max = mat % cdf(n_k, i_e)
! Sample reduced photon energy from the tabulated PDFs
do
! Generate a random number r and determine the index i for which
! cdf(i) <= r*cdf,max <= cdf(i+1)
r = prn()
i_k = binary_search(mat % cdf(:, i_e), n_k, r*c_max)
! Get interpolation bounds
k_l = ttb_k_grid(i_k)
k_r = ttb_k_grid(i_k+1)
x_l = mat % dcs(i_k, i_e)
x_r = mat % dcs(i_k+1, i_e)
if (k_l < k_c) then
x_l = x_l + (k_c - k_l) * (x_r - x_l) / (k_r - k_l)
k_l = k_c
end if
! Sample the reduced photon energy k from the distribution 1/k on the
! interval (k(i), k(i+1))
k = k_l * (k_r / k_l)**r
! Get the interpolated DCS
x = x_l + (k - k_l) * (x_r - x_l) / (k_r - k_l)
! Determine whether to deliver k
if (prn() * max(x_l, x_r) < x) exit
end do
w = k * p % E
if (w < energy_cutoff(PHOTON)) cycle
! Create secondary photon
call p % create_secondary(p % coord(1) % uvw, w, PHOTON, run_ce=.true.)
E_lost = E_lost + w
end do
end subroutine thick_target_bremsstrahlung

1
src/physics.F90
View file

@ -369,7 +369,6 @@ contains
! TODO: create reaction types
if (electron_treatment == ELECTRON_TTB) then
! TODO: implement thick-target bremsstrahlung model for positrons
call thick_target_bremsstrahlung(p, E_lost)
end if