diff --git a/openmc/data/photon.py b/openmc/data/photon.py index d2750d9344..c918a4dd4c 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -448,26 +448,40 @@ class IncidentPhoton(EqualityMixin): filename = os.path.join(os.path.dirname(__file__), 'BREMX.DAT') brem = open(filename, 'r').read().split() + # Incident electron kinetic energy grid + _BREMSSTRAHLUNG['electron_energy'] = np.logspace(-3, 3, 200) + log_energy = np.log(_BREMSSTRAHLUNG['electron_energy']) + # Get number of tabulated electron and photon energy values n = int(brem[37]) k = int(brem[38]) # Index in data - j = 39 + p = 39 - # Get incident electron kinetic energy values - _BREMSSTRAHLUNG['electron_energy'] = np.fromiter(brem[j:j+n], float, n) - j += n + # Get log of incident electron kinetic energy values, used for cubic + # spline interpolation in log energy + logx = np.log(np.fromiter(brem[p:p+n], float, n)) + p += n # Get reduced photon energy values - _BREMSSTRAHLUNG['photon_energy'] = np.fromiter(brem[j:j+k], float, k) - j += k + _BREMSSTRAHLUNG['photon_energy'] = np.fromiter(brem[p:p+k], float, k) + p += k for i in range(1, 101): + dcs = np.empty([len(log_energy), k]) + # Get the scaled cross section values for each electron energy and # reduced photon energy for this Z - dcs = np.reshape(np.fromiter(brem[j:j+n*k], float, n*k), (n, k)) - j += k*n + logy = np.log(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]) + + # Get scaled DCS values (millibarns) on new energy grid + dcs[:,j] = np.exp(cs(log_energy)) _BREMSSTRAHLUNG[i] = {'dcs': dcs} diff --git a/src/global.F90 b/src/global.F90 index 30d824a6d2..44d092c294 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -15,7 +15,7 @@ module global use mesh_header, only: RegularMesh use mgxs_header, only: Mgxs, MgxsContainer use nuclide_header - use photon_header, only: PhotonInteraction, ElementMicroXS + use photon_header, only: PhotonInteraction, Bremsstrahlung, ElementMicroXS use plot_header, only: ObjectPlot use sab_header, only: SAlphaBeta use set_header, only: SetInt @@ -97,6 +97,7 @@ module global type(Nuclide), allocatable, target :: nuclides(:) ! Nuclide cross-sections type(PhotonInteraction), allocatable :: elements(:) ! Photon cross sections type(SAlphaBeta), allocatable, target :: sab_tables(:) ! S(a,b) tables + type(Bremsstrahlung), allocatable :: ttb(:) ! Bremsstrahlung cross sections integer :: n_sab_tables ! Number of S(a,b) thermal scattering tables @@ -498,6 +499,8 @@ contains if (allocated(sab_tables)) deallocate(sab_tables) if (allocated(micro_xs)) deallocate(micro_xs) + if (allocated(ttb)) deallocate(ttb) + ! Deallocate external source if (allocated(external_source)) deallocate(external_source) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 84c70b87b3..541af42eb4 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -5417,6 +5417,7 @@ contains allocate(nuclides(n_nuclides_total)) allocate(elements(n_elements)) allocate(sab_tables(n_sab_tables)) + if (electron_treatment == ELECTRON_TTB) allocate(ttb(n_materials)) ! Read cross sections do i = 1, size(materials) @@ -5491,6 +5492,11 @@ contains materials(i) % fissionable = .true. end if end do + + ! Generate material bremsstrahlung data + if (photon_transport .and. electron_treatment == ELECTRON_TTB) then + call ttb(i) % init(i) + end if end do ! Set up logarithmic grid for nuclides diff --git a/src/photon_header.F90 b/src/photon_header.F90 index a5cdd9791b..bb610b058f 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -2,12 +2,16 @@ module photon_header use hdf5, only: HID_T, HSIZE_T, SIZE_T + use algorithm, only: binary_search use constants, only: ZERO, HALF, SUBSHELLS use dict_header, only: DictIntInt use endf_header, only: Tabulated1D use hdf5_interface real(8), allocatable :: compton_profile_pz(:) + real(8), allocatable :: ttb_energy_electron(:) ! incident electron energy grid + real(8), allocatable :: ttb_energy_photon(:) ! reduced photon energy grid + real(8) :: ttb_energy_cutoff type ElectronSubshell integer :: index_subshell ! index in SUBSHELLS @@ -53,10 +57,29 @@ module photon_header real(8), allocatable :: binding_energy(:) real(8), allocatable :: electron_pdf(:) + ! Stopping power data + real(8) :: density + real(8), allocatable :: stopping_power_collision(:) + real(8), allocatable :: stopping_power_radiative(:) + + ! Bremsstrahlung scaled DCS + real(8), allocatable :: dcs(:,:) + contains procedure :: from_hdf5 => photon_from_hdf5 end type PhotonInteraction + type Bremsstrahlung + integer :: i_material ! Index in materials array + + real(8), allocatable :: yield(:,:) ! Photon number yield + real(8), allocatable :: dcs(:,:) ! Scaled bremsstrahlung DCS + real(8), allocatable :: cdf(:,:) ! Bremsstrahlung energy CDF + + contains + procedure :: init => bremsstrahlung_init + end type Bremsstrahlung + !=============================================================================== ! ELEMENTMICROXS contains cached microscopic photon cross sections for a ! particular element at the current energy @@ -88,6 +111,8 @@ contains integer :: n_shell integer :: n_profile integer :: n_transition + integer :: n_k + integer :: n_e character(3), allocatable :: designators(:) real(8) :: c real(8), allocatable :: matrix(:,:) @@ -266,6 +291,43 @@ contains this % pair_production_total(:) = this % pair_production_nuclear + & this % pair_production_electron + if (electron_treatment == ELECTRON_TTB) then + ! TODO: read in + ttb_energy_cutoff = 1.e-3 + + ! Read bremsstrahlung scaled DCS + rgroup = open_group(group_id, 'bremsstrahlung') + dset_id = open_dataset(rgroup, 'dcs') + call get_shape(dset_id, dims2) + n_k = int(dims2(1), 4) + n_e = int(dims2(2), 4) + allocate(this % dcs(n_k, n_e)) + call read_dataset(this % dcs, dset_id) + call close_dataset(dset_id) + + ! Get energy grids used for bremsstrahlung DCS and for stopping powers + if (.not. allocated(ttb_energy_electron)) then + allocate(ttb_energy_electron(n_e)) + call read_dataset(ttb_energy_electron, rgroup, 'electron_energy') + end if + if (.not. allocated(ttb_energy_photon)) then + allocate(ttb_energy_photon(n_k)) + call read_dataset(ttb_energy_photon, rgroup, 'photon_energy') + end if + call close_group(rgroup) + + ! Read stopping power data + if (this % Z < 99) then + rgroup = open_group(group_id, 'stopping_powers') + allocate(this % stopping_power_collision(n_e)) + 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 close_group(rgroup) + end if + end if + ! Take logarithm of energies and cross sections since they are log-log ! interpolated this % energy = log(this % energy) @@ -294,6 +356,123 @@ contains this % pair_production_total = -500.0_8 end where + if (electron_treatment == ELECTRON_TTB) then + ttb_energy_electron = log(ttb_energy_electron) + end if + end subroutine photon_from_hdf5 + subroutine bremsstrahlung_init(this, i_material) + class(Bremsstrahlung), intent(inout) :: this + integer, intent(in) :: i_material + + integer :: i, j + integer :: i_k + integer :: n_e, n_k + real(8) :: e + real(8) :: c + real(8) :: k, k_l, k_r, k_c + real(8) :: x_l, x_r, x_c + real(8) :: Z_eq_sq + real(8) :: beta + real(8), allocatable :: mass_fraction(:) + real(8), allocatable :: stopping_power(:) + real(8), allocatable :: mfp_inv(:) + type(Material), pointer :: mat + type(PhotonInteraction), pointer :: elm + + ! Get pointer to this material + mat => materials(i_material) + this % i_material = i_material + + ! Allocate and initialize arrays + n_k = size(ttb_energy_photon) + n_e = size(ttb_energy_electron) + allocate(mass_fraction(size(mat % element))) + 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)) + stopping_power(:) = ZERO + mfp_inv(:) = ZERO + this % dcs(:,:) = ZERO + this % cdf(:,:) = ZERO + + ! TODO + ! Calculate the "equivalent" atomic number Zeq and get the mass fraction of + ! each element + Z_eq_sq = 0 + do i = 1, mat % n_nuclides + ! Zeq**2 = (atomic fraction of x)*Zx**2 + (atomic fraction of y)*Zy**2 + ... + end do + + ! 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. + do i = 1, size(mat % element) + ! Get pointer to current element + elm => mat % element(i) + + ! TODO: n_atoms + ! Accumulate material DCS + this % dcs = this % dcs + n_atoms(i) * elm % Z**2 / Z_eq_sq * elm % dcs + + ! Accumulate material total stopping power + stopping_power = stopping_power + mass_fraction(i) * elm % density * & + (elm % stopping_power_collision + elm % stopping_power_radiative) + end do + + ! Calculate inverse bremsstrahlung mean free path + do i = 1, n_e + e = exp(ttb_energy_electron(i)) + if (e <= ttb_energy_cutoff) cycle + + ! Ratio of the velocity of the charged particle to the speed of light + beta = sqrt(e*(e + TWO*MASS_ELECTRON/1.e6)) / (e + MASS_ELECTRON/1.e6) + + ! Integration lower bound + k_c = ttb_energy_cutoff / e + + ! Find the upper bounding index of the reduced photon cutoff energy + i_k = binary_search(ttb_energy_photon, n_k, k_c) + 1 + + ! Get the interpolation bounds + k_l = ttb_energy_photon(i_k-1) + k_r = ttb_energy_photon(i_k) + x_l = this % dcs(i_k-1, i) + x_r = this % dcs(i_k, i) + + ! 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) + + ! Integrate 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_energy_photon(j+1)) - & + log(ttb_energy_photon(j))) * (this % dcs(j,i) + this % dcs(j+1,i)) + this % cdf(j+1,i) = c + + ! TODO: density in atom/cm^3 + mfp_inv(i) = c * mat % density * Z_eq_sq / beta**2 * 1.0e-27_8 + end do + + ! TODO: + ! Calculate photon number yield + ! cs = cubic_spline(exp(ttb_energy_electron), mfp_inv / stopping_power) + ! do i = 1, n_e + ! yield(i) = cs % integrate(ttb_energy_cutoff, exp(ttb_energy_electron(i))) + ! Use logarithm of number yield since it is log-log interpolated + ! yield = log(yield) + + end subroutine bremsstrahlung_init + end module photon_header