diff --git a/openmc/data/photon.py b/openmc/data/photon.py index c6202117c0..4a75849241 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 973bb0e4dc..ad3efb305c 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 @@ -4383,10 +4393,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 @@ -4429,9 +4439,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 70bf4bd21c..bf49f81528 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 162b4d526a..1af03c6b62 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 4b81e2598f..60d1c09b42 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