diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 4a75849241..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 @@ -106,23 +105,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,17 +366,11 @@ 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), - '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`. @@ -439,14 +415,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) @@ -620,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 @@ -665,15 +631,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} @@ -784,10 +750,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') @@ -810,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]) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index ad3efb305c..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 @@ -4385,9 +4386,17 @@ 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 + ! 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 bf49f81528..f1e716f9fe 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -699,165 +699,221 @@ 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 - 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) :: 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(:) - real(8), allocatable :: mfp_inv(:) + 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 - n_k = size(ttb_k_grid) - n_e = size(ttb_e_grid) - 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(z(n_e)) - stopping_power(:) = ZERO - mfp_inv(:) = ZERO - this % dcs(:,:) = ZERO - this % cdf(:,:) = 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 - ! 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 + ! Determine whether we are generating electron or positron data + if (particle == POSITRON) then + positron_ = .true. else - density = mat % density - density_gpcc = mat % density * (mass_sum / atom_sum) * MASS_NEUTRON / & - N_AVOGADRO + positron_ = .false. end if - Z_eq_sq = Z_eq_sq / atom_sum - atom_fraction = atom_fraction / atom_sum - mass_fraction = mass_fraction / mass_sum + ! Get the size of the energy grids + n_k = size(ttb_k_grid) + n_e = size(ttb_e_grid) + + ! 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(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), source=ZERO) + allocate(f(n_e)) + allocate(z(n_e)) + + Z_eq_sq = ZERO + 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)) - ! TODO: for molecular DCS, atom_fraction should actually be the number of - ! atoms in the molecule. + 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 - this % dcs = this % dcs + atom_fraction(i) * elm % Z**2 / Z_eq_sq * elm % dcs + dcs = dcs + atom_density * elm % Z**2 * 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_density & + * MASS_NEUTRON / N_AVOGADRO * elm % stopping_power_collision + + ! Accumulate material radiative stopping power + 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 inverse bremsstrahlung mean free path - do i = 1, n_e - e = ttb_e_grid(i) - if (e <= energy_cutoff(PHOTON)) cycle + ! 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 - ! Ratio of the velocity of the charged particle to the speed of light - beta = sqrt(e*(e + TWO*MASS_ELECTRON)) / (e + MASS_ELECTRON) + ! Total material stopping power + stopping_power = stopping_power_collision + stopping_power_radiative - ! Integration lower bound - k_c = energy_cutoff(PHOTON) / e + ! Loop over photon energies + do i = 1, n_e - 1 + w = ttb_e_grid(i) - ! Find the upper bounding index of the reduced photon cutoff energy - i_k = binary_search(ttb_k_grid, n_k, k_c) + 1 + ! Loop over incident particle energies + do j = i, n_e + e = ttb_e_grid(j) - ! 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) + ! Reduced photon energy + k = w / e - ! 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 lower bounding index of the reduced photon energy + i_k = binary_search(ttb_k_grid, n_k, k) - ! 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 + ! 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) + + ! 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) + + ! 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) = (1.0e-3_8 * 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) = exp(-500.0_8) + + ! 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 + + ! Set photon number yield + 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) + elsewhere + this % yield = -500.0_8 end where - deallocate(atom_fraction, mass_fraction, stopping_power, mfp_inv, z) + deallocate(stopping_power_collision, stopping_power_radiative, & + 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 1af03c6b62..c180f95a83 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 @@ -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(:,:) @@ -65,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(:) @@ -77,12 +72,15 @@ module photon_header procedure :: calculate_xs => photon_calculate_xs end type PhotonInteraction - type Bremsstrahlung - integer :: i_material ! Index in materials array - - real(8), allocatable :: yield(:) ! Photon number yield - real(8), allocatable :: dcs(:,:) ! Bremsstrahlung scaled DCS + 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 +88,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 @@ -129,9 +127,12 @@ contains integer :: n_k integer :: n_e 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) @@ -281,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) @@ -340,10 +325,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') @@ -357,9 +340,50 @@ 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 + + ! 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 @@ -481,7 +505,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 60d1c09b42..cba7726cab 100644 --- a/src/photon_physics.F90 +++ b/src/photon_physics.F90 @@ -3,8 +3,8 @@ module photon_physics use algorithm, only: binary_search use constants use particle_header, only: Particle - use photon_header, only: PhotonInteraction, Bremsstrahlung, & - compton_profile_pz, ttb_e_grid, ttb_k_grid, ttb + use photon_header, only: PhotonInteraction, BremsstrahlungData, & + compton_profile_pz, ttb_e_grid, ttb use random_lcg, only: prn use settings @@ -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)) @@ -520,32 +545,33 @@ 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 - type(Bremsstrahlung), pointer :: mat + real(8) :: w, w_l, w_r + real(8) :: p_l, p_r + real(8) :: c, c_l, c_max + type(BremsstrahlungData), pointer :: mat if (p % E < energy_cutoff(PHOTON)) return - ! Get bremsstrahlung data for this material - mat => ttb(p % material) - - k_c = energy_cutoff(PHOTON) / p % E + ! 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) - 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) + if (j == n_e) j = j - 1 ! Get the interpolation bounds e_l = ttb_e_grid(j) @@ -556,65 +582,56 @@ 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 + + ! Sample index of the tabulated PDF in the energy grid, j or j+1 + 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-1, i_e) + p_r = mat % pdf(i_e, i_e) + 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 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 + ! 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 + + 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 - - ! 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 - call p % create_secondary(p % coord(1) % uvw, w, PHOTON, run_ce=.true.) - E_lost = E_lost + w end do end subroutine thick_target_bremsstrahlung 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