diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 8644f8e7cc..5ae81f25da 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -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} diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 973bb0e4dc..4defed19a2 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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 diff --git a/src/material_header.F90 b/src/material_header.F90 index 70bf4bd21c..7a506abf3b 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -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 diff --git a/src/math.F90 b/src/math.F90 index 5e4c8f81b1..127f134581 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -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 diff --git a/src/photon_header.F90 b/src/photon_header.F90 index 15203e4ff7..d6ba13abec 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -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 diff --git a/src/photon_physics.F90 b/src/photon_physics.F90 index aea6a65884..ee377ae488 100644 --- a/src/photon_physics.F90 +++ b/src/photon_physics.F90 @@ -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