crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 13:45:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Generate Compton electrons; bug fixes

Browse source
This commit is contained in:
amandalund 2018-05-18 10:46:37 -05:00
parent 54f51dd033
commit 0c997d9492
10 changed files with 112 additions and 67 deletions
Download patch file Download diff file Expand all files Collapse all files

2
openmc/data/photon.py
View file

@ -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,

4
src/api.F90
View file

@ -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

2
src/cmfd_input.F90
View file

@ -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

28
src/input_xml.F90
View file

@ -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

2
src/material_header.F90
View file

@ -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)

8
src/nuclide_header.F90
View file

@ -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 &

30
src/photon_header.F90
View file

@ -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(&

42
src/photon_physics.F90
View file

@ -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

48
src/physics.F90
View file

@ -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

13
src/source_header.F90
View file

@ -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