crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 21:55:41 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

In the process of updating bremsstrahlung CDF and PDF

Browse source
This commit is contained in:
amandalund 2018-03-22 15:37:57 -05:00
parent e9bd7a203e
commit a0015bd747
6 changed files with 142 additions and 112 deletions
Download patch file Download diff file Expand all files Collapse all files

8
openmc/data/photon.py
View file

@ -636,15 +636,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}

3
src/input_xml.F90
View file

@ -4375,9 +4375,10 @@ 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
! Take logarithm of energies since they are log-log interpolated
ttb_e_grid = log(ttb_e_grid)
end if

142
src/material_header.F90
View file

@ -705,11 +705,12 @@ contains
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) :: awr
real(8) :: density
real(8) :: density_gpcc
@ -720,7 +721,8 @@ contains
real(8), allocatable :: atom_fraction(:)
real(8), allocatable :: mass_fraction(:)
real(8), allocatable :: stopping_power(:)
real(8), allocatable :: mfp_inv(:)
real(8), allocatable :: dcs(:,:)
real(8), allocatable :: f(:)
real(8), allocatable :: z(:)
type(Material), pointer :: mat
type(PhotonInteraction), pointer :: elm
@ -732,18 +734,19 @@ contains
! Allocate and initialize arrays
n_k = size(ttb_k_grid)
n_e = size(ttb_e_grid)
allocate(this % pdf(n_e, n_e))
allocate(this % cdf(n_e, n_e))
allocate(this % yield(n_e))
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(dcs(n_k, n_e))
allocate(f(n_e))
allocate(z(n_e))
stopping_power(:) = ZERO
mfp_inv(:) = ZERO
this % dcs(:,:) = ZERO
this % pdf(:,:) = ZERO
this % cdf(:,:) = ZERO
stopping_power(:) = ZERO
dcs(:,:) = 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
@ -800,64 +803,99 @@ contains
! TODO: for molecular DCS, atom_fraction should actually be the number of
! atoms in the molecule.
! Accumulate material DCS
this % dcs = this % dcs + atom_fraction(i) * elm % Z**2 / Z_eq_sq * elm % dcs
dcs = dcs + atom_fraction(i) * elm % Z**2 / Z_eq_sq * elm % dcs
! Accumulate material total stopping power
stopping_power = stopping_power + mass_fraction(i) * density_gpcc * &
(elm % stopping_power_collision + elm % stopping_power_radiative)
end do
! Calculate inverse bremsstrahlung mean free path
do i = 1, n_e
e = ttb_e_grid(i)
if (e <= energy_cutoff(PHOTON)) cycle
! Loop over photon energies
do i = 1, n_e - 1
w = ttb_e_grid(i)
! Ratio of the velocity of the charged particle to the speed of light
beta = sqrt(e*(e + TWO*MASS_ELECTRON)) / (e + MASS_ELECTRON)
! Loop over incident particle energies
do j = i, n_e
e = ttb_e_grid(j)
! Integration lower bound
k_c = energy_cutoff(PHOTON) / e
! Reduced photon energy
k = w / e
! Find the upper bounding index of the reduced photon cutoff energy
i_k = binary_search(ttb_k_grid, n_k, k_c) + 1
! 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-1)
k_r = ttb_k_grid(i_k)
x_l = this % dcs(i_k-1, i)
x_r = this % dcs(i_k, i)
! 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)
! 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 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)
! 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
! 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) = (density * 1.0e-3_8 * Z_eq_sq * 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) = 1.0e-9_8 * this % pdf(j-1,j)
! 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
! Use logarithm of number yield since it is log-log interpolated
if (c > ZERO) then
c = log(c)
end if
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)
end where
deallocate(atom_fraction, mass_fraction, stopping_power, mfp_inv, z)
deallocate(atom_fraction, mass_fraction, 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

12
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
@ -60,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(:)
@ -75,9 +75,9 @@ module photon_header
type Bremsstrahlung
integer :: i_material ! Index in materials array
real(8), allocatable :: yield(:) ! Photon number yield
real(8), allocatable :: dcs(:,:) ! Bremsstrahlung scaled DCS
real(8), allocatable :: pdf(:,:) ! Bremsstrahlung energy PDF
real(8), allocatable :: cdf(:,:) ! Bremsstrahlung energy CDF
real(8), allocatable :: yield(:) ! Photon number yield
end type Bremsstrahlung
type(PhotonInteraction), allocatable, target :: elements(:) ! Photon cross sections
@ -334,7 +334,7 @@ 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
end if

87
src/photon_physics.F90
View file

@ -387,17 +387,16 @@ 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
real(8) :: w, w_l, w_r
real(8) :: p_l, p_r
real(8) :: c, c_l, c_max
type(Bremsstrahlung), pointer :: mat
if (p % E < energy_cutoff(PHOTON)) return
@ -405,11 +404,8 @@ contains
! Get bremsstrahlung data for this material
mat => ttb(p % material)
k_c = energy_cutoff(PHOTON) / p % E
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)
@ -433,50 +429,45 @@ contains
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) then
i_e = j
! Maximum value of the CDF
c_max = mat % cdf(i_e, i_e)
else
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, i_e-1)
p_r = mat % pdf(i_e, i_e)
c_l = mat % cdf(i_e, i_e-1)
write(*,*) "p_r: ", p_r, "p_l: ", p_l, "p_r/p_l: ", p_r/p_l
write(*,*) "e_r: ", e_r, "e_l: ", e_l, "e_r/e_l: ", e_r/e_l
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)
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
end if
! 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)
! 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