From a0015bd74739d7306e65d08255e94fad40fedabf Mon Sep 17 00:00:00 2001 From: amandalund Date: Thu, 22 Mar 2018 15:37:57 -0500 Subject: [PATCH 01/10] In the process of updating bremsstrahlung CDF and PDF --- openmc/data/photon.py | 8 +-- src/input_xml.F90 | 3 +- src/material_header.F90 | 142 +++++++++++++++++++++++++--------------- src/math.F90 | 2 +- src/photon_header.F90 | 12 ++-- src/photon_physics.F90 | 87 +++++++++++------------- 6 files changed, 142 insertions(+), 112 deletions(-) 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 From 03e0d1145f06495dd8522e2fd86d1d4de3ea1009 Mon Sep 17 00:00:00 2001 From: amandalund Date: Thu, 29 Mar 2018 11:48:58 -0500 Subject: [PATCH 02/10] Debugging and dumping data --- src/material_header.F90 | 42 ++++++++++++++++++++++--- src/photon_physics.F90 | 68 ++++++++++++++++++++++++++--------------- 2 files changed, 81 insertions(+), 29 deletions(-) diff --git a/src/material_header.F90 b/src/material_header.F90 index 7a506abf3b..c6c6bbce84 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -872,7 +872,7 @@ contains 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) + this % pdf(j,j) = 1.0e-6_8 * this % pdf(j-1,j) ! Loop over photon energies c = ZERO @@ -888,13 +888,45 @@ contains 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 + ! Set photon number yield 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) + + ! Set small non-zero value at lowest energy + this % yield(1) = 1.0e-6_8 * this % yield(2) + + 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) + end where + deallocate(atom_fraction, mass_fraction, stopping_power, dcs, f, z) end subroutine bremsstrahlung_init diff --git a/src/photon_physics.F90 b/src/photon_physics.F90 index ee377ae488..85ed977e6e 100644 --- a/src/photon_physics.F90 +++ b/src/photon_physics.F90 @@ -399,7 +399,15 @@ contains real(8) :: c, c_l, c_max type(Bremsstrahlung), pointer :: mat - if (p % E < energy_cutoff(PHOTON)) return + real(8) :: photon_energies(100) + + !p % E = 100.0e6_8 + + !if (p % E < energy_cutoff(PHOTON)) return + if (p % E < energy_cutoff(PHOTON)) then + write(13,*) p % E, 0 + return + end if ! Get bremsstrahlung data for this material mat => ttb(p % material) @@ -409,6 +417,7 @@ contains ! 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) @@ -419,36 +428,38 @@ 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 + !if (n == 0) return + if (n == 0) then + write(13,*) p % E, n + return + end if ! 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 + 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, i_e-1) + p_l = mat % pdf(i_e-1, i_e) 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) + 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 @@ -464,17 +475,26 @@ write(*,*) "e_r: ", e_r, "e_l: ", e_l, "e_r/e_l: ", e_r/e_l 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) + 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 - - 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 + 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.) + E_lost = E_lost + w + 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 From b804181e60df301e21c51c3ad9006099a781a3ef Mon Sep 17 00:00:00 2001 From: amandalund Date: Thu, 26 Apr 2018 18:44:39 -0500 Subject: [PATCH 03/10] Accurately simulate pair production --- openmc/data/photon.py | 59 ++++++++++++++---- src/photon_header.F90 | 22 +++++++ src/photon_physics.F90 | 137 +++++++++++++++++++++++++++++++++++++++++ src/physics.F90 | 48 ++++++++------- 4 files changed, 230 insertions(+), 36 deletions(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 5ae81f25da..d6937a5990 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -106,6 +106,22 @@ _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. @@ -351,19 +367,6 @@ class IncidentPhoton(EqualityMixin): Number of protons in the target nucleus atomic_relaxation : openmc.data.AtomicRelaxation or None Atomic relaxation data - compton_profiles : dict - Dictionary of Compton profile data with keys 'num_electrons' (number of - electrons in each subshell), 'binding_energy' (ionization potential of - each subshell), and 'J' (Hartree-Fock Compton profile as a function of - the projection of the electron momentum on the scattering vector, - :math:`p_z` for each subshell). Note that subshell occupancies may not - match the atomic relaxation data. - 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). bremsstrahlung : dict Dictionary of bremsstrahlung DCS data with keys 'electron_energy' (incident electron kinetic energy values in eV), 'photon_energy' @@ -371,9 +374,27 @@ class IncidentPhoton(EqualityMixin): kinetic energy), and 'dcs' (cross sectin values in mb). The cross sections are in scaled form: :math:`(\beta^2/Z^2) E_k (d\sigma/dE_k)`, where :math:`E_k` is the energy of the emitted photon. + compton_profiles : dict + Dictionary of Compton profile data with keys 'num_electrons' (number of + electrons in each subshell), 'binding_energy' (ionization potential of + each subshell), and 'J' (Hartree-Fock Compton profile as a function of + the projection of the electron momentum on the scattering vector, + :math:`p_z` for each subshell). Note that subshell occupancies may not + match the atomic relaxation data. 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). summed_reactions : collections.OrderedDict Contains summed cross sections. The keys are MT values and the values are instances of :class:`PhotonReaction`. @@ -418,6 +439,14 @@ 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) @@ -755,6 +784,10 @@ 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') diff --git a/src/photon_header.F90 b/src/photon_header.F90 index d6ba13abec..557263bbb7 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -53,6 +53,11 @@ 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(:,:) @@ -124,6 +129,7 @@ contains integer :: n_k integer :: n_e character(3), allocatable :: designators(:) + real(8) :: a real(8) :: c real(8), allocatable :: matrix(:,:) @@ -274,6 +280,22 @@ 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) diff --git a/src/photon_physics.F90 b/src/photon_physics.F90 index 85ed977e6e..c181196167 100644 --- a/src/photon_physics.F90 +++ b/src/photon_physics.F90 @@ -378,6 +378,143 @@ contains end subroutine atomic_relaxation +!=============================================================================== +! PAIR_PRODUCTION samples the kinetic energy and direction of the electron and +! positron created when a photon is absorbed near an atomic nucleus. The +! simulation procedure follows the semiempirical model outlined in 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). +!=============================================================================== + + subroutine pair_production(elm, alpha, E_electron, E_positron, uvw_electron, & + uvw_positron) + type(PhotonInteraction), intent(in) :: elm + real(8), intent(in) :: alpha + real(8), intent(out) :: E_electron + real(8), intent(out) :: E_positron + real(8), intent(out) :: uvw_electron(3) + real(8), intent(out) :: uvw_positron(3) + + integer :: i + real(8) :: f + real(8) :: a + real(8) :: b + real(8) :: r + 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 + + ! The reduced screening radius r is the ratio of the screening radius to + ! the Compton wavelength of the electron, where the screening radius is + ! obtained under the assumption that the Coulomb field of the nucleus is + ! 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 + + ! 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 + + ! 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 + 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 + 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 + + ! To aid sampling, the unnormalized PDF can be expressed as + ! p = u_1*U_1(e)*pi_1(e) + u_2*U_2(e)*pi_2(e), where pi_1 and pi_2 are + ! normalized PDFs on the interval (e_min, e_max) from which values of e can + ! be sampled using the inverse transform method, and + ! U_1 = phi_1(e)/phi_1(1/2) and U_2 = phi_2(e)/phi_2(1/2) are valid + ! rejection functions. The reduced energy can now be sampled using a + ! combination of the composition and rejection methods. + u1 = TWO/THREE*(HALF - ONE/alpha)**2*phi1_max + u2 = phi2_max + do + rn = prn() + + ! Sample the index i in (1, 2) using the point probabilities + ! p(1) = u_1/(u_1 + u_2) and p(2) = u_2/(u_1 + u_2) + if (prn() < u1/(u1 + u2)) then + i = 1 + + ! Sample e from pi_1 using the inverse transform method + if (rn >= HALF) then + e = HALF + (HALF - ONE/alpha)*(TWO*rn - ONE)**(ONE/THREE) + else + e = HALF - (HALF - ONE/alpha)*(ONE - TWO*rn)**(ONE/THREE) + end if + else + i = 2 + + ! Sample e from pi_2 using the inverse transform method + e = ONE/alpha + (HALF - ONE/alpha)*TWO*rn + end if + + ! Calculate phi_i(e) and deliver e if rn <= U_i(e) + b = r/(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)) + if (i == 1) then + phi1 = 7.0_8/THREE - t1 - 6.0_8*t2 - t3 + t4 + if (prn() <= phi1/phi1_max) exit + else + phi2 = 11.0_8/6.0_8 - t1 - THREE*t2 + HALF*t3 + t4 + if (prn() <= phi2/phi2_max) exit + end if + end do + + ! Compute the kinetic energy of the electron and the positron + E_electron = (alpha*e - ONE)*MASS_ELECTRON + E_positron = (alpha*(ONE - e) - ONE)*MASS_ELECTRON + + ! Sample the direction of the electron. The cosine of the polar angle of + ! the direction relative to the incident photon is sampled from + ! p(mu) = C/(1 - beta*mu)^2 using the inverse transform method. + beta = sqrt(E_electron*(E_electron + TWO*MASS_ELECTRON)) & + / (E_electron + MASS_ELECTRON) + rn = TWO*prn() - ONE + mu = (rn + beta)/(rn*beta + ONE) + phi = TWO*PI*prn() + uvw_electron(1) = mu + uvw_electron(2) = sqrt(ONE - mu*mu)*cos(phi) + uvw_electron(3) = sqrt(ONE - mu*mu)*sin(phi) + + ! Sample the direction of the positron + beta = sqrt(E_positron*(E_positron + TWO*MASS_ELECTRON)) & + / (E_positron + MASS_ELECTRON) + rn = TWO*prn() - ONE + mu = (rn + beta)/(rn*beta + ONE) + phi = TWO*PI*prn() + uvw_positron(1) = mu + uvw_positron(2) = sqrt(ONE - mu*mu)*cos(phi) + uvw_positron(3) = sqrt(ONE - mu*mu)*sin(phi) + + end subroutine pair_production + !=============================================================================== ! THICK_TARGET_BREMSSTRAHLUNG !=============================================================================== diff --git a/src/physics.F90 b/src/physics.F90 index 9825848442..204158f059 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -12,7 +12,7 @@ module physics use particle_header, only: Particle use photon_header use photon_physics, only: rayleigh_scatter, compton_scatter, & - atomic_relaxation, & + atomic_relaxation, pair_production, & thick_target_bremsstrahlung use physics_common use random_lcg, only: prn, advance_prn_seed, prn_set_stream @@ -172,10 +172,15 @@ contains real(8) :: alpha ! photon energy divided by electron rest mass real(8) :: alpha_out ! outgoing photon energy over electron rest mass real(8) :: mu ! scattering cosine + real(8) :: mu_electron ! electron scattering cosine + real(8) :: mu_positron ! positron scattering cosine real(8) :: phi ! azimuthal angle - real(8) :: E_electron ! electron energy real(8) :: uvw(3) ! new direction real(8) :: rel_vel ! relative velocity of electron + real(8) :: E_electron ! electron energy + real(8) :: E_positron ! positron energy + real(8) :: uvw_electron(3) ! new electron direction + real(8) :: uvw_positron(3) ! new positron direction ! Kill photon if below energy cutoff -- an extra check is made here because ! photons with energy below the cutoff may have been produced by neutrons @@ -273,29 +278,26 @@ contains end do end if prob = prob_after + + ! Pair production + prob = prob + micro_photon_xs(i_element) % pair_production + if (prob > cutoff) then + + call pair_production(elm, alpha, E_electron, E_positron, uvw_electron, & + uvw_positron) + + ! Create secondary electron + call p % create_secondary(uvw_electron, E_electron, ELECTRON, .true.) + + ! Create secondary positron + call p % create_secondary(uvw_positron, E_positron, POSITRON, .true.) + + p % event_MT = PAIR_PROD + p % alive = .false. + p % E = ZERO + end if end associate - ! Pair production - prob = prob + micro_photon_xs(i_element) % pair_production - if (prob > cutoff) then - ! Sample angle isotropically - mu = TWO*prn() - ONE - phi = TWO*PI*prn() - uvw(1) = mu - uvw(2) = sqrt(ONE - mu*mu)*cos(phi) - uvw(3) = sqrt(ONE - mu*mu)*sin(phi) - - ! Compute the kinetic energy of each particle - E_electron = HALF * (p % E - 2 * MASS_ELECTRON) - - ! Create electron-positron pair traveling in opposite directions - call p % create_secondary( uvw, E_electron, ELECTRON, .true.) - call p % create_secondary(-uvw, E_electron, POSITRON, .true.) - p % event_MT = PAIR_PROD - p % alive = .false. - p % E = ZERO - end if - end subroutine sample_photon_reaction !=============================================================================== From 54f51dd03320a28ab7210d183ff49790bd4b05ae Mon Sep 17 00:00:00 2001 From: amandalund Date: Fri, 27 Apr 2018 14:23:33 -0500 Subject: [PATCH 04/10] Remove unnecessary lines --- src/material_header.F90 | 4 ---- src/photon_header.F90 | 1 - src/photon_physics.F90 | 6 +++++- 3 files changed, 5 insertions(+), 6 deletions(-) diff --git a/src/material_header.F90 b/src/material_header.F90 index c6c6bbce84..7fcdcd081b 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -911,10 +911,6 @@ contains end do close(17) - open(unit=14, file="yield.txt", action="write") - write(14,*) this % yield - close(14) - ! Set small non-zero value at lowest energy this % yield(1) = 1.0e-6_8 * this % yield(2) diff --git a/src/photon_header.F90 b/src/photon_header.F90 index 557263bbb7..1e56759f9a 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -469,7 +469,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 diff --git a/src/photon_physics.F90 b/src/photon_physics.F90 index c181196167..3f892b227e 100644 --- a/src/photon_physics.F90 +++ b/src/photon_physics.F90 @@ -4,7 +4,7 @@ module photon_physics use constants use particle_header, only: Particle use photon_header, only: PhotonInteraction, Bremsstrahlung, & - compton_profile_pz, ttb_e_grid, ttb_k_grid, ttb + compton_profile_pz, ttb_e_grid, ttb use random_lcg, only: prn use settings @@ -542,7 +542,9 @@ contains !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 @@ -577,7 +579,9 @@ 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 From 0c997d9492825d2b747e39d5d87b22b4bec65dd3 Mon Sep 17 00:00:00 2001 From: amandalund Date: Fri, 18 May 2018 10:46:37 -0500 Subject: [PATCH 05/10] Generate Compton electrons; bug fixes --- openmc/data/photon.py | 2 +- src/api.F90 | 4 ++-- src/cmfd_input.F90 | 2 +- src/input_xml.F90 | 28 ++++++++++++++++-------- src/material_header.F90 | 2 +- src/nuclide_header.F90 | 8 +++---- src/photon_header.F90 | 30 ++++++++++++++++++-------- src/photon_physics.F90 | 42 ++++++++++++++++-------------------- src/physics.F90 | 48 +++++++++++++++++++++++++++++++---------- src/source_header.F90 | 13 +++++------ 10 files changed, 112 insertions(+), 67 deletions(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index d6937a5990..ed12355eef 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -599,7 +599,7 @@ class IncidentPhoton(EqualityMixin): for i in range(1, 101): group = f['{:03}'.format(i)] num_electrons = group['num_electrons'].value - binding_energy = group['binding_energy'].value + binding_energy = group['binding_energy'].value*EV_PER_MEV J = group['J'].value _COMPTON_PROFILES[i] = {'num_electrons': num_electrons, 'binding_energy': binding_energy, diff --git a/src/api.F90 b/src/api.F90 index 8f33e85e32..6acb2497ce 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -118,8 +118,8 @@ contains create_fission_neutrons = .true. electron_treatment = ELECTRON_LED energy_cutoff(:) = [ZERO, 1000.0_8, ZERO, ZERO] - energy_max_neutron = INFINITY - energy_min_neutron = ZERO + energy_max(:) = [INFINITY, INFINITY] + energy_min(:) = [ZERO, ZERO] entropy_on = .false. gen_per_batch = 1 index_entropy_mesh = -1 diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index dbfabb2540..aca9194ca2 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -116,7 +116,7 @@ contains end if else if(.not.allocated(cmfd % egrid)) allocate(cmfd % egrid(2)) - cmfd % egrid = [ ZERO, energy_max_neutron ] + cmfd % egrid = [ ZERO, energy_max(NEUTRON) ] cmfd % indices(4) = 1 ! one energy group end if diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 4defed19a2..fa077549c3 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -4073,8 +4073,8 @@ contains end do ! Get the minimum and maximum energies - energy_min_neutron = energy_bins(num_energy_groups + 1) - energy_max_neutron = energy_bins(1) + energy_min(NEUTRON) = energy_bins(num_energy_groups + 1) + energy_max(NEUTRON) = energy_bins(1) ! Get the datasets present in the library call get_groups(file_id, names) @@ -4316,9 +4316,9 @@ contains ! Determine if minimum/maximum energy for this nuclide is greater/less ! than the previous if (size(nuclides(i_nuclide) % grid) >= 1) then - energy_min_neutron = max(energy_min_neutron, & + energy_min(NEUTRON) = max(energy_min(NEUTRON), & nuclides(i_nuclide) % grid(1) % energy(1)) - energy_max_neutron = min(energy_max_neutron, nuclides(i_nuclide) % & + energy_max(NEUTRON) = min(energy_max(NEUTRON), nuclides(i_nuclide) % & grid(1) % energy(size(nuclides(i_nuclide) % grid(1) % energy))) end if @@ -4345,6 +4345,16 @@ contains 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)))) + end if + ! Add element to set call element_already_read % add(element) end if @@ -4384,10 +4394,10 @@ contains ! Set up logarithmic grid for nuclides do i = 1, size(nuclides) - call nuclides(i) % init_grid(energy_min_neutron, & - energy_max_neutron, n_log_bins) + call nuclides(i) % init_grid(energy_min(NEUTRON), & + energy_max(NEUTRON), n_log_bins) end do - log_spacing = log(energy_max_neutron/energy_min_neutron) / n_log_bins + log_spacing = log(energy_max(NEUTRON)/energy_min(NEUTRON)) / n_log_bins do i = 1, size(materials) ! Skip materials with no S(a,b) tables @@ -4430,9 +4440,9 @@ contains ! grid has not been allocated if (size(nuclides(i) % grid) > 0) then if (nuclides(i) % grid(1) % energy(size(nuclides(i) % grid(1) % energy)) & - == energy_max_neutron) then + == energy_max(NEUTRON)) then call write_message("Maximum neutron transport energy: " // & - trim(to_str(energy_max_neutron)) // " eV for " // & + trim(to_str(energy_max(NEUTRON))) // " eV for " // & trim(adjustl(nuclides(i) % name)), 7) exit end if diff --git a/src/material_header.F90 b/src/material_header.F90 index 7fcdcd081b..460103b65a 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -299,7 +299,7 @@ contains logical :: check_sab ! should we check for S(a,b) table? ! Find energy index on energy grid - i_grid = int(log(p % E/energy_min_neutron)/log_spacing) + i_grid = int(log(p % E/energy_min(NEUTRON))/log_spacing) ! Determine if this material has S(a,b) tables check_sab = (this % n_sab > 0) diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index ef1ef693c1..27ad4c48dc 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -208,8 +208,8 @@ module nuclide_header !$omp threadprivate(micro_xs, material_xs) ! Minimum/maximum energies - real(8) :: energy_min_neutron = ZERO - real(8) :: energy_max_neutron = INFINITY + real(8) :: energy_min(2) = [ZERO, ZERO] + real(8) :: energy_max(2) = [INFINITY, INFINITY] contains @@ -1702,8 +1702,8 @@ contains if (res_scat_on) call nuclides(n) % assign_0K_elastic_scattering() ! Initialize nuclide grid - call nuclides(n) % init_grid(energy_min_neutron, & - energy_max_neutron, n_log_bins) + call nuclides(n) % init_grid(energy_min(NEUTRON), & + energy_max(NEUTRON), n_log_bins) else err = E_DATA call set_errmsg("Nuclide '" // trim(name_) // "' is not present & diff --git a/src/photon_header.F90 b/src/photon_header.F90 index 1e56759f9a..35b506a7f7 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -253,7 +253,8 @@ contains this % shells(i) % transition_subshells(:,:) = int(matrix(1:2, :), 4) this % shells(i) % transition_energy(:) = matrix(3, :) - this % shells(i) % transition_probability(:) = matrix(4, :) + this % shells(i) % transition_probability(:) = matrix(4, :) & + / sum(matrix(4, :)) deallocate(matrix) end if call close_dataset(dset_id) @@ -398,14 +399,16 @@ contains !=============================================================================== subroutine photon_calculate_xs(this, E, xs) - class(PhotonInteraction), intent(in) :: this ! index into nuclides array + class(PhotonInteraction), intent(in) :: this ! index into elements array real(8), intent(in) :: E ! energy type(ElementMicroXS), intent(inout) :: xs - integer :: i_grid ! index on nuclide energy grid - integer :: n_grid ! number of grid points - real(8) :: f ! interp factor on nuclide energy grid - real(8) :: log_E ! logarithm of the energy + integer :: i_grid ! index on element energy grid + integer :: i_shell ! index in subshells + integer :: i_start ! threshold index + integer :: n_grid ! number of grid points + real(8) :: f ! interp factor on element energy grid + real(8) :: log_E ! logarithm of the energy ! Perform binary search on the element energy grid in order to determine ! which points to interpolate between @@ -438,9 +441,18 @@ contains f*(this % incoherent(i_grid+1) - this % incoherent(i_grid))) ! Calculate microscopic photoelectric cross section - xs % photoelectric = exp(this % photoelectric_total(& - i_grid) + f*(this % photoelectric_total(i_grid+1) - & - this % photoelectric_total(i_grid))) + xs % photoelectric = ZERO + do i_shell = 1, size(this % shells) + ! Check threshold of reaction + i_start = this % shells(i_shell) % threshold + if (i_grid <= i_start) cycle + + ! Evaluation subshell photoionization cross section + xs % photoelectric = xs % photoelectric + & + exp(this % shells(i_shell) % cross_section(i_grid-i_start) + & + f*(this % shells(i_shell) % cross_section(i_grid+1-i_start) - & + this % shells(i_shell) % cross_section(i_grid-i_start))) + end do ! Calculate microscopic pair production cross section xs % pair_production = exp(& diff --git a/src/photon_physics.F90 b/src/photon_physics.F90 index 3f892b227e..88c68dd186 100644 --- a/src/photon_physics.F90 +++ b/src/photon_physics.F90 @@ -79,11 +79,12 @@ contains ! COMPTON_SCATTER !=============================================================================== - subroutine compton_scatter(el, alpha, alpha_out, mu, use_doppler) + subroutine compton_scatter(el, alpha, alpha_out, mu, i_shell, use_doppler) type(PhotonInteraction), intent(in) :: el real(8), intent(in) :: alpha real(8), intent(out) :: alpha_out real(8), intent(out) :: mu + integer, intent(out) :: i_shell logical, intent(in), optional :: use_doppler real(8) :: x @@ -118,8 +119,10 @@ contains ! Perform rejection on form factor if (prn() < form_factor_x / form_factor_xmax) then if (use_doppler_) then - call compton_doppler(el, alpha, mu, e_out) + call compton_doppler(el, alpha, mu, e_out, i_shell) alpha_out = e_out/MASS_ELECTRON + else + i_shell = 0 end if exit end if @@ -131,13 +134,14 @@ contains ! COMPTON_DOPPLER !=============================================================================== - subroutine compton_doppler(el, alpha, mu, e_out) + subroutine compton_doppler(el, alpha, mu, e_out, i_shell) type(PhotonInteraction), intent(in) :: el real(8), intent(in) :: alpha real(8), intent(in) :: mu real(8), intent(out) :: e_out + integer, intent(out) :: i_shell - integer :: i, i_shell + integer :: i integer :: n real(8) :: rn, m real(8) :: c, c_l, c_max @@ -335,9 +339,9 @@ contains ! Sample transition rn = prn() c = ZERO - do i_transition = 1, elm % shells(i_shell) % n_transitions - 1 + do i_transition = 1, elm % shells(i_shell) % n_transitions c = c + elm % shells(i_shell) % & - transition_probability(i_transition + 1) + transition_probability(i_transition) if (rn < c) exit end do @@ -387,14 +391,14 @@ contains ! Issy-les-Moulineaux, France (2011). !=============================================================================== - subroutine pair_production(elm, alpha, E_electron, E_positron, uvw_electron, & - uvw_positron) + subroutine pair_production(elm, alpha, E_electron, E_positron, mu_electron, & + mu_positron) type(PhotonInteraction), intent(in) :: elm real(8), intent(in) :: alpha real(8), intent(out) :: E_electron real(8), intent(out) :: E_positron - real(8), intent(out) :: uvw_electron(3) - real(8), intent(out) :: uvw_positron(3) + real(8), intent(out) :: mu_electron + real(8), intent(out) :: mu_positron integer :: i real(8) :: f @@ -491,27 +495,19 @@ contains E_electron = (alpha*e - ONE)*MASS_ELECTRON E_positron = (alpha*(ONE - e) - ONE)*MASS_ELECTRON - ! Sample the direction of the electron. The cosine of the polar angle of - ! the direction relative to the incident photon is sampled from + ! Sample the scattering angle of the electron. The cosine of the polar + ! angle of the direction relative to the incident photon is sampled from ! p(mu) = C/(1 - beta*mu)^2 using the inverse transform method. beta = sqrt(E_electron*(E_electron + TWO*MASS_ELECTRON)) & / (E_electron + MASS_ELECTRON) rn = TWO*prn() - ONE - mu = (rn + beta)/(rn*beta + ONE) - phi = TWO*PI*prn() - uvw_electron(1) = mu - uvw_electron(2) = sqrt(ONE - mu*mu)*cos(phi) - uvw_electron(3) = sqrt(ONE - mu*mu)*sin(phi) + mu_electron = (rn + beta)/(rn*beta + ONE) - ! Sample the direction of the positron + ! Sample the scattering angle of the positron beta = sqrt(E_positron*(E_positron + TWO*MASS_ELECTRON)) & / (E_positron + MASS_ELECTRON) rn = TWO*prn() - ONE - mu = (rn + beta)/(rn*beta + ONE) - phi = TWO*PI*prn() - uvw_positron(1) = mu - uvw_positron(2) = sqrt(ONE - mu*mu)*cos(phi) - uvw_positron(3) = sqrt(ONE - mu*mu)*sin(phi) + mu_positron = (rn + beta)/(rn*beta + ONE) end subroutine pair_production diff --git a/src/physics.F90 b/src/physics.F90 index 204158f059..624c8ed9b4 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -177,10 +177,9 @@ contains real(8) :: phi ! azimuthal angle real(8) :: uvw(3) ! new direction real(8) :: rel_vel ! relative velocity of electron + real(8) :: e_b ! binding energy of electron real(8) :: E_electron ! electron energy real(8) :: E_positron ! positron energy - real(8) :: uvw_electron(3) ! new electron direction - real(8) :: uvw_positron(3) ! new positron direction ! Kill photon if below energy cutoff -- an extra check is made here because ! photons with energy below the cutoff may have been produced by neutrons @@ -216,9 +215,34 @@ contains ! Incoherent (Compton) scattering prob = prob + micro_photon_xs(i_element) % incoherent if (prob > cutoff) then - call compton_scatter(elm, alpha, alpha_out, mu, .true.) + call compton_scatter(elm, alpha, alpha_out, mu, i_shell, .true.) + + ! Determine binding energy of shell. The binding energy is zero if + ! doppler broadening is not used. + if (i_shell == 0) then + e_b = ZERO + else + e_b = elm % binding_energy(i_shell) + end if + + ! Create Compton electron + E_electron = (alpha - alpha_out)*MASS_ELECTRON - e_b + mu_electron = (alpha - alpha_out*mu) & + / sqrt(alpha**2 + alpha_out**2 - TWO*alpha*alpha_out*mu) + phi = TWO*PI*prn() + uvw = rotate_angle(p % coord(1) % uvw, mu_electron, phi) + call p % create_secondary(uvw, E_electron, ELECTRON, .true.) + + ! TODO: Compton subshell data does not match atomic relaxation data + ! Allow electrons to fill orbital and produce auger electrons + ! and fluorescent photons + if (i_shell > 0) then + call atomic_relaxation(p, elm, i_shell) + end if + + phi = phi + PI p % E = alpha_out*MASS_ELECTRON - p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu) + p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu, phi) p % event_MT = INCOHERENT return end if @@ -282,20 +306,22 @@ contains ! Pair production prob = prob + micro_photon_xs(i_element) % pair_production if (prob > cutoff) then - - call pair_production(elm, alpha, E_electron, E_positron, uvw_electron, & - uvw_positron) + call pair_production(elm, alpha, E_electron, E_positron, mu_electron, & + mu_positron) ! Create secondary electron - call p % create_secondary(uvw_electron, E_electron, ELECTRON, .true.) + uvw = rotate_angle(p % coord(1) % uvw, mu_electron) + call p % create_secondary(uvw, E_electron, ELECTRON, .true.) ! Create secondary positron - call p % create_secondary(uvw_positron, E_positron, POSITRON, .true.) + uvw = rotate_angle(p % coord(1) % uvw, mu_positron) + call p % create_secondary(uvw, E_positron, POSITRON, .true.) p % event_MT = PAIR_PROD p % alive = .false. p % E = ZERO end if + end associate end subroutine sample_photon_reaction @@ -1578,7 +1604,7 @@ contains call rxn % products(1 + group) % sample(E_in, site % E, mu) ! resample if energy is greater than maximum neutron energy - if (site % E < energy_max_neutron) exit + if (site % E < energy_max(NEUTRON)) exit ! check for large number of resamples n_sample = n_sample + 1 @@ -1602,7 +1628,7 @@ contains call rxn % products(1) % sample(E_in, site % E, mu) ! resample if energy is greater than maximum neutron energy - if (site % E < energy_max_neutron) exit + if (site % E < energy_max(NEUTRON)) exit ! check for large number of resamples n_sample = n_sample + 1 diff --git a/src/source_header.F90 b/src/source_header.F90 index 1c739d1d5a..3cac26aa7c 100644 --- a/src/source_header.F90 +++ b/src/source_header.F90 @@ -9,7 +9,7 @@ module source_header use error use geometry, only: find_cell use material_header, only: materials - use nuclide_header, only: energy_min_neutron, energy_max_neutron + use nuclide_header, only: energy_min, energy_max use particle_header, only: Particle use settings, only: photon_transport use string, only: to_lower @@ -295,13 +295,13 @@ contains ! Sample angle site % uvw(:) = this % angle % sample() - ! Check for monoenergetic source above maximum neutron energy + ! Check for monoenergetic source above maximum particle energy select type (energy => this % energy) type is (Discrete) - if (any(energy % x > energy_max_neutron)) then + if (any(energy % x > energy_max(this % particle))) then call fatal_error("Source energy above range of energies of at least & &one cross section table") - else if (any(energy % x < energy_min_neutron)) then + else if (any(energy % x < energy_min(this % particle))) then call fatal_error("Source energy below range of energies of at least & &one cross section table") end if @@ -311,8 +311,9 @@ contains ! Sample energy spectrum site % E = this % energy % sample() - ! Resample if energy falls outside minimum or maximum neutron energy - if (site % E < energy_max_neutron .and. site % E > energy_min_neutron) exit + ! Resample if energy falls outside minimum or maximum particle energy + if (site % E < energy_max(this % particle) .and. & + site % E > energy_min(this % particle)) exit end do ! Set delayed group From d4237668d20b88e9cc213cfa14844cd9f1f7246f Mon Sep 17 00:00:00 2001 From: amandalund Date: Thu, 31 May 2018 20:52:54 -0500 Subject: [PATCH 06/10] Allow TTB implementation to handle any photon cutoff energy --- src/input_xml.F90 | 7 ++++++ src/material_header.F90 | 2 +- src/photon_header.F90 | 51 +++++++++++++++++++++++++++++++++++++---- src/photon_physics.F90 | 34 +++++++++++++-------------- 4 files changed, 72 insertions(+), 22 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index fa077549c3..4d29cdc811 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -4388,6 +4388,13 @@ contains 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))) + end if + ! 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 460103b65a..6cdff7397c 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -872,7 +872,7 @@ contains 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-6_8 * this % pdf(j-1,j) + this % pdf(j,j) = exp(-500.0_8) ! Loop over photon energies c = ZERO diff --git a/src/photon_header.F90 b/src/photon_header.F90 index 35b506a7f7..bb48c7f40a 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -131,7 +131,11 @@ contains 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) @@ -340,10 +344,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') @@ -360,6 +362,47 @@ contains 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 diff --git a/src/photon_physics.F90 b/src/photon_physics.F90 index 88c68dd186..1a2e12379c 100644 --- a/src/photon_physics.F90 +++ b/src/photon_physics.F90 @@ -536,13 +536,13 @@ contains !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 + 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 mat => ttb(p % material) @@ -573,13 +573,13 @@ contains n = int(y + prn()) 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 + 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 @@ -628,9 +628,9 @@ 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) + !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 From 7fc8c32295def1c83ee4f7b25e2b679ce7561b7d Mon Sep 17 00:00:00 2001 From: amandalund Date: Fri, 8 Jun 2018 09:56:59 -0500 Subject: [PATCH 07/10] Added TTB positron treatment --- src/input_xml.F90 | 5 +-- src/material_header.F90 | 73 +++++++++++++++++++++++++++++++---------- src/photon_header.F90 | 11 ++++--- src/photon_physics.F90 | 16 +++++---- 4 files changed, 76 insertions(+), 29 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 4d29cdc811..6370fc9fb6 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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 diff --git a/src/material_header.F90 b/src/material_header.F90 index 6cdff7397c..4afff94c8a 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -699,9 +699,10 @@ 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 @@ -711,6 +712,8 @@ contains 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 @@ -720,33 +723,44 @@ contains real(8) :: mass_sum real(8), allocatable :: atom_fraction(:) real(8), allocatable :: mass_fraction(:) + 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) - this % i_material = i_material - ! Allocate and initialize arrays + ! Determine whether we are generating electron or positron data + if (particle == POSITRON) then + positron_ = .true. + else + positron_ = .false. + end if + + ! Get the size of the energy grids 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 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(atom_fraction(mat % n_nuclides)) allocate(mass_fraction(mat % n_nuclides)) + 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)) + allocate(dcs(n_k, n_e), source=ZERO) allocate(f(n_e)) allocate(z(n_e)) - 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 @@ -768,6 +782,7 @@ contains 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) @@ -800,16 +815,39 @@ contains ! 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. ! Accumulate material 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) + ! Accumulate material collision stopping power + stopping_power_collision = stopping_power_collision + & + mass_fraction(i) * density_gpcc * elm % stopping_power_collision + + ! Accumulate material radiative stopping power + stopping_power_radiative = stopping_power_radiative + & + mass_fraction(i) * density_gpcc * elm % stopping_power_radiative end do + ! 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 + + ! 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) @@ -923,7 +961,8 @@ contains this % yield = log(this % yield) end where - deallocate(atom_fraction, mass_fraction, stopping_power, dcs, f, z) + deallocate(atom_fraction, mass_fraction, stopping_power_collision, & + stopping_power_radiative, stopping_power, dcs, f, z) end subroutine bremsstrahlung_init diff --git a/src/photon_header.F90 b/src/photon_header.F90 index bb48c7f40a..f3f2a3c8da 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -77,12 +77,15 @@ module photon_header procedure :: calculate_xs => photon_calculate_xs end type PhotonInteraction - type Bremsstrahlung - integer :: i_material ! Index in materials array - + 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 +93,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 diff --git a/src/photon_physics.F90 b/src/photon_physics.F90 index 1a2e12379c..29724bb636 100644 --- a/src/photon_physics.F90 +++ b/src/photon_physics.F90 @@ -3,7 +3,7 @@ module photon_physics use algorithm, only: binary_search use constants use particle_header, only: Particle - use photon_header, only: PhotonInteraction, Bremsstrahlung, & + use photon_header, only: PhotonInteraction, BremsstrahlungData, & compton_profile_pz, ttb_e_grid, ttb use random_lcg, only: prn use settings @@ -530,9 +530,9 @@ contains real(8) :: w, w_l, w_r real(8) :: p_l, p_r real(8) :: c, c_l, c_max - type(Bremsstrahlung), pointer :: mat + type(BremsstrahlungData), pointer :: mat - real(8) :: photon_energies(100) + !real(8) :: photon_energies(100) !p % E = 100.0e6_8 @@ -544,8 +544,12 @@ contains ! return !end if - ! Get bremsstrahlung data for this material - mat => ttb(p % material) + ! 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) @@ -620,7 +624,7 @@ contains w = exp(w_l)*(a*(c - c_l)/(exp(w_l)*p_l) + ONE)**(ONE/a) end if - photon_energies(i) = w + !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.) From 3fce470f781eb63c00ea814db8969447ec80857c Mon Sep 17 00:00:00 2001 From: amandalund Date: Mon, 11 Jun 2018 14:17:57 -0500 Subject: [PATCH 08/10] Simplify/clean up calculation of material TTB data --- src/material_header.F90 | 91 +++++++++++++++-------------------------- 1 file changed, 32 insertions(+), 59 deletions(-) diff --git a/src/material_header.F90 b/src/material_header.F90 index 4afff94c8a..eec95164fd 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -715,14 +715,11 @@ contains 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(:) @@ -753,8 +750,6 @@ contains allocate(this % yield(n_e)) ! Allocate temporary arrays - allocate(atom_fraction(mat % n_nuclides)) - allocate(mass_fraction(mat % n_nuclides)) allocate(stopping_power_collision(n_e), source=ZERO) allocate(stopping_power_radiative(n_e), source=ZERO) allocate(stopping_power(n_e)) @@ -762,45 +757,8 @@ contains allocate(f(n_e)) allocate(z(n_e)) - ! 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 - else - density = mat % density - density_gpcc = mat % density * (mass_sum / atom_sum) * MASS_NEUTRON / & - N_AVOGADRO - end if - - Z_eq_sq = Z_eq_sq / atom_sum - atom_fraction = atom_fraction / atom_sum - mass_fraction = mass_fraction / mass_sum + 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 @@ -815,17 +773,34 @@ contains ! 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_fraction(i) * elm % Z**2 / Z_eq_sq * elm % dcs + dcs = dcs + atom_density * elm % Z**2 * elm % dcs ! Accumulate material collision stopping power - stopping_power_collision = stopping_power_collision + & - mass_fraction(i) * density_gpcc * elm % stopping_power_collision + 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_fraction(i) * density_gpcc * elm % stopping_power_radiative + 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 @@ -876,8 +851,7 @@ contains 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) + f(j) = (1.0e-3_8 * x) / (beta**2 * stopping_power(j) * w) end do ! Number of points to integrate @@ -949,9 +923,6 @@ contains end do close(17) - ! Set small non-zero value at lowest energy - this % yield(1) = 1.0e-6_8 * this % yield(2) - open(unit=14, file="yield.txt", action="write") write(14,*) this % yield close(14) @@ -959,10 +930,12 @@ contains ! 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_collision, & - stopping_power_radiative, stopping_power, dcs, f, z) + deallocate(stopping_power_collision, stopping_power_radiative, & + stopping_power, dcs, f, z) end subroutine bremsstrahlung_init From 492a4e950fffcf607d8d8592903b686e59aa1332 Mon Sep 17 00:00:00 2001 From: amandalund Date: Mon, 25 Jun 2018 13:20:23 -0500 Subject: [PATCH 09/10] Clean up --- openmc/data/photon.py | 34 ---------------- src/material_header.F90 | 40 +++++-------------- src/photon_header.F90 | 22 ----------- src/photon_physics.F90 | 86 ++++++++++++++++++++--------------------- src/physics.F90 | 1 - 5 files changed, 51 insertions(+), 132 deletions(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index ed12355eef..2449b0e387 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -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') diff --git a/src/material_header.F90 b/src/material_header.F90 index eec95164fd..f1e716f9fe 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -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) diff --git a/src/photon_header.F90 b/src/photon_header.F90 index f3f2a3c8da..c180f95a83 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -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) diff --git a/src/photon_physics.F90 b/src/photon_physics.F90 index 29724bb636..cba7726cab 100644 --- a/src/photon_physics.F90 +++ b/src/photon_physics.F90 @@ -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 diff --git a/src/physics.F90 b/src/physics.F90 index 624c8ed9b4..02872ef03e 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -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 From 057cd9a92a10112cd1b35605eef3e90d21dd8c04 Mon Sep 17 00:00:00 2001 From: amandalund Date: Tue, 26 Jun 2018 11:43:28 -0500 Subject: [PATCH 10/10] Get rid of unnecessary stopping power related data --- openmc/data/photon.py | 22 +++++++++------------- scripts/openmc-make-stopping-powers | 8 +++----- 2 files changed, 12 insertions(+), 18 deletions(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 2449b0e387..a723b114f1 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -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 @@ -368,11 +367,10 @@ class IncidentPhoton(EqualityMixin): Contains the cross sections for each photon reaction. The keys are MT values and the values are instances of :class:`PhotonReaction`. 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`. @@ -590,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 @@ -776,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) diff --git a/scripts/openmc-make-stopping-powers b/scripts/openmc-make-stopping-powers index 8ee9602efa..76ac6686c8 100755 --- a/scripts/openmc-make-stopping-powers +++ b/scripts/openmc-make-stopping-powers @@ -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])