mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-26 21:25:36 -04:00
Clean up
This commit is contained in:
amandalund
parent
3fce470f78
commit
492a4e950f
5 changed files with 51 additions and 132 deletions
|
|
@ -106,23 +106,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,11 +367,6 @@ 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),
|
||||
|
|
@ -439,14 +417,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)
|
||||
|
|
@ -784,10 +754,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')
|
||||
|
|
|
|||
|
|
@ -761,14 +761,15 @@ contains
|
|||
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))
|
||||
|
|
@ -904,29 +905,6 @@ contains
|
|||
this % yield(j) = c
|
||||
end do
|
||||
|
||||
open(unit=13, file="energies.txt", action="write", status="replace")
|
||||
close(13)
|
||||
|
||||
open(unit=15, file="e_grid.txt", action="write")
|
||||
write(15,*) ttb_e_grid
|
||||
close(15)
|
||||
|
||||
open(unit=16, file="pdf.txt", action="write")
|
||||
do i = 1, n_e
|
||||
write(16,*) this % pdf(:,i)
|
||||
end do
|
||||
close(16)
|
||||
|
||||
open(unit=17, file="cdf.txt", action="write")
|
||||
do i = 1, n_e
|
||||
write(17,*) this % cdf(:,i)
|
||||
end do
|
||||
close(17)
|
||||
|
||||
open(unit=14, file="yield.txt", action="write")
|
||||
write(14,*) this % yield
|
||||
close(14)
|
||||
|
||||
! Use logarithm of number yield since it is log-log interpolated
|
||||
where (this % yield > ZERO)
|
||||
this % yield = log(this % yield)
|
||||
|
|
|
|||
|
|
@ -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(:,:)
|
||||
|
|
@ -132,7 +127,6 @@ contains
|
|||
integer :: n_k
|
||||
integer :: n_e
|
||||
character(3), allocatable :: designators(:)
|
||||
real(8) :: a
|
||||
real(8) :: c
|
||||
real(8) :: f
|
||||
real(8) :: y
|
||||
|
|
@ -288,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)
|
||||
|
|
|
|||
|
|
@ -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))
|
||||
|
|
@ -532,17 +557,7 @@ contains
|
|||
real(8) :: c, c_l, c_max
|
||||
type(BremsstrahlungData), pointer :: mat
|
||||
|
||||
!real(8) :: photon_energies(100)
|
||||
|
||||
!p % E = 100.0e6_8
|
||||
|
||||
if (p % E < energy_cutoff(PHOTON)) return
|
||||
!if (p % E < energy_cutoff(PHOTON)) then
|
||||
! open(unit=13, file="energies.txt", action="write", position="append")
|
||||
! write(13,*) p % E, 0
|
||||
! close(13)
|
||||
! return
|
||||
!end if
|
||||
|
||||
! Get bremsstrahlung data for this material and particle type
|
||||
if (p % type == POSITRON) then
|
||||
|
|
@ -578,12 +593,6 @@ contains
|
|||
|
||||
E_lost = ZERO
|
||||
if (n == 0) return
|
||||
!if (n == 0) then
|
||||
! open(unit=13, file="energies.txt", action="write", position="append")
|
||||
! write(13,*) p % E, n
|
||||
! close(13)
|
||||
! return
|
||||
!end if
|
||||
|
||||
! Sample index of the tabulated PDF in the energy grid, j or j+1
|
||||
if (prn() <= f .or. j == 1) then
|
||||
|
|
@ -617,14 +626,7 @@ contains
|
|||
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
|
||||
! Temporary fix
|
||||
if (i_w == i_e - 1) then
|
||||
w = exp(w_l)
|
||||
else
|
||||
w = exp(w_l)*(a*(c - c_l)/(exp(w_l)*p_l) + ONE)**(ONE/a)
|
||||
end if
|
||||
|
||||
!photon_energies(i) = w
|
||||
if (w > energy_cutoff(PHOTON)) then
|
||||
! Create secondary photon
|
||||
call p % create_secondary(p % coord(1) % uvw, w, PHOTON, run_ce=.true.)
|
||||
|
|
@ -632,10 +634,6 @@ contains
|
|||
end if
|
||||
end do
|
||||
|
||||
!open(unit=13, file="energies.txt", action="write", position="append")
|
||||
!write(13,*) p % E, n, photon_energies(:n)
|
||||
!close(13)
|
||||
|
||||
end subroutine thick_target_bremsstrahlung
|
||||
|
||||
end module photon_physics
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue