mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 13:45:36 -04:00
Initial implementation of photon physics
This commit is contained in:
Paul Romano
parent
2622489b3b
commit
94dd50422f
13 changed files with 732 additions and 23 deletions
|
|
@ -342,6 +342,7 @@ set(LIBOPENMC_FORTRAN_SRC
|
|||
src/particle_restart.F90
|
||||
src/particle_restart_write.F90
|
||||
src/photon_header.F90
|
||||
src/photon_physics.F90
|
||||
src/physics_common.F90
|
||||
src/physics.F90
|
||||
src/physics_mg.F90
|
||||
|
|
|
|||
|
|
@ -128,7 +128,8 @@ class Source(object):
|
|||
"""
|
||||
element = ET.Element("source")
|
||||
element.set("strength", str(self.strength))
|
||||
element.set("particle", self.particle)
|
||||
if self.particle != 'neutron':
|
||||
element.set("particle", self.particle)
|
||||
if self.file is not None:
|
||||
element.set("file", self.file)
|
||||
if self.space is not None:
|
||||
|
|
|
|||
|
|
@ -16,6 +16,7 @@ module bank_header
|
|||
real(C_DOUBLE) :: uvw(3) ! diretional cosines
|
||||
real(C_DOUBLE) :: E ! energy / energy group if in MG mode.
|
||||
integer(C_INT) :: delayed_group ! delayed group
|
||||
integer(C_INT) :: particle ! particle type (neutron, photon, etc.)
|
||||
end type Bank
|
||||
|
||||
end module bank_header
|
||||
|
|
|
|||
|
|
@ -20,11 +20,26 @@ module cross_section
|
|||
contains
|
||||
|
||||
!===============================================================================
|
||||
! CALCULATE_XS determines the macroscopic cross sections for the material the
|
||||
! particle is currently traveling through.
|
||||
! CALCULATE_XS determines the macroscopic neutron and/or photon cross sections
|
||||
! for the material the particle is currently traveling through.
|
||||
!===============================================================================
|
||||
|
||||
subroutine calculate_xs(p)
|
||||
type(Particle), intent(inout) :: p
|
||||
|
||||
if (p % type == NEUTRON) then
|
||||
call calculate_neutron_xs(p)
|
||||
elseif (p % type == PHOTON) then
|
||||
call calculate_photon_xs(p)
|
||||
end if
|
||||
end subroutine calculate_xs
|
||||
|
||||
!===============================================================================
|
||||
! CALCULATE_NEUTRON_XS determines the macroscopic cross sections for the
|
||||
! material the particle is currently traveling through.
|
||||
!===============================================================================
|
||||
|
||||
subroutine calculate_neutron_xs(p)
|
||||
|
||||
type(Particle), intent(inout) :: p
|
||||
|
||||
|
|
@ -125,7 +140,7 @@ contains
|
|||
end do
|
||||
end associate
|
||||
|
||||
end subroutine calculate_xs
|
||||
end subroutine calculate_neutron_xs
|
||||
|
||||
!===============================================================================
|
||||
! CALCULATE_NUCLIDE_XS determines microscopic cross sections for a nuclide of a
|
||||
|
|
@ -557,6 +572,140 @@ contains
|
|||
|
||||
end subroutine calculate_urr_xs
|
||||
|
||||
!===============================================================================
|
||||
! CALCULATE_PHOTON_XS determines the macroscopic photon cross sections for the
|
||||
! material the particle is currently traveling through.
|
||||
!===============================================================================
|
||||
|
||||
subroutine calculate_photon_xs(p)
|
||||
type(Particle), intent(inout) :: p
|
||||
|
||||
integer :: i ! loop index over nuclides
|
||||
integer :: i_element ! index into elements array
|
||||
real(8) :: atom_density ! atom density of a nuclide
|
||||
|
||||
! Set all material macroscopic cross sections to zero
|
||||
material_xs % total = ZERO
|
||||
material_xs % coherent = ZERO
|
||||
material_xs % incoherent = ZERO
|
||||
material_xs % photoelectric = ZERO
|
||||
material_xs % pair_production = ZERO
|
||||
|
||||
! Exit subroutine if material is void
|
||||
if (p % material == MATERIAL_VOID) return
|
||||
|
||||
associate (mat => materials(p % material))
|
||||
! Add contribution from each nuclide in material
|
||||
do i = 1, mat % n_nuclides
|
||||
! ========================================================================
|
||||
! CALCULATE MICROSCOPIC CROSS SECTION
|
||||
|
||||
! Determine microscopic cross sections for this nuclide
|
||||
i_element = mat % element(i)
|
||||
|
||||
! Calculate microscopic cross section for this nuclide
|
||||
if (p % E /= micro_photon_xs(i_element) % last_E) then
|
||||
call calculate_element_xs(i_element, p % E)
|
||||
end if
|
||||
|
||||
! ========================================================================
|
||||
! ADD TO MACROSCOPIC CROSS SECTION
|
||||
|
||||
! Copy atom density of nuclide in material
|
||||
atom_density = mat % atom_density(i)
|
||||
|
||||
! Add contributions to material macroscopic total cross section
|
||||
material_xs % total = material_xs % total + &
|
||||
atom_density * micro_photon_xs(i_element) % total
|
||||
|
||||
! Add contributions to material macroscopic coherent cross section
|
||||
material_xs % coherent = material_xs % coherent + &
|
||||
atom_density * micro_photon_xs(i_element) % coherent
|
||||
|
||||
! Add contributions to material macroscopic incoherent cross section
|
||||
material_xs % incoherent = material_xs % incoherent + &
|
||||
atom_density * micro_photon_xs(i_element) % incoherent
|
||||
|
||||
! Add contributions to material macroscopic photoelectric cross section
|
||||
material_xs % photoelectric = material_xs % photoelectric + &
|
||||
atom_density * micro_photon_xs(i_element) % photoelectric
|
||||
|
||||
! Add contributions to material macroscopic pair production cross section
|
||||
material_xs % pair_production = material_xs % pair_production + &
|
||||
atom_density * micro_photon_xs(i_element) % pair_production
|
||||
end do
|
||||
end associate
|
||||
|
||||
end subroutine calculate_photon_xs
|
||||
|
||||
!===============================================================================
|
||||
! CALCULATE_ELEMENT_XS determines microscopic photon cross sections for an
|
||||
! element of a given index in the elements array at the energy of the given
|
||||
! particle
|
||||
!===============================================================================
|
||||
|
||||
subroutine calculate_element_xs(i_element, E)
|
||||
integer, intent(in) :: i_element ! index into nuclides array
|
||||
real(8), intent(in) :: E ! energy
|
||||
|
||||
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
|
||||
|
||||
associate (elm => elements(i_element))
|
||||
! Perform binary search on the element energy grid in order to determine
|
||||
! which points to interpolate between
|
||||
n_grid = size(elm % energy)
|
||||
log_E = log(E)
|
||||
if (log_E <= elm % energy(1)) then
|
||||
i_grid = 1
|
||||
elseif (log_E > elm % energy(n_grid)) then
|
||||
i_grid = n_grid - 1
|
||||
else
|
||||
i_grid = binary_search(elm % energy, n_grid, log_E)
|
||||
end if
|
||||
|
||||
! check for case where two energy points are the same
|
||||
if (elm % energy(i_grid) == elm % energy(i_grid+1)) i_grid = i_grid + 1
|
||||
|
||||
! calculate interpolation factor
|
||||
f = (log_E - elm%energy(i_grid))/(elm%energy(i_grid+1) - elm%energy(i_grid))
|
||||
|
||||
micro_photon_xs(i_element) % index_grid = i_grid
|
||||
micro_photon_xs(i_element) % interp_factor = f
|
||||
|
||||
! Calculate microscopic coherent cross section
|
||||
micro_photon_xs(i_element) % coherent = exp(elm % coherent(i_grid) + f * &
|
||||
(elm % coherent(i_grid+1) - elm % coherent(i_grid)))
|
||||
|
||||
! Calculate microscopic incoherent cross section
|
||||
micro_photon_xs(i_element) % incoherent = exp(elm % incoherent(i_grid) + &
|
||||
f*(elm % incoherent(i_grid+1) - elm % incoherent(i_grid)))
|
||||
|
||||
! Calculate microscopic photoelectric cross section
|
||||
micro_photon_xs(i_element) % photoelectric = exp(elm % photoelectric_total(&
|
||||
i_grid) + f*(elm % photoelectric_total(i_grid+1) - &
|
||||
elm % photoelectric_total(i_grid)))
|
||||
|
||||
! Calculate microscopic pair production cross section
|
||||
micro_photon_xs(i_element) % pair_production = exp(&
|
||||
elm % pair_production_total(i_grid) + f*(&
|
||||
elm % pair_production_total(i_grid+1) - &
|
||||
elm % pair_production_total(i_grid)))
|
||||
|
||||
! Calculate microscopic total cross section
|
||||
micro_photon_xs(i_element) % total = &
|
||||
micro_photon_xs(i_element) % coherent + &
|
||||
micro_photon_xs(i_element) % incoherent + &
|
||||
micro_photon_xs(i_element) % photoelectric + &
|
||||
micro_photon_xs(i_element) % pair_production
|
||||
|
||||
micro_photon_xs(i_element) % last_E = E
|
||||
end associate
|
||||
|
||||
end subroutine calculate_element_xs
|
||||
|
||||
!===============================================================================
|
||||
! FIND_ENERGY_INDEX determines the index on the union energy grid at a certain
|
||||
! energy
|
||||
|
|
|
|||
|
|
@ -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
|
||||
use photon_header, only: PhotonInteraction, ElementMicroXS
|
||||
use plot_header, only: ObjectPlot
|
||||
use sab_header, only: SAlphaBeta
|
||||
use set_header, only: SetInt
|
||||
|
|
@ -79,6 +79,7 @@ module global
|
|||
|
||||
! Cross section caches
|
||||
type(NuclideMicroXS), allocatable :: micro_xs(:) ! Cache for each nuclide
|
||||
type(ElementMicroXS), allocatable :: micro_photon_xs(:) ! Cache for each element
|
||||
type(MaterialMacroXS) :: material_xs ! Cache for current material
|
||||
|
||||
! Dictionaries to look up cross sections and listings
|
||||
|
|
@ -452,8 +453,8 @@ module global
|
|||
real(8) :: res_scat_energy_max = 1000.0_8
|
||||
character(10), allocatable :: res_scat_nuclides(:)
|
||||
|
||||
!$omp threadprivate(micro_xs, material_xs, fission_bank, n_bank, &
|
||||
!$omp& trace, thread_id, current_work, filter_matches)
|
||||
!$omp threadprivate(micro_xs, micro_photon_xs, material_xs, fission_bank, &
|
||||
!$omp& n_bank, trace, thread_id, current_work, filter_matches)
|
||||
|
||||
contains
|
||||
|
||||
|
|
|
|||
|
|
@ -163,13 +163,13 @@ contains
|
|||
integer, intent(in) :: intracomm ! MPI intracommunicator
|
||||
#endif
|
||||
|
||||
integer :: bank_blocks(5) ! Count for each datatype
|
||||
integer :: bank_blocks(6) ! Count for each datatype
|
||||
#ifdef MPIF08
|
||||
type(MPI_Datatype) :: bank_types(5)
|
||||
type(MPI_Datatype) :: bank_types(6)
|
||||
#else
|
||||
integer :: bank_types(5) ! Datatypes
|
||||
integer :: bank_types(6) ! Datatypes
|
||||
#endif
|
||||
integer(MPI_ADDRESS_KIND) :: bank_disp(5) ! Displacements
|
||||
integer(MPI_ADDRESS_KIND) :: bank_disp(6) ! Displacements
|
||||
logical :: init_called
|
||||
type(Bank) :: b
|
||||
|
||||
|
|
@ -201,14 +201,16 @@ contains
|
|||
call MPI_GET_ADDRESS(b % uvw, bank_disp(3), mpi_err)
|
||||
call MPI_GET_ADDRESS(b % E, bank_disp(4), mpi_err)
|
||||
call MPI_GET_ADDRESS(b % delayed_group, bank_disp(5), mpi_err)
|
||||
call MPI_GET_ADDRESS(b % particle, bank_disp(6), mpi_err)
|
||||
|
||||
! Adjust displacements
|
||||
bank_disp = bank_disp - bank_disp(1)
|
||||
|
||||
! Define MPI_BANK for fission sites
|
||||
bank_blocks = (/ 1, 3, 3, 1, 1 /)
|
||||
bank_types = (/ MPI_REAL8, MPI_REAL8, MPI_REAL8, MPI_REAL8, MPI_INTEGER /)
|
||||
call MPI_TYPE_CREATE_STRUCT(5, bank_blocks, bank_disp, &
|
||||
bank_blocks = (/ 1, 3, 3, 1, 1, 1 /)
|
||||
bank_types = (/ MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, &
|
||||
MPI_INT, MPI_INT /)
|
||||
call MPI_TYPE_CREATE_STRUCT(6, bank_blocks, bank_disp, &
|
||||
bank_types, MPI_BANK, mpi_err)
|
||||
call MPI_TYPE_COMMIT(MPI_BANK, mpi_err)
|
||||
|
||||
|
|
@ -245,6 +247,8 @@ contains
|
|||
c_loc(tmpb(1)%E)), H5T_NATIVE_DOUBLE, hdf5_err)
|
||||
call h5tinsert_f(hdf5_bank_t, "delayed_group", h5offsetof(c_loc(tmpb(1)), &
|
||||
c_loc(tmpb(1)%delayed_group)), H5T_NATIVE_INTEGER, hdf5_err)
|
||||
call h5tinsert_f(hdf5_bank_t, "particle", h5offsetof(c_loc(tmpb(1)), &
|
||||
c_loc(tmpb(1)%particle)), H5T_NATIVE_INTEGER, hdf5_err)
|
||||
|
||||
! Determine type for integer(8)
|
||||
hdf5_integer8_t = h5kind_to_type(8, H5_INTEGER_KIND)
|
||||
|
|
|
|||
|
|
@ -143,6 +143,12 @@ module nuclide_header
|
|||
real(8) :: absorption ! macroscopic absorption xs
|
||||
real(8) :: fission ! macroscopic fission xs
|
||||
real(8) :: nu_fission ! macroscopic production xs
|
||||
|
||||
! Photon cross sections
|
||||
real(8) :: coherent ! macroscopic coherent xs
|
||||
real(8) :: incoherent ! macroscopic incoherent xs
|
||||
real(8) :: photoelectric ! macroscopic photoelectric xs
|
||||
real(8) :: pair_production ! macroscopic pair production xs
|
||||
end type MaterialMacroXS
|
||||
|
||||
!===============================================================================
|
||||
|
|
@ -155,7 +161,7 @@ module nuclide_header
|
|||
character(MAX_FILE_LEN) :: path
|
||||
end type Library
|
||||
|
||||
contains
|
||||
contains
|
||||
|
||||
!===============================================================================
|
||||
! NUCLIDE_CLEAR resets and deallocates data in Nuclide
|
||||
|
|
|
|||
|
|
@ -199,6 +199,7 @@ contains
|
|||
call this % initialize()
|
||||
|
||||
! copy attributes from source bank site
|
||||
this % type = src % particle
|
||||
this % wgt = src % wgt
|
||||
this % last_wgt = src % wgt
|
||||
this % coord(1) % xyz = src % xyz
|
||||
|
|
@ -222,9 +223,10 @@ contains
|
|||
! the secondary bank and increments the number of sites in the secondary bank.
|
||||
!===============================================================================
|
||||
|
||||
subroutine create_secondary(this, uvw, type, run_CE)
|
||||
subroutine create_secondary(this, uvw, E, type, run_CE)
|
||||
class(Particle), intent(inout) :: this
|
||||
real(8), intent(in) :: uvw(3)
|
||||
real(8), intent(in) :: E
|
||||
integer, intent(in) :: type
|
||||
logical, intent(in) :: run_CE
|
||||
|
||||
|
|
@ -237,14 +239,15 @@ contains
|
|||
end if
|
||||
|
||||
n = this % n_secondary + 1
|
||||
this % secondary_bank(n) % wgt = this % wgt
|
||||
this % secondary_bank(n) % particle = type
|
||||
this % secondary_bank(n) % wgt = this % wgt
|
||||
this % secondary_bank(n) % xyz(:) = this % coord(1) % xyz
|
||||
this % secondary_bank(n) % uvw(:) = uvw
|
||||
this % n_secondary = n
|
||||
this % secondary_bank(this % n_secondary) % E = this % E
|
||||
this % secondary_bank(n) % E = E
|
||||
if (.not. run_CE) then
|
||||
this % secondary_bank(this % n_secondary) % E = real(this % g, 8)
|
||||
this % secondary_bank(n) % E = real(this % g, 8)
|
||||
end if
|
||||
this % n_secondary = n
|
||||
|
||||
end subroutine create_secondary
|
||||
|
||||
|
|
|
|||
|
|
@ -56,6 +56,22 @@ module photon_header
|
|||
procedure :: from_hdf5 => photon_from_hdf5
|
||||
end type PhotonInteraction
|
||||
|
||||
!===============================================================================
|
||||
! ELEMENTMICROXS contains cached microscopic photon cross sections for a
|
||||
! particular element at the current energy
|
||||
!===============================================================================
|
||||
|
||||
type ElementMicroXS
|
||||
integer :: index_grid ! index on element energy grid
|
||||
real(8) :: last_E = ZERO ! last evaluated energy
|
||||
real(8) :: interp_factor ! interpolation factor on energy grid
|
||||
real(8) :: total ! microscropic total photon xs
|
||||
real(8) :: coherent ! microscopic coherent xs
|
||||
real(8) :: incoherent ! microscopic incoherent xs
|
||||
real(8) :: photoelectric ! microscopic photoelectric xs
|
||||
real(8) :: pair_production ! microscopic pair production xs
|
||||
end type ElementMicroXS
|
||||
|
||||
contains
|
||||
|
||||
subroutine photon_from_hdf5(this, group_id)
|
||||
|
|
|
|||
365
src/photon_physics.F90
Normal file
365
src/photon_physics.F90
Normal file
|
|
@ -0,0 +1,365 @@
|
|||
module photon_physics
|
||||
|
||||
use algorithm, only: binary_search
|
||||
use constants
|
||||
use particle_header, only: Particle
|
||||
use photon_header, only: PhotonInteraction, compton_profile_pz
|
||||
use random_lcg, only: prn
|
||||
|
||||
contains
|
||||
|
||||
!===============================================================================
|
||||
! KLEIN_NISHINA
|
||||
!===============================================================================
|
||||
|
||||
subroutine klein_nishina(alpha, alpha_out, mu)
|
||||
real(8), intent(in) :: alpha
|
||||
real(8), intent(out) :: alpha_out
|
||||
real(8), intent(out) :: mu
|
||||
|
||||
real(8) :: beta ! 1 + 2a
|
||||
real(8) :: t ! (1 + 2a)/(9 + 2a)
|
||||
real(8) :: r, s, x
|
||||
real(8) :: gamma
|
||||
|
||||
beta = ONE + TWO*alpha
|
||||
if (alpha < THREE) then
|
||||
! Kahn's rejection method
|
||||
t = beta/(beta + 8.0_8)
|
||||
do
|
||||
if (prn() < t) then
|
||||
! Left branch of flow chart
|
||||
r = TWO*prn()
|
||||
x = ONE + alpha*r
|
||||
if (prn() < FOUR/x*(ONE - ONE/x)) then
|
||||
mu = 1 - r
|
||||
exit
|
||||
end if
|
||||
else
|
||||
! Right branch of flow chart
|
||||
x = beta/(ONE + TWO*alpha*prn())
|
||||
mu = ONE + (ONE - x)/alpha
|
||||
if (prn() < HALF*(mu**2 + ONE/x)) exit
|
||||
end if
|
||||
end do
|
||||
alpha_out = alpha/x
|
||||
|
||||
else
|
||||
! Koblinger's direct method
|
||||
gamma = ONE - beta**(-2)
|
||||
s = prn()*(FOUR/alpha + HALF*gamma + &
|
||||
(ONE - (ONE + beta)/alpha**2)*log(beta))
|
||||
if (s <= 2./alpha) then
|
||||
! For first term, x = 1 + 2ar
|
||||
! Therefore, a' = a/(1 + 2ar)
|
||||
alpha_out = alpha/(ONE + TWO*alpha*prn())
|
||||
elseif (s <= FOUR/alpha) then
|
||||
! For third term, x = beta/(1 + 2ar)
|
||||
! Therefore, a' = a(1 + 2ar)/beta
|
||||
alpha_out = alpha*(ONE + TWO*alpha*prn())/beta
|
||||
elseif (s <= FOUR/alpha + HALF*gamma) then
|
||||
! For fourth term, x = 1/sqrt(1 - gamma*r)
|
||||
! Therefore, a' = a*sqrt(1 - gamma*r)
|
||||
alpha_out = alpha*sqrt(ONE - gamma*prn())
|
||||
else
|
||||
! For third term, x = beta^r
|
||||
! Therefore, a' = a/beta^r
|
||||
alpha_out = alpha/beta**prn()
|
||||
end if
|
||||
|
||||
! Calculate cosine of scattering angle based on basic relation
|
||||
mu = ONE + ONE/alpha - ONE/alpha_out
|
||||
end if
|
||||
|
||||
end subroutine klein_nishina
|
||||
|
||||
!===============================================================================
|
||||
! COMPTON_SCATTER
|
||||
!===============================================================================
|
||||
|
||||
subroutine compton_scatter(el, alpha, alpha_out, mu, use_doppler)
|
||||
type(PhotonInteraction), intent(in) :: el
|
||||
real(8), intent(in) :: alpha
|
||||
real(8), intent(out) :: alpha_out
|
||||
real(8), intent(out) :: mu
|
||||
logical, intent(in), optional :: use_doppler
|
||||
|
||||
real(8) :: x
|
||||
real(8) :: form_factor_xmax
|
||||
real(8) :: form_factor_x
|
||||
real(8) :: e_out
|
||||
logical :: use_doppler_
|
||||
|
||||
if (present(use_doppler)) then
|
||||
use_doppler_ = use_doppler
|
||||
else
|
||||
use_doppler_ = .false.
|
||||
end if
|
||||
|
||||
form_factor_xmax = ZERO
|
||||
do
|
||||
! Sample Klein-Nishina distribution for trial energy and angle
|
||||
call klein_nishina(alpha, alpha_out, mu)
|
||||
|
||||
! Note that the parameter used here does not correspond exactly to the
|
||||
! momentum transfer q in ENDF-102 Eq. (27.2). Rather, this is the
|
||||
! parameter as defined by Hubbell, where the actual data comes from
|
||||
x = MASS_ELECTRON/PLANCK_C*alpha*sqrt(HALF*(ONE - mu))
|
||||
|
||||
! Calculate S(x, Z) and S(x_max, Z)
|
||||
form_factor_x = el % incoherent_form_factor % evaluate(x)
|
||||
if (form_factor_xmax == ZERO) then
|
||||
form_factor_xmax = el % incoherent_form_factor % evaluate(&
|
||||
MASS_ELECTRON/PLANCK_C*alpha)
|
||||
end if
|
||||
|
||||
! 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)
|
||||
alpha_out = e_out/MASS_ELECTRON
|
||||
end if
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
|
||||
end subroutine compton_scatter
|
||||
|
||||
!===============================================================================
|
||||
! COMPTON_DOPPLER
|
||||
!===============================================================================
|
||||
|
||||
subroutine compton_doppler(el, alpha, mu, e_out)
|
||||
type(PhotonInteraction), intent(in) :: el
|
||||
real(8), intent(in) :: alpha
|
||||
real(8), intent(in) :: mu
|
||||
real(8), intent(out) :: e_out
|
||||
|
||||
integer :: i, i_shell
|
||||
integer :: n
|
||||
real(8) :: rn, m
|
||||
real(8) :: c, c_l, c_max
|
||||
real(8) :: pz_l, pz_r, pz, pz_max
|
||||
real(8) :: p_l, p_r
|
||||
real(8) :: e, e_b
|
||||
real(8) :: e_out1, e_out2
|
||||
real(8) :: a, b, quad
|
||||
real(8) :: f
|
||||
real(8) :: momentum_sq
|
||||
|
||||
n = size(compton_profile_pz)
|
||||
|
||||
do
|
||||
! Sample electron shell
|
||||
rn = prn()
|
||||
c = ZERO
|
||||
do i_shell = 1, size(el % electron_pdf) - 1
|
||||
c = c + el % electron_pdf(i_shell + 1)
|
||||
if (rn < c) exit
|
||||
end do
|
||||
|
||||
! Determine binding energy of shell
|
||||
e_b = el % binding_energy(i_shell)
|
||||
|
||||
! Determine p_z,max
|
||||
e = alpha*MASS_ELECTRON
|
||||
if (e < e_b) then
|
||||
e_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON
|
||||
exit
|
||||
end if
|
||||
|
||||
pz_max = -FINE_STRUCTURE*(e_b - (e - e_b)*alpha*(ONE - mu)) / &
|
||||
sqrt(TWO*e*(e - e_b)*(ONE - mu) + e_b**2)
|
||||
if (pz_max < ZERO) then
|
||||
e_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON
|
||||
exit
|
||||
end if
|
||||
|
||||
! Determine profile cdf value corresponding to p_z,max
|
||||
if (pz_max > compton_profile_pz(n)) then
|
||||
c_max = el % profile_cdf(n, i_shell)
|
||||
else
|
||||
i = binary_search(compton_profile_pz, n, pz_max)
|
||||
pz_l = compton_profile_pz(i)
|
||||
pz_r = compton_profile_pz(i + 1)
|
||||
p_l = el % profile_pdf(i, i_shell)
|
||||
p_r = el % profile_pdf(i + 1, i_shell)
|
||||
c_l = el % profile_cdf(i, i_shell)
|
||||
if (pz_l == pz_r) then
|
||||
c_max = c_l
|
||||
elseif (p_l == p_r) then
|
||||
c_max = c_l + (pz_max - pz_l)*p_l
|
||||
else
|
||||
m = (p_l - p_r)/(pz_l - pz_r)
|
||||
c_max = c_l + ((m*(pz_max - pz_l) + p_l)**2 - p_l**2)/(TWO*m)
|
||||
end if
|
||||
end if
|
||||
|
||||
! Sample value on bounded cdf
|
||||
c = prn()*c_max
|
||||
|
||||
! Determine pz corresponding to sampled cdf value
|
||||
i = binary_search(el % profile_cdf(:, i_shell), n, c)
|
||||
pz_l = compton_profile_pz(i)
|
||||
pz_r = compton_profile_pz(i + 1)
|
||||
p_l = el % profile_pdf(i, i_shell)
|
||||
p_r = el % profile_pdf(i + 1, i_shell)
|
||||
c_l = el % profile_cdf(i, i_shell)
|
||||
if (pz_l == pz_r) then
|
||||
pz = pz_l
|
||||
elseif (p_l == p_r) then
|
||||
pz = pz_l + (c - c_l)/p_l
|
||||
else
|
||||
m = (p_l - p_r)/(pz_l - pz_r)
|
||||
pz = pz_l + (sqrt(p_l**2 + TWO*m*(c - c_l)) - p_l)/m
|
||||
end if
|
||||
|
||||
! Determine outgoing photon energy corresponding to electron momentum
|
||||
momentum_sq = (pz/FINE_STRUCTURE)**2
|
||||
f = ONE + alpha*(ONE - mu)
|
||||
a = momentum_sq - f*f
|
||||
b = TWO*e*(f - momentum_sq*mu)
|
||||
c = e**2*(momentum_sq - ONE)
|
||||
|
||||
quad = b**2 - FOUR*a*c
|
||||
if (quad < 0) then
|
||||
e_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON
|
||||
exit
|
||||
end if
|
||||
quad = sqrt(quad)
|
||||
e_out1 = -(b + quad)/(TWO*a)
|
||||
e_out2 = -(b - quad)/(TWO*a)
|
||||
|
||||
if (e_out1 > ZERO) then
|
||||
if (e_out2 > ZERO) then
|
||||
if (prn() < HALF) then
|
||||
e_out = e_out1
|
||||
else
|
||||
e_out = e_out2
|
||||
end if
|
||||
else
|
||||
e_out = e_out1
|
||||
end if
|
||||
else
|
||||
if (e_out2 > ZERO) e_out = e_out2
|
||||
end if
|
||||
if (e_out < e - e_b) exit
|
||||
end do
|
||||
|
||||
end subroutine compton_doppler
|
||||
|
||||
!===============================================================================
|
||||
! RAYLEIGH_SCATTER
|
||||
!===============================================================================
|
||||
|
||||
subroutine rayleigh_scatter(el, alpha, mu)
|
||||
type(PhotonInteraction), intent(in) :: el
|
||||
real(8), intent(in) :: alpha
|
||||
real(8), intent(out) :: mu
|
||||
|
||||
integer :: i
|
||||
real(8) :: F
|
||||
real(8) :: F_max
|
||||
real(8) :: x2
|
||||
real(8) :: x2_max
|
||||
real(8) :: r
|
||||
|
||||
do
|
||||
! Determine maximum value of x^2
|
||||
x2_max = (MASS_ELECTRON/PLANCK_C*alpha)**2
|
||||
|
||||
! Determine F(x^2_max, Z)
|
||||
F_max = el % coherent_int_form_factor % evaluate(x2_max)
|
||||
|
||||
! Sample cumulative distribution
|
||||
F = prn()*F_max
|
||||
|
||||
! Determine x^2 corresponding to F
|
||||
i = binary_search(el%coherent_int_form_factor%y, &
|
||||
size(el%coherent_int_form_factor%y), F)
|
||||
r = (F - el%coherent_int_form_factor%y(i)) / &
|
||||
(el%coherent_int_form_factor%y(i+1) - el%coherent_int_form_factor%y(i))
|
||||
x2 = el%coherent_int_form_factor%x(i) + r*(el%coherent_int_form_factor%x(i+1) - &
|
||||
el%coherent_int_form_factor%x(i))
|
||||
|
||||
! Calculate mu
|
||||
mu = ONE - TWO*x2/x2_max
|
||||
|
||||
if (prn() < HALF*(ONE + mu**2)) exit
|
||||
end do
|
||||
|
||||
end subroutine rayleigh_scatter
|
||||
|
||||
!===============================================================================
|
||||
! ATOMIC_RELAXATION
|
||||
!===============================================================================
|
||||
|
||||
recursive subroutine atomic_relaxation(p, elm, i_shell)
|
||||
type(Particle), intent(inout) :: p
|
||||
type(PhotonInteraction), intent(in) :: elm
|
||||
integer, intent(in) :: i_shell
|
||||
|
||||
integer :: i_hole
|
||||
integer :: i_transition
|
||||
integer :: primary
|
||||
integer :: secondary
|
||||
real(8) :: c
|
||||
real(8) :: rn
|
||||
real(8) :: E
|
||||
real(8) :: mu
|
||||
real(8) :: phi
|
||||
real(8) :: uvw(3)
|
||||
|
||||
! Check for no transitions
|
||||
if (elm % shells(i_shell) % n_transitions == 0) return
|
||||
|
||||
! Sample transition
|
||||
rn = prn()
|
||||
c = ZERO
|
||||
do i_transition = 1, elm % shells(i_shell) % n_transitions - 1
|
||||
c = c + elm % shells(i_shell) % &
|
||||
transition_probability(i_transition + 1)
|
||||
if (rn < c) exit
|
||||
end do
|
||||
|
||||
! Get primary and secondary subshell designators
|
||||
primary = elm % shells(i_shell) % transition_subshells(1, i_transition)
|
||||
secondary = elm % shells(i_shell) % transition_subshells(2, i_transition)
|
||||
|
||||
if (secondary == 0) then
|
||||
! Non-radiative trnasition -- Auger/Coster-Kronig effect
|
||||
|
||||
! TODO: Create electron
|
||||
! E_electron = transition_energy(i_transition)
|
||||
|
||||
! Fill secondary (higher) hole first
|
||||
if (elm % shell_dict % has_key(secondary)) then
|
||||
i_hole = elm % shell_dict % get_key(secondary)
|
||||
call atomic_relaxation(p, elm, i_hole)
|
||||
end if
|
||||
else
|
||||
! Radiative transition -- get X-ray energy
|
||||
E = elm % shells(i_shell) % transition_energy(i_transition)
|
||||
|
||||
if (E > ZERO) then
|
||||
! Sample angle isotropically for X-ray
|
||||
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)
|
||||
|
||||
! Create X-ray
|
||||
call p % create_secondary(uvw, E, PHOTON, run_ce=.true.)
|
||||
end if
|
||||
end if
|
||||
|
||||
! Fill primary hole
|
||||
if (elm % shell_dict % has_key(primary)) then
|
||||
i_hole = elm % shell_dict % get_key(primary)
|
||||
call atomic_relaxation(p, elm, i_hole)
|
||||
end if
|
||||
|
||||
end subroutine atomic_relaxation
|
||||
|
||||
end module photon_physics
|
||||
163
src/physics.F90
163
src/physics.F90
|
|
@ -14,6 +14,8 @@ module physics
|
|||
use output, only: write_message
|
||||
use particle_header, only: Particle
|
||||
use particle_restart_write, only: write_particle_restart
|
||||
use photon_physics, only: rayleigh_scatter, compton_scatter, &
|
||||
atomic_relaxation
|
||||
use physics_common
|
||||
use random_lcg, only: prn, advance_prn_seed, prn_set_stream
|
||||
use reaction_header, only: Reaction
|
||||
|
|
@ -41,8 +43,12 @@ contains
|
|||
! Add to collision counter for particle
|
||||
p % n_collision = p % n_collision + 1
|
||||
|
||||
! Sample nuclide/reaction for the material the particle is in
|
||||
call sample_reaction(p)
|
||||
! Sample reaction for the material the particle is in
|
||||
if (p % type == NEUTRON) then
|
||||
call sample_reaction(p)
|
||||
elseif (p % type == PHOTON) then
|
||||
call sample_photon_reaction(p)
|
||||
end if
|
||||
|
||||
! Display information about collision
|
||||
if (verbosity >= 10 .or. trace) then
|
||||
|
|
@ -137,6 +143,110 @@ contains
|
|||
|
||||
end subroutine sample_reaction
|
||||
|
||||
!===============================================================================
|
||||
! SAMPLE_PHOTON_REACTION samples an element based on the macroscopic cross
|
||||
! sections for each nuclide within a material and then samples a reaction for
|
||||
! that element and calls the appropriate routine to process the physics.
|
||||
!===============================================================================
|
||||
|
||||
subroutine sample_photon_reaction(p)
|
||||
type(Particle), intent(inout) :: p
|
||||
|
||||
integer :: i_shell ! index in subshells
|
||||
integer :: i_grid ! index on energy grid
|
||||
integer :: i_element ! index in nuclides array
|
||||
integer :: i_start ! threshold index
|
||||
real(8) :: prob ! cumulative probability
|
||||
real(8) :: cutoff ! sampled total cross section
|
||||
real(8) :: f ! interpolation factor
|
||||
real(8) :: xs ! photoionization cross section
|
||||
real(8) :: prob_after
|
||||
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) :: phi ! azimuthal angle
|
||||
|
||||
! Sample element within material
|
||||
i_element = sample_element(p)
|
||||
p % event_nuclide = i_element
|
||||
|
||||
! Calculate photon energy over electron rest mass equivalent
|
||||
alpha = p % E/MASS_ELECTRON
|
||||
|
||||
! For tallying purposes, this routine might be called directly. In that
|
||||
! case, we need to sample a reaction via the cutoff variable
|
||||
prob = ZERO
|
||||
cutoff = prn() * micro_photon_xs(i_element) % total
|
||||
|
||||
associate (elm => elements(i_element))
|
||||
! Coherent (Rayleigh) scattering
|
||||
prob = prob + micro_photon_xs(i_element) % coherent
|
||||
if (prob > cutoff) then
|
||||
call rayleigh_scatter(elm, alpha, mu)
|
||||
p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu)
|
||||
return
|
||||
end if
|
||||
|
||||
! Incoherent (Compton) scattering
|
||||
prob = prob + micro_photon_xs(i_element) % incoherent
|
||||
if (prob > cutoff) then
|
||||
call compton_scatter(elm, alpha, alpha_out, mu, .true.)
|
||||
p % E = alpha_out*MASS_ELECTRON
|
||||
p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu)
|
||||
return
|
||||
end if
|
||||
|
||||
! Photoelectric effect
|
||||
prob_after = prob + micro_photon_xs(i_element) % photoelectric
|
||||
if (prob_after > cutoff) then
|
||||
do i_shell = 1, size(elm % shells)
|
||||
! Get grid index and interpolation factor
|
||||
i_grid = micro_photon_xs(i_element) % index_grid
|
||||
f = micro_photon_xs(i_element) % interp_factor
|
||||
|
||||
! Check threshold of reaction
|
||||
i_start = elm % shells(i_shell) % threshold
|
||||
if (i_grid <= i_start) cycle
|
||||
|
||||
! Evaluation subshell photoionization cross section
|
||||
xs = exp(elm % shells(i_shell) % cross_section(i_grid - i_start) + &
|
||||
f*(elm % shells(i_shell) % cross_section(i_grid + 1 - i_start) - &
|
||||
elm % shells(i_shell) % cross_section(i_grid - i_start)))
|
||||
|
||||
prob = prob + xs
|
||||
if (prob > cutoff) then
|
||||
! TODO: Create electron
|
||||
! E_electron = p % E - elm % shells(i_shell) % binding_energy
|
||||
|
||||
call atomic_relaxation(p, elm, i_shell)
|
||||
p % alive = .false.
|
||||
return
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
prob = prob_after
|
||||
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()
|
||||
p % coord(1) % uvw(1) = mu
|
||||
p % coord(1) % uvw(2) = sqrt(ONE - mu*mu)*cos(phi)
|
||||
p % coord(1) % uvw(3) = sqrt(ONE - mu*mu)*sin(phi)
|
||||
|
||||
! Set energy
|
||||
p % E = MASS_ELECTRON
|
||||
|
||||
! Create photon in opposite direction
|
||||
call p % create_secondary(-p % coord(1) % uvw, MASS_ELECTRON, &
|
||||
PHOTON, .true.)
|
||||
end if
|
||||
|
||||
end subroutine sample_photon_reaction
|
||||
|
||||
!===============================================================================
|
||||
! SAMPLE_NUCLIDE
|
||||
!===============================================================================
|
||||
|
|
@ -199,6 +309,49 @@ contains
|
|||
|
||||
end subroutine sample_nuclide
|
||||
|
||||
!===============================================================================
|
||||
! SAMPLE_ELEMENT
|
||||
!===============================================================================
|
||||
|
||||
function sample_element(p) result(i_element)
|
||||
type(Particle), intent(in) :: p
|
||||
integer :: i_element
|
||||
|
||||
integer :: i
|
||||
real(8) :: prob
|
||||
real(8) :: cutoff
|
||||
real(8) :: atom_density ! atom density of nuclide in atom/b-cm
|
||||
real(8) :: sigma ! microscopic total xs for nuclide
|
||||
|
||||
associate (mat => materials(p % material))
|
||||
! Sample cumulative distribution function
|
||||
cutoff = prn() * material_xs % total
|
||||
|
||||
i = 0
|
||||
prob = ZERO
|
||||
do while (prob < cutoff)
|
||||
i = i + 1
|
||||
|
||||
! Check to make sure that a nuclide was sampled
|
||||
if (i > mat % n_nuclides) then
|
||||
call write_particle_restart(p)
|
||||
call fatal_error("Did not sample any element during collision.")
|
||||
end if
|
||||
|
||||
! Find atom density
|
||||
i_element = mat % element(i)
|
||||
atom_density = mat % atom_density(i)
|
||||
|
||||
! Determine microscopic cross section
|
||||
sigma = atom_density * micro_photon_xs(i_element) % total
|
||||
|
||||
! Increment probability to compare to cutoff
|
||||
prob = prob + sigma
|
||||
end do
|
||||
end associate
|
||||
|
||||
end function sample_element
|
||||
|
||||
!===============================================================================
|
||||
! SAMPLE_FISSION
|
||||
!===============================================================================
|
||||
|
|
@ -1140,6 +1293,9 @@ contains
|
|||
! Bank source neutrons by copying particle data
|
||||
bank_array(i) % xyz = p % coord(1) % xyz
|
||||
|
||||
! Set particle as neutron
|
||||
bank_array(i) % particle = NEUTRON
|
||||
|
||||
! Set weight of fission bank site
|
||||
bank_array(i) % wgt = ONE/weight
|
||||
|
||||
|
|
@ -1329,7 +1485,8 @@ contains
|
|||
if (mod(yield, ONE) == ZERO) then
|
||||
! If yield is integral, create exactly that many secondary particles
|
||||
do i = 1, nint(yield) - 1
|
||||
call p % create_secondary(p % coord(1) % uvw, NEUTRON, run_CE=.true.)
|
||||
call p % create_secondary(p % coord(1) % uvw, p % E, &
|
||||
NEUTRON, run_CE=.true.)
|
||||
end do
|
||||
else
|
||||
! Otherwise, change weight of particle based on yield
|
||||
|
|
|
|||
|
|
@ -378,6 +378,7 @@ contains
|
|||
|
||||
!$omp parallel
|
||||
allocate(micro_xs(n_nuclides_total))
|
||||
allocate(micro_photon_xs(n_elements))
|
||||
!$omp end parallel
|
||||
|
||||
if (.not. restart_run) call initialize_source()
|
||||
|
|
@ -403,6 +404,7 @@ contains
|
|||
|
||||
!$omp parallel
|
||||
deallocate(micro_xs)
|
||||
deallocate(micro_photon_xs)
|
||||
!$omp end parallel
|
||||
|
||||
! Increment total number of generations
|
||||
|
|
|
|||
|
|
@ -134,6 +134,9 @@ contains
|
|||
! Set particle defaults
|
||||
call p % initialize()
|
||||
|
||||
! Set particle type
|
||||
site % particle = external_source(i) % particle
|
||||
|
||||
! Sample spatial distribution
|
||||
site % xyz(:) = external_source(i) % space % sample()
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue