mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-29 06:35:48 -04:00
implemented sampling of secondary photons from neutron collisons
This commit is contained in:
samuelshaner
parent
686362248b
commit
d1b8efc5aa
4 changed files with 225 additions and 27 deletions
|
|
@ -53,11 +53,12 @@ contains
|
|||
logical :: check_sab ! should we check for S(a,b) table?
|
||||
|
||||
! Set all material macroscopic cross sections to zero
|
||||
material_xs % total = ZERO
|
||||
material_xs % elastic = ZERO
|
||||
material_xs % absorption = ZERO
|
||||
material_xs % fission = ZERO
|
||||
material_xs % nu_fission = ZERO
|
||||
material_xs % total = ZERO
|
||||
material_xs % elastic = ZERO
|
||||
material_xs % absorption = ZERO
|
||||
material_xs % fission = ZERO
|
||||
material_xs % nu_fission = ZERO
|
||||
material_xs % nu_photon_total = ZERO
|
||||
|
||||
! Exit subroutine if material is void
|
||||
if (p % material == MATERIAL_VOID) return
|
||||
|
|
@ -137,6 +138,10 @@ contains
|
|||
! Add contributions to material macroscopic nu-fission cross section
|
||||
material_xs % nu_fission = material_xs % nu_fission + &
|
||||
atom_density * micro_xs(i_nuclide) % nu_fission
|
||||
|
||||
! Add contributions to material macroscopic nu-fission cross section
|
||||
material_xs % nu_photon_total = material_xs % nu_photon_total + &
|
||||
atom_density * micro_xs(i_nuclide) % nu_photon_total
|
||||
end do
|
||||
end associate
|
||||
|
||||
|
|
@ -257,8 +262,9 @@ contains
|
|||
micro_xs(i_nuclide) % interp_factor = f
|
||||
|
||||
! Initialize nuclide cross-sections to zero
|
||||
micro_xs(i_nuclide) % fission = ZERO
|
||||
micro_xs(i_nuclide) % nu_fission = ZERO
|
||||
micro_xs(i_nuclide) % fission = ZERO
|
||||
micro_xs(i_nuclide) % nu_fission = ZERO
|
||||
micro_xs(i_nuclide) % nu_photon_total = ZERO
|
||||
|
||||
! Calculate microscopic nuclide total cross section
|
||||
micro_xs(i_nuclide) % total = (ONE - f) * xs % total(i_grid) &
|
||||
|
|
@ -272,6 +278,10 @@ contains
|
|||
micro_xs(i_nuclide) % absorption = (ONE - f) * xs % absorption( &
|
||||
i_grid) + f * xs % absorption(i_grid + 1)
|
||||
|
||||
! Calculate microscopic nuclide nu-photon total cross section
|
||||
micro_xs(i_nuclide) % nu_photon_total = (ONE - f) * xs % &
|
||||
nu_photon_total(i_grid) + f * xs % nu_photon_total(i_grid + 1)
|
||||
|
||||
if (nuc % fissionable) then
|
||||
! Calculate microscopic nuclide total cross section
|
||||
micro_xs(i_nuclide) % fission = (ONE - f) * xs % fission(i_grid) &
|
||||
|
|
@ -442,6 +452,7 @@ contains
|
|||
integer :: i_energy ! index for energy
|
||||
integer :: i_low ! band index at lower bounding energy
|
||||
integer :: i_up ! band index at upper bounding energy
|
||||
integer :: i_grid ! index on nuclide energy grid
|
||||
real(8) :: f ! interpolation factor
|
||||
real(8) :: r ! pseudo-random number
|
||||
real(8) :: elastic ! elastic cross section
|
||||
|
|
@ -563,6 +574,11 @@ contains
|
|||
micro_xs(i_nuclide) % fission = fission
|
||||
micro_xs(i_nuclide) % total = elastic + inelastic + capture + fission
|
||||
|
||||
! Set the nu-photon production cross section
|
||||
i_grid = int(log(E/energy_min_neutron)/log_spacing)
|
||||
micro_xs(i_nuclide) % nu_photon_total = &
|
||||
nuc % compute_nu_photon_total(E, i_temp, i_grid)
|
||||
|
||||
! Determine nu-fission cross section
|
||||
if (nuc % fissionable) then
|
||||
micro_xs(i_nuclide) % nu_fission = nuc % nu(E, EMISSION_TOTAL) * &
|
||||
|
|
|
|||
|
|
@ -333,8 +333,8 @@ contains
|
|||
c_k = c_k1
|
||||
end do
|
||||
|
||||
! Check to make sure k is <= NP - 1
|
||||
k = min(k, n_energy_out - 1)
|
||||
! Check to make sure 1 <= k <= NP - 1
|
||||
k = max(1, min(k, n_energy_out - 1))
|
||||
|
||||
E_l_k = this%distribution(l)%e_out(k)
|
||||
p_l_k = this%distribution(l)%p(k)
|
||||
|
|
@ -361,7 +361,7 @@ contains
|
|||
end if
|
||||
|
||||
! Now interpolate between incident energy bins i and i + 1
|
||||
if (.not. histogram_interp) then
|
||||
if (.not. histogram_interp .and. n_energy_out > 1) then
|
||||
if (l == i) then
|
||||
E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1)
|
||||
else
|
||||
|
|
|
|||
|
|
@ -5,7 +5,7 @@ module nuclide_header
|
|||
|
||||
use hdf5, only: HID_T, HSIZE_T, SIZE_T
|
||||
|
||||
use algorithm, only: sort, find
|
||||
use algorithm, only: sort, find, binary_search
|
||||
use constants
|
||||
use dict_header, only: DictIntInt
|
||||
use endf, only: reaction_name, is_fission, is_disappearance
|
||||
|
|
@ -37,12 +37,13 @@ module nuclide_header
|
|||
end type EnergyGrid
|
||||
|
||||
type SumXS
|
||||
real(8), allocatable :: total(:) ! total cross section
|
||||
real(8), allocatable :: elastic(:) ! elastic scattering
|
||||
real(8), allocatable :: fission(:) ! fission
|
||||
real(8), allocatable :: nu_fission(:) ! neutron production
|
||||
real(8), allocatable :: absorption(:) ! absorption (MT > 100)
|
||||
real(8), allocatable :: heating(:) ! heating
|
||||
real(8), allocatable :: total(:) ! total cross section
|
||||
real(8), allocatable :: elastic(:) ! elastic scattering
|
||||
real(8), allocatable :: fission(:) ! fission
|
||||
real(8), allocatable :: nu_fission(:) ! neutron production
|
||||
real(8), allocatable :: absorption(:) ! absorption (MT > 100)
|
||||
real(8), allocatable :: heating(:) ! heating
|
||||
real(8), allocatable :: nu_photon_total(:) ! photon production
|
||||
end type SumXS
|
||||
|
||||
type :: Nuclide
|
||||
|
|
@ -99,6 +100,7 @@ module nuclide_header
|
|||
procedure :: from_hdf5 => nuclide_from_hdf5
|
||||
procedure :: init_grid => nuclide_init_grid
|
||||
procedure :: nu => nuclide_nu
|
||||
procedure :: compute_nu_photon_total => compute_nuclide_nu_photon_total
|
||||
procedure, private :: create_derived => nuclide_create_derived
|
||||
end type Nuclide
|
||||
|
||||
|
|
@ -117,6 +119,7 @@ module nuclide_header
|
|||
real(8) :: absorption ! microscopic absorption xs
|
||||
real(8) :: fission ! microscopic fission xs
|
||||
real(8) :: nu_fission ! microscopic production xs
|
||||
real(8) :: nu_photon_total ! microscopic photon production xs
|
||||
|
||||
! Information for S(a,b) use
|
||||
integer :: index_sab ! index in sab_tables (zero means no table)
|
||||
|
|
@ -138,11 +141,12 @@ module nuclide_header
|
|||
!===============================================================================
|
||||
|
||||
type MaterialMacroXS
|
||||
real(8) :: total ! macroscopic total xs
|
||||
real(8) :: elastic ! macroscopic elastic scattering xs
|
||||
real(8) :: absorption ! macroscopic absorption xs
|
||||
real(8) :: fission ! macroscopic fission xs
|
||||
real(8) :: nu_fission ! macroscopic production xs
|
||||
real(8) :: total ! macroscopic total xs
|
||||
real(8) :: elastic ! macroscopic elastic scattering xs
|
||||
real(8) :: absorption ! macroscopic absorption xs
|
||||
real(8) :: fission ! macroscopic fission xs
|
||||
real(8) :: nu_fission ! macroscopic production xs
|
||||
real(8) :: nu_photon_total ! macroscopic photon production xs
|
||||
|
||||
! Photon cross sections
|
||||
real(8) :: coherent ! macroscopic coherent xs
|
||||
|
|
@ -481,7 +485,7 @@ contains
|
|||
subroutine nuclide_create_derived(this)
|
||||
class(Nuclide), intent(inout) :: this
|
||||
|
||||
integer :: i, j, k
|
||||
integer :: i, j, k, l
|
||||
integer :: t
|
||||
integer :: m
|
||||
integer :: n
|
||||
|
|
@ -501,11 +505,13 @@ contains
|
|||
allocate(this % sum_xs(i) % fission(n_grid))
|
||||
allocate(this % sum_xs(i) % nu_fission(n_grid))
|
||||
allocate(this % sum_xs(i) % absorption(n_grid))
|
||||
allocate(this % sum_xs(i) % nu_photon_total(n_grid))
|
||||
this % sum_xs(i) % total(:) = ZERO
|
||||
this % sum_xs(i) % elastic(:) = ZERO
|
||||
this % sum_xs(i) % fission(:) = ZERO
|
||||
this % sum_xs(i) % nu_fission(:) = ZERO
|
||||
this % sum_xs(i) % absorption(:) = ZERO
|
||||
this % sum_xs(i) % nu_photon_total(:) = ZERO
|
||||
end do
|
||||
|
||||
i_fission = 0
|
||||
|
|
@ -539,6 +545,18 @@ contains
|
|||
this % sum_xs(t) % total(j:j+n-1) = this % sum_xs(t) % total(j:j+n-1) + &
|
||||
rx % xs(t) % value
|
||||
|
||||
! Calculate nu-photon total cross section
|
||||
do k = 1, size(rx % products)
|
||||
if (rx % products(k) % particle == PHOTON) then
|
||||
do l = 1, n
|
||||
this % sum_xs(t) % nu_photon_total(l+j-1) = &
|
||||
this % sum_xs(t) % nu_photon_total(l+j-1) + &
|
||||
rx % xs(t) % value(l) * rx % products(k) % &
|
||||
yield % evaluate(this % grid(t) % energy(l+j-1))
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
|
||||
! Add contribution to absorption cross section
|
||||
if (is_disappearance(rx % MT)) then
|
||||
this % sum_xs(t) % absorption(j:j+n-1) = this % sum_xs(t) % &
|
||||
|
|
@ -565,10 +583,6 @@ contains
|
|||
this % sum_xs(t) % absorption(j:j+n-1) = this % sum_xs(t) % &
|
||||
absorption(j:j+n-1) + rx % xs(t) % value
|
||||
|
||||
! If total fission reaction is present, there's no need to store the
|
||||
! reaction cross-section since it was copied to this % fission
|
||||
if (rx % MT == N_FISSION) deallocate(rx % xs(t) % value)
|
||||
|
||||
! Keep track of this reaction for easy searching later
|
||||
if (t == 1) then
|
||||
i_fission = i_fission + 1
|
||||
|
|
@ -690,6 +704,72 @@ contains
|
|||
|
||||
end function nuclide_nu
|
||||
|
||||
!===============================================================================
|
||||
! COMPUTE_NUCLIDE_NU_PHOTON_TOTAL is an interface to compute the number of
|
||||
! photons produced
|
||||
!===============================================================================
|
||||
|
||||
pure function compute_nuclide_nu_photon_total(this, E, t, i_log_union) result(nu_photon_total)
|
||||
class(Nuclide), intent(in) :: this
|
||||
real(8), intent(in) :: E
|
||||
integer, intent(in) :: t
|
||||
integer, intent(in) :: i_log_union
|
||||
real(8) :: rx_xs
|
||||
real(8) :: nu_photon_total
|
||||
real(8) :: f
|
||||
integer :: m, j, k
|
||||
integer :: i_grid, i_low, i_high
|
||||
|
||||
associate (grid => this % grid(t), xs => this % sum_xs(t))
|
||||
! Determine the energy grid index using a logarithmic mapping to
|
||||
! reduce the energy range over which a binary search needs to be
|
||||
! performed
|
||||
|
||||
if (E < grid % energy(1)) then
|
||||
i_grid = 1
|
||||
elseif (E > grid % energy(size(grid % energy))) then
|
||||
i_grid = size(grid % energy) - 1
|
||||
else
|
||||
! Determine bounding indices based on which equal log-spaced
|
||||
! interval the energy is in
|
||||
i_low = grid % grid_index(i_log_union)
|
||||
i_high = grid % grid_index(i_log_union + 1) + 1
|
||||
|
||||
! Perform binary search over reduced range
|
||||
i_grid = binary_search(grid % energy(i_low:i_high), &
|
||||
i_high - i_low + 1, E) + i_low - 1
|
||||
end if
|
||||
|
||||
! check for rare case where two energy points are the same
|
||||
if (grid % energy(i_grid) == grid % energy(i_grid + 1)) &
|
||||
i_grid = i_grid + 1
|
||||
|
||||
! calculate interpolation factor
|
||||
f = (E - grid % energy(i_grid)) / &
|
||||
(grid % energy(i_grid + 1) - grid % energy(i_grid))
|
||||
end associate
|
||||
|
||||
nu_photon_total = ZERO
|
||||
|
||||
! Calculate nu-photon total cross section
|
||||
do m = 1, size(this % reactions)
|
||||
associate (rx => this % reactions(m))
|
||||
j = rx % xs(t) % threshold
|
||||
do k = 1, size(rx % products)
|
||||
if (rx % products(k) % particle == PHOTON) then
|
||||
if (i_grid >= j) then
|
||||
rx_xs = (ONE - f) * rx % xs(t) % value(i_grid - j + 1) &
|
||||
+ f * rx % xs(t) % value(i_grid - j + 2)
|
||||
nu_photon_total = nu_photon_total + rx_xs * &
|
||||
rx % products(k) % yield % evaluate(E)
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
end associate
|
||||
end do
|
||||
|
||||
end function compute_nuclide_nu_photon_total
|
||||
|
||||
subroutine nuclide_init_grid(this, E_min, E_max, M)
|
||||
class(Nuclide), intent(inout) :: this
|
||||
real(8), intent(in) :: E_min ! Minimum energy in MeV
|
||||
|
|
|
|||
102
src/physics.F90
102
src/physics.F90
|
|
@ -113,6 +113,9 @@ contains
|
|||
end if
|
||||
end if
|
||||
|
||||
! Create secondary photons
|
||||
call sample_secondary_photons(p, i_nuclide)
|
||||
|
||||
! If survival biasing is being used, the following subroutine adjusts the
|
||||
! weight of the particle. Otherwise, it checks to see if absorption occurs
|
||||
|
||||
|
|
@ -421,6 +424,59 @@ contains
|
|||
|
||||
end subroutine sample_fission
|
||||
|
||||
!===============================================================================
|
||||
! SAMPLE_PHOTON_PRODUCT
|
||||
!===============================================================================
|
||||
|
||||
subroutine sample_photon_product(i_nuclide, E, i_reaction, i_product)
|
||||
integer, intent(in) :: i_nuclide ! index in nuclides array
|
||||
real(8), intent(in) :: E ! energy of neutron
|
||||
integer, intent(out) :: i_reaction ! index in nuc % reactions array
|
||||
integer, intent(out) :: i_product ! index in nuc % reactions array
|
||||
|
||||
integer :: i_grid
|
||||
integer :: i_temp
|
||||
integer :: threshold
|
||||
real(8) :: f
|
||||
real(8) :: prob
|
||||
real(8) :: cutoff
|
||||
real(8) :: yield
|
||||
type(Nuclide), pointer :: nuc
|
||||
|
||||
! Get pointer to nuclide
|
||||
nuc => nuclides(i_nuclide)
|
||||
|
||||
! Get grid index and interpolation factor and sample proton production cdf
|
||||
i_temp = micro_xs(i_nuclide) % index_temp
|
||||
i_grid = micro_xs(i_nuclide) % index_grid
|
||||
f = micro_xs(i_nuclide) % interp_factor
|
||||
cutoff = prn() * micro_xs(i_nuclide) % nu_photon_total
|
||||
prob = ZERO
|
||||
|
||||
! Loop through each reaction type
|
||||
REACTION_LOOP: do i_reaction = 1, size(nuc % reactions)
|
||||
associate (rx => nuc % reactions(i_reaction))
|
||||
do i_product = 1, size(rx % products)
|
||||
if (rx % products(i_product) % particle == PHOTON) then
|
||||
|
||||
threshold = rx % xs(i_temp) % threshold
|
||||
|
||||
! if energy is below threshold for this reaction, skip it
|
||||
if (i_grid < threshold) cycle
|
||||
|
||||
! add to cumulative probability
|
||||
yield = rx % products(i_product) % yield % evaluate(E)
|
||||
prob = prob + ((ONE - f) * rx % xs(i_temp) % value(i_grid - threshold + 1) &
|
||||
+ f*(rx % xs(i_temp) % value(i_grid - threshold + 2))) * yield
|
||||
|
||||
if (prob > cutoff) exit REACTION_LOOP
|
||||
end if
|
||||
end do
|
||||
end associate
|
||||
end do REACTION_LOOP
|
||||
|
||||
end subroutine sample_photon_product
|
||||
|
||||
!===============================================================================
|
||||
! ABSORPTION
|
||||
!===============================================================================
|
||||
|
|
@ -1498,4 +1554,50 @@ contains
|
|||
|
||||
end subroutine inelastic_scatter
|
||||
|
||||
!===============================================================================
|
||||
! SAMPLE_SECONDARY_PHOTONS
|
||||
!===============================================================================
|
||||
|
||||
subroutine sample_secondary_photons(p, i_nuclide)
|
||||
type(Particle), intent(inout) :: p
|
||||
integer, intent(in) :: i_nuclide
|
||||
|
||||
integer :: i_reaction ! index in nuc % reactions array
|
||||
integer :: i_product ! index in nuc % reactions % products array
|
||||
type(Reaction), pointer :: rx
|
||||
|
||||
real(8) :: nu_t
|
||||
real(8) :: mu
|
||||
real(8) :: E
|
||||
real(8) :: uvw(3)
|
||||
integer :: nu
|
||||
integer :: i
|
||||
|
||||
! Sample the number of photons produced
|
||||
nu_t = micro_xs(i_nuclide) % nu_photon_total / micro_xs(i_nuclide) % total
|
||||
if (prn() > nu_t - int(nu_t)) then
|
||||
nu = int(nu_t)
|
||||
else
|
||||
nu = int(nu_t) + 1
|
||||
end if
|
||||
|
||||
! Sample each secondary photon
|
||||
do i = 1, nu
|
||||
|
||||
! Sample the reaction and product
|
||||
call sample_photon_product(i_nuclide, p % E, i_reaction, i_product)
|
||||
rx => nuclides(i_nuclide) % reactions(i_reaction)
|
||||
|
||||
! Sample the outgoing energy and angle
|
||||
call rx % products(i_product) % sample(p % E, E, mu)
|
||||
|
||||
! Sample the new direction
|
||||
uvw = rotate_angle(p % coord(1) % uvw, mu)
|
||||
|
||||
! Create the secondary photon
|
||||
!call p % create_secondary(uvw, E, PHOTON, run_CE=.true.)
|
||||
end do
|
||||
|
||||
end subroutine sample_secondary_photons
|
||||
|
||||
end module physics
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue