mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 22:26:08 -04:00
Moved cross section calculation routines into separate module.
This commit is contained in:
Paul Romano
parent
7599ee7e31
commit
1602fe672a
4 changed files with 422 additions and 392 deletions
|
|
@ -14,6 +14,16 @@ ace.o: string.o
|
|||
ace_header.o: constants.o
|
||||
ace_header.o: endf_header.o
|
||||
|
||||
cross_section.o: ace_header.o
|
||||
cross_section.o: constants.o
|
||||
cross_section.o: error.o
|
||||
cross_section.o: fission.o
|
||||
cross_section.o: global.o
|
||||
cross_section.o: material_header.o
|
||||
cross_section.o: particle_header.o
|
||||
cross_section.o: random_lcg.o
|
||||
cross_section.o: search.o
|
||||
|
||||
datatypes.o: datatypes_header.o
|
||||
|
||||
doppler.o: constants.o
|
||||
|
|
@ -147,6 +157,7 @@ particle_header.o: geometry_header.o
|
|||
|
||||
physics.o: ace_header.o
|
||||
physics.o: constants.o
|
||||
physics.o: cross_section.o
|
||||
physics.o: endf.o
|
||||
physics.o: error.o
|
||||
physics.o: fission.o
|
||||
|
|
|
|||
|
|
@ -2,6 +2,7 @@ objects = \
|
|||
ace.o \
|
||||
ace_header.o \
|
||||
bank_header.o \
|
||||
cross_section.o \
|
||||
datatypes.o \
|
||||
datatypes_header.o \
|
||||
doppler.o \
|
||||
|
|
|
|||
408
src/cross_section.F90
Normal file
408
src/cross_section.F90
Normal file
|
|
@ -0,0 +1,408 @@
|
|||
module cross_section
|
||||
|
||||
use ace_header, only: Nuclide, SAB_Table, Reaction, UrrData
|
||||
use constants
|
||||
use error, only: fatal_error
|
||||
use fission, only: nu_total
|
||||
use global
|
||||
use material_header, only: Material
|
||||
use particle_header, only: Particle
|
||||
use random_lcg, only: prn
|
||||
use search, only: binary_search
|
||||
|
||||
implicit none
|
||||
|
||||
contains
|
||||
|
||||
!===============================================================================
|
||||
! CALCULATE_XS determines the macroscopic cross sections for the material the
|
||||
! particle is currently traveling through.
|
||||
!===============================================================================
|
||||
|
||||
subroutine calculate_xs(p)
|
||||
|
||||
type(Particle), pointer :: p
|
||||
|
||||
integer :: i ! loop index over nuclides
|
||||
integer :: index_nuclide ! index into nuclides array
|
||||
integer :: index_sab ! index into sab_tables array
|
||||
real(8) :: atom_density ! atom density of a nuclide
|
||||
real(8) :: sab_threshold ! threshold for S(a,b) table
|
||||
type(Material), pointer :: mat => null() ! current material
|
||||
|
||||
! 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
|
||||
|
||||
mat => materials(p % material)
|
||||
|
||||
! Find energy index on unionized grid
|
||||
call find_energy_index(p)
|
||||
|
||||
! Check if there's an S(a,b) table for this material
|
||||
if (mat % has_sab_table) then
|
||||
sab_threshold = sab_tables(mat % sab_table) % threshold_inelastic
|
||||
else
|
||||
sab_threshold = ZERO
|
||||
end if
|
||||
|
||||
! Add contribution from each nuclide in material
|
||||
do i = 1, mat % n_nuclides
|
||||
! Determine microscopic cross sections for this nuclide
|
||||
index_nuclide = mat % nuclide(i)
|
||||
|
||||
! Determine whether to use S(a,b) based on energy of particle and whether
|
||||
! the nuclide matches the S(a,b) table
|
||||
if (p % E < sab_threshold .and. i == mat % sab_nuclide) then
|
||||
index_sab = mat % sab_table
|
||||
else
|
||||
index_sab = 0
|
||||
end if
|
||||
|
||||
! Calculate microscopic cross section for this nuclide
|
||||
if (p % E /= micro_xs(index_nuclide) % last_E) then
|
||||
call calculate_nuclide_xs(p, index_nuclide, index_sab)
|
||||
end if
|
||||
|
||||
! 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_xs(index_nuclide) % total
|
||||
|
||||
! Add contributions to material macroscopic scattering cross section
|
||||
material_xs % elastic = material_xs % elastic + &
|
||||
atom_density * micro_xs(index_nuclide) % elastic
|
||||
|
||||
! Add contributions to material macroscopic absorption cross section
|
||||
material_xs % absorption = material_xs % absorption + &
|
||||
atom_density * micro_xs(index_nuclide) % absorption
|
||||
|
||||
! Add contributions to material macroscopic fission cross section
|
||||
material_xs % fission = material_xs % fission + &
|
||||
atom_density * micro_xs(index_nuclide) % fission
|
||||
|
||||
! Add contributions to material macroscopic nu-fission cross section
|
||||
material_xs % nu_fission = material_xs % nu_fission + &
|
||||
atom_density * micro_xs(index_nuclide) % nu_fission
|
||||
end do
|
||||
|
||||
end subroutine calculate_xs
|
||||
|
||||
!===============================================================================
|
||||
! CALCULATE_NUCLIDE_XS determines microscopic cross sections for a nuclide of a
|
||||
! given index in the nuclides array at the energy of the given particle
|
||||
!===============================================================================
|
||||
|
||||
subroutine calculate_nuclide_xs(p, index_nuclide, index_sab)
|
||||
|
||||
type(Particle), pointer :: p
|
||||
integer, intent(in) :: index_nuclide ! index into nuclides array
|
||||
integer, intent(in) :: index_sab ! index into sab_tables array
|
||||
|
||||
integer :: IE ! index on nuclide energy grid
|
||||
real(8) :: f ! interp factor on nuclide energy grid
|
||||
type(Nuclide), pointer :: nuc => null()
|
||||
|
||||
! Set pointer to nuclide
|
||||
nuc => nuclides(index_nuclide)
|
||||
|
||||
! TODO: If not using unionized energy grid, we need to find the index on the
|
||||
! nuclide energy grid using lethargy mapping or whatever other technique
|
||||
|
||||
! get index on nuclide energy grid
|
||||
IE = nuc % grid_index(p % IE)
|
||||
|
||||
! check for rare case where two energy points are the same
|
||||
if (nuc % energy(IE) == nuc % energy(IE+1)) IE = IE + 1
|
||||
|
||||
! calculate interpolation factor
|
||||
f = (p%E - nuc%energy(IE))/(nuc%energy(IE+1) - nuc%energy(IE))
|
||||
|
||||
micro_xs(index_nuclide) % index_grid = IE
|
||||
micro_xs(index_nuclide) % interp_factor = f
|
||||
|
||||
! Initialize sab treatment to false
|
||||
micro_xs(index_nuclide) % use_sab = .false.
|
||||
micro_xs(index_nuclide) % elastic_sab = ZERO
|
||||
|
||||
! Initialize nuclide cross-sections to zero
|
||||
micro_xs(index_nuclide) % fission = ZERO
|
||||
micro_xs(index_nuclide) % nu_fission = ZERO
|
||||
|
||||
! Calculate microscopic nuclide total cross section
|
||||
micro_xs(index_nuclide) % total = &
|
||||
(ONE-f) * nuc % total(IE) + f * nuc % total(IE+1)
|
||||
|
||||
! Calculate microscopic nuclide total cross section
|
||||
micro_xs(index_nuclide) % elastic = &
|
||||
(ONE-f) * nuc % elastic(IE) + f * nuc % elastic(IE+1)
|
||||
|
||||
! Calculate microscopic nuclide absorption cross section
|
||||
micro_xs(index_nuclide) % absorption = &
|
||||
(ONE-f) * nuc % absorption(IE) + f * nuc % absorption(IE+1)
|
||||
|
||||
if (nuc % fissionable) then
|
||||
! Calculate microscopic nuclide total cross section
|
||||
micro_xs(index_nuclide) % fission = &
|
||||
(ONE-f) * nuc % fission(IE) + f * nuc % fission(IE+1)
|
||||
|
||||
! Calculate microscopic nuclide nu-fission cross section
|
||||
micro_xs(index_nuclide) % nu_fission = &
|
||||
(ONE-f) * nuc % nu_fission(IE) + f * nuc % nu_fission(IE+1)
|
||||
end if
|
||||
|
||||
! If there is S(a,b) data for this nuclide, we need to do a few
|
||||
! things. Since the total cross section was based on non-S(a,b) data, we
|
||||
! need to correct it by subtracting the non-S(a,b) elastic cross section and
|
||||
! then add back in the calculated S(a,b) elastic+inelastic cross section.
|
||||
|
||||
if (index_sab > 0) call calculate_sab_xs(p, index_nuclide, index_sab)
|
||||
|
||||
! if the particle is in the unresolved resonance range and there are
|
||||
! probability tables, we need to determine cross sections from the table
|
||||
|
||||
if (urr_ptables_on .and. nuc % urr_present) then
|
||||
if (p % E > nuc % urr_data % energy(1) .and. &
|
||||
p % E < nuc % urr_data % energy(nuc % urr_data % n_energy)) then
|
||||
call calculate_urr_xs(p, index_nuclide)
|
||||
end if
|
||||
end if
|
||||
|
||||
! Set last evaluated energy
|
||||
micro_xs(index_nuclide) % last_E = p % E
|
||||
|
||||
end subroutine calculate_nuclide_xs
|
||||
|
||||
!===============================================================================
|
||||
! CALCULATE_SAB_XS determines the elastic and inelastic scattering
|
||||
! cross-sections in the thermal energy range. These cross sections replace
|
||||
! whatever data were taken from the normal Nuclide table.
|
||||
!===============================================================================
|
||||
|
||||
subroutine calculate_sab_xs(p, index_nuclide, index_sab)
|
||||
|
||||
type(Particle), pointer :: p
|
||||
integer, intent(in) :: index_nuclide ! index into nuclides array
|
||||
integer, intent(in) :: index_sab ! index into sab_tables array
|
||||
|
||||
integer :: IE_sab ! index on S(a,b) energy grid
|
||||
real(8) :: f_sab ! interp factor on S(a,b) energy grid
|
||||
real(8) :: inelastic ! S(a,b) inelastic cross section
|
||||
real(8) :: elastic ! S(a,b) elastic cross section
|
||||
type(SAB_Table), pointer :: sab => null()
|
||||
|
||||
! Set flag that S(a,b) treatment should be used for scattering
|
||||
micro_xs(index_nuclide) % use_sab = .true.
|
||||
|
||||
! Get pointer to S(a,b) table
|
||||
sab => sab_tables(index_sab)
|
||||
|
||||
! Get index and interpolation factor for inelastic grid
|
||||
if (p%E < sab % inelastic_e_in(1)) then
|
||||
IE_sab = 1
|
||||
f_sab = ZERO
|
||||
else
|
||||
IE_sab = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, p%E)
|
||||
f_sab = (p%E - sab%inelastic_e_in(IE_sab)) / &
|
||||
(sab%inelastic_e_in(IE_sab+1) - sab%inelastic_e_in(IE_sab))
|
||||
end if
|
||||
|
||||
! Calculate S(a,b) inelastic scattering cross section
|
||||
inelastic = (ONE-f_sab) * sab % inelastic_sigma(IE_sab) + f_sab * &
|
||||
sab % inelastic_sigma(IE_sab + 1)
|
||||
|
||||
! Check for elastic data
|
||||
if (p % E < sab % threshold_elastic) then
|
||||
! Determine whether elastic scattering is given in the coherent or
|
||||
! incoherent approximation. For coherent, the cross section is
|
||||
! represented as P/E whereas for incoherent, it is simply P
|
||||
|
||||
if (sab % elastic_mode == SAB_ELASTIC_EXACT) then
|
||||
if (p % E < sab % elastic_e_in(1)) then
|
||||
! If energy is below that of the lowest Bragg peak, the elastic
|
||||
! cross section will be zero
|
||||
elastic = ZERO
|
||||
else
|
||||
IE_sab = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E)
|
||||
elastic = sab % elastic_P(IE_sab) / p % E
|
||||
end if
|
||||
else
|
||||
! Determine index on elastic energy grid
|
||||
if (p % E < sab % elastic_e_in(1)) then
|
||||
IE_sab = 1
|
||||
else
|
||||
IE_sab = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E)
|
||||
end if
|
||||
|
||||
! Get interpolation factor for elastic grid
|
||||
f_sab = (p%E - sab%elastic_e_in(IE_sab))/(sab%elastic_e_in(IE_sab+1) - &
|
||||
sab%elastic_e_in(IE_sab))
|
||||
|
||||
! Calculate S(a,b) elastic scattering cross section
|
||||
elastic = (ONE-f_sab) * sab % elastic_P(IE_sab) + f_sab * &
|
||||
sab % elastic_P(IE_sab + 1)
|
||||
end if
|
||||
else
|
||||
! No elastic data
|
||||
elastic = ZERO
|
||||
end if
|
||||
|
||||
! Correct total and elastic cross sections
|
||||
micro_xs(index_nuclide) % total = micro_xs(index_nuclide) % total - &
|
||||
micro_xs(index_nuclide) % elastic + inelastic + elastic
|
||||
micro_xs(index_nuclide) % elastic = inelastic + elastic
|
||||
|
||||
! Store S(a,b) elastic cross section for sampling later
|
||||
micro_xs(index_nuclide) % elastic_sab = elastic
|
||||
|
||||
end subroutine calculate_sab_xs
|
||||
|
||||
!===============================================================================
|
||||
! CALCULATE_URR_XS determines cross sections in the unresolved resonance range
|
||||
! from probability tables
|
||||
!===============================================================================
|
||||
|
||||
subroutine calculate_urr_xs(p, index_nuclide)
|
||||
|
||||
type(Particle), pointer :: p
|
||||
integer, intent(in) :: index_nuclide ! index into nuclides array
|
||||
|
||||
integer :: i_energy ! index for energy
|
||||
integer :: i_table ! index for table
|
||||
real(8) :: f ! interpolation factor
|
||||
real(8) :: r ! pseudo-random number
|
||||
real(8) :: elastic ! smooth elastic cross section
|
||||
real(8) :: absorption ! smooth absorption cross section
|
||||
real(8) :: fission ! smooth fission cross section
|
||||
real(8) :: inelastic ! smooth inelastic cross section
|
||||
type(UrrData), pointer :: urr => null()
|
||||
type(Nuclide), pointer :: nuc => null()
|
||||
type(Reaction), pointer :: rxn => null()
|
||||
|
||||
! copy cross-sections already calculated
|
||||
elastic = micro_xs(index_nuclide) % elastic
|
||||
absorption = micro_xs(index_nuclide) % absorption
|
||||
fission = micro_xs(index_nuclide) % fission
|
||||
|
||||
! get pointer to probability table
|
||||
nuc => nuclides(index_nuclide)
|
||||
urr => nuc % urr_data
|
||||
|
||||
! determine energy table
|
||||
i_energy = 1
|
||||
do
|
||||
if (p % E > urr % energy(i_energy)) exit
|
||||
i_energy = i_energy + 1
|
||||
end do
|
||||
|
||||
! determine interpolation factor on table
|
||||
f = (p % E - urr % energy(i_energy)) / &
|
||||
(urr % energy(i_energy + 1) - urr % energy(i_energy))
|
||||
|
||||
! sample probability table using the cumulative distribution
|
||||
r = prn()
|
||||
i_table = 1
|
||||
do
|
||||
if (urr % prob(i_energy, URR_CUM_PROB, i_table) > r) exit
|
||||
i_table = i_table + 1
|
||||
end do
|
||||
|
||||
! determine elastic, fission, and capture cross sections from probability
|
||||
! table
|
||||
if (urr % interp == LINEAR_LINEAR) then
|
||||
micro_xs(index_nuclide) % elastic = &
|
||||
(ONE - f) * urr % prob(i_energy, URR_ELASTIC, i_table) + &
|
||||
f * urr % prob(i_energy + 1, URR_ELASTIC, i_table)
|
||||
micro_xs(index_nuclide) % fission = &
|
||||
(ONE - f) * urr % prob(i_energy, URR_FISSION, i_table) + &
|
||||
f * urr % prob(i_energy + 1, URR_FISSION, i_table)
|
||||
micro_xs(index_nuclide) % absorption = &
|
||||
micro_xs(index_nuclide) % fission + &
|
||||
(ONE - f) * urr % prob(i_energy, URR_N_GAMMA, i_table) + &
|
||||
f * urr % prob(i_energy + 1, URR_N_GAMMA, i_table)
|
||||
elseif (urr % interp == LOG_LOG) then
|
||||
message = "Log-log interpolation on probability table not yet supported."
|
||||
call fatal_error()
|
||||
end if
|
||||
|
||||
! Determine treatment of inelastic scattering
|
||||
inelastic = ZERO
|
||||
if (urr % inelastic_flag > 0) then
|
||||
! Get pointer to inelastic scattering reaction
|
||||
rxn => nuc % reactions(nuc % urr_inelastic)
|
||||
|
||||
! Get index on energy grid and interpolation factor
|
||||
i_energy = micro_xs(index_nuclide) % index_grid
|
||||
f = micro_xs(index_nuclide) % interp_factor
|
||||
|
||||
! Determine inelastic scattering cross section
|
||||
if (i_energy >= rxn % IE) then
|
||||
inelastic = (ONE - f) * rxn % sigma(i_energy - rxn%IE + 1) + &
|
||||
f * rxn % sigma(i_energy - rxn%IE + 2)
|
||||
end if
|
||||
end if
|
||||
|
||||
! Multiply by smooth cross-section if needed
|
||||
if (urr % multiply_smooth) then
|
||||
micro_xs(index_nuclide) % elastic = elastic * &
|
||||
micro_xs(index_nuclide) % elastic
|
||||
micro_xs(index_nuclide) % absorption = absorption * &
|
||||
micro_xs(index_nuclide) % absorption
|
||||
micro_xs(index_nuclide) % fission = fission * &
|
||||
micro_xs(index_nuclide) % fission
|
||||
end if
|
||||
|
||||
! Determine nu-fission cross section
|
||||
if (nuc % fissionable) then
|
||||
micro_xs(index_nuclide) % nu_fission = nu_total(nuc, p % E) * &
|
||||
micro_xs(index_nuclide) % fission
|
||||
end if
|
||||
|
||||
! Calculate total cross section as sum of partials
|
||||
micro_xs(index_nuclide) % total = micro_xs(index_nuclide) % elastic + &
|
||||
inelastic + micro_xs(index_nuclide) % absorption
|
||||
|
||||
end subroutine calculate_urr_xs
|
||||
|
||||
!===============================================================================
|
||||
! FIND_ENERGY_INDEX determines the index on the union energy grid and the
|
||||
! interpolation factor for a particle at a certain energy
|
||||
!===============================================================================
|
||||
|
||||
subroutine find_energy_index(p)
|
||||
|
||||
type(Particle), pointer :: p
|
||||
|
||||
integer :: IE ! index on union energy grid
|
||||
real(8) :: E ! energy of particle
|
||||
real(8) :: interp ! interpolation factor
|
||||
|
||||
! copy particle's energy
|
||||
E = p % E
|
||||
|
||||
! if particle's energy is outside of energy grid range, set to first or last
|
||||
! index. Otherwise, do a binary search through the union energy grid.
|
||||
if (E < e_grid(1)) then
|
||||
IE = 1
|
||||
elseif (E > e_grid(n_grid)) then
|
||||
IE = n_grid - 1
|
||||
else
|
||||
IE = binary_search(e_grid, n_grid, E)
|
||||
end if
|
||||
|
||||
! calculate the interpolation factor -- note this will be outside of [0,1)
|
||||
! for a particle outside the energy range of the union grid
|
||||
interp = (E - e_grid(IE))/(e_grid(IE+1) - e_grid(IE))
|
||||
|
||||
! set particle attributes
|
||||
p % IE = IE
|
||||
p % interp = interp
|
||||
|
||||
end subroutine find_energy_index
|
||||
|
||||
end module cross_section
|
||||
394
src/physics.F90
394
src/physics.F90
|
|
@ -1,7 +1,8 @@
|
|||
module physics
|
||||
|
||||
use ace_header, only: Nuclide, Reaction, DistEnergy, UrrData
|
||||
use ace_header, only: Nuclide, Reaction, DistEnergy
|
||||
use constants
|
||||
use cross_section, only: calculate_xs, find_energy_index
|
||||
use endf, only: reaction_name, is_fission, is_scatter
|
||||
use error, only: fatal_error, warning
|
||||
use fission, only: nu_total, nu_prompt, nu_delayed
|
||||
|
|
@ -135,397 +136,6 @@ contains
|
|||
|
||||
end subroutine transport
|
||||
|
||||
!===============================================================================
|
||||
! CALCULATE_XS determines the macroscopic cross sections for the material the
|
||||
! particle is currently traveling through.
|
||||
!===============================================================================
|
||||
|
||||
subroutine calculate_xs(p)
|
||||
|
||||
type(Particle), pointer :: p
|
||||
|
||||
integer :: i ! loop index over nuclides
|
||||
integer :: index_nuclide ! index into nuclides array
|
||||
integer :: index_sab ! index into sab_tables array
|
||||
real(8) :: atom_density ! atom density of a nuclide
|
||||
real(8) :: sab_threshold ! threshold for S(a,b) table
|
||||
type(Material), pointer :: mat => null() ! current material
|
||||
|
||||
! 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
|
||||
|
||||
mat => materials(p % material)
|
||||
|
||||
! Find energy index on unionized grid
|
||||
call find_energy_index(p)
|
||||
|
||||
! Check if there's an S(a,b) table for this material
|
||||
if (mat % has_sab_table) then
|
||||
sab_threshold = sab_tables(mat % sab_table) % threshold_inelastic
|
||||
else
|
||||
sab_threshold = ZERO
|
||||
end if
|
||||
|
||||
! Add contribution from each nuclide in material
|
||||
do i = 1, mat % n_nuclides
|
||||
! Determine microscopic cross sections for this nuclide
|
||||
index_nuclide = mat % nuclide(i)
|
||||
|
||||
! Determine whether to use S(a,b) based on energy of particle and whether
|
||||
! the nuclide matches the S(a,b) table
|
||||
if (p % E < sab_threshold .and. i == mat % sab_nuclide) then
|
||||
index_sab = mat % sab_table
|
||||
else
|
||||
index_sab = 0
|
||||
end if
|
||||
|
||||
! Calculate microscopic cross section for this nuclide
|
||||
if (p % E /= micro_xs(index_nuclide) % last_E) then
|
||||
call calculate_nuclide_xs(p, index_nuclide, index_sab)
|
||||
end if
|
||||
|
||||
! 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_xs(index_nuclide) % total
|
||||
|
||||
! Add contributions to material macroscopic scattering cross section
|
||||
material_xs % elastic = material_xs % elastic + &
|
||||
atom_density * micro_xs(index_nuclide) % elastic
|
||||
|
||||
! Add contributions to material macroscopic absorption cross section
|
||||
material_xs % absorption = material_xs % absorption + &
|
||||
atom_density * micro_xs(index_nuclide) % absorption
|
||||
|
||||
! Add contributions to material macroscopic fission cross section
|
||||
material_xs % fission = material_xs % fission + &
|
||||
atom_density * micro_xs(index_nuclide) % fission
|
||||
|
||||
! Add contributions to material macroscopic nu-fission cross section
|
||||
material_xs % nu_fission = material_xs % nu_fission + &
|
||||
atom_density * micro_xs(index_nuclide) % nu_fission
|
||||
end do
|
||||
|
||||
end subroutine calculate_xs
|
||||
|
||||
!===============================================================================
|
||||
! CALCULATE_NUCLIDE_XS determines microscopic cross sections for a nuclide of a
|
||||
! given index in the nuclides array at the energy of the given particle
|
||||
!===============================================================================
|
||||
|
||||
subroutine calculate_nuclide_xs(p, index_nuclide, index_sab)
|
||||
|
||||
type(Particle), pointer :: p
|
||||
integer, intent(in) :: index_nuclide ! index into nuclides array
|
||||
integer, intent(in) :: index_sab ! index into sab_tables array
|
||||
|
||||
integer :: IE ! index on nuclide energy grid
|
||||
real(8) :: f ! interp factor on nuclide energy grid
|
||||
type(Nuclide), pointer :: nuc => null()
|
||||
|
||||
! Set pointer to nuclide
|
||||
nuc => nuclides(index_nuclide)
|
||||
|
||||
! TODO: If not using unionized energy grid, we need to find the index on the
|
||||
! nuclide energy grid using lethargy mapping or whatever other technique
|
||||
|
||||
! get index on nuclide energy grid
|
||||
IE = nuc % grid_index(p % IE)
|
||||
|
||||
! check for rare case where two energy points are the same
|
||||
if (nuc % energy(IE) == nuc % energy(IE+1)) IE = IE + 1
|
||||
|
||||
! calculate interpolation factor
|
||||
f = (p%E - nuc%energy(IE))/(nuc%energy(IE+1) - nuc%energy(IE))
|
||||
|
||||
micro_xs(index_nuclide) % index_grid = IE
|
||||
micro_xs(index_nuclide) % interp_factor = f
|
||||
|
||||
! Initialize sab treatment to false
|
||||
micro_xs(index_nuclide) % use_sab = .false.
|
||||
micro_xs(index_nuclide) % elastic_sab = ZERO
|
||||
|
||||
! Initialize nuclide cross-sections to zero
|
||||
micro_xs(index_nuclide) % fission = ZERO
|
||||
micro_xs(index_nuclide) % nu_fission = ZERO
|
||||
|
||||
! Calculate microscopic nuclide total cross section
|
||||
micro_xs(index_nuclide) % total = &
|
||||
(ONE-f) * nuc % total(IE) + f * nuc % total(IE+1)
|
||||
|
||||
! Calculate microscopic nuclide total cross section
|
||||
micro_xs(index_nuclide) % elastic = &
|
||||
(ONE-f) * nuc % elastic(IE) + f * nuc % elastic(IE+1)
|
||||
|
||||
! Calculate microscopic nuclide absorption cross section
|
||||
micro_xs(index_nuclide) % absorption = &
|
||||
(ONE-f) * nuc % absorption(IE) + f * nuc % absorption(IE+1)
|
||||
|
||||
if (nuc % fissionable) then
|
||||
! Calculate microscopic nuclide total cross section
|
||||
micro_xs(index_nuclide) % fission = &
|
||||
(ONE-f) * nuc % fission(IE) + f * nuc % fission(IE+1)
|
||||
|
||||
! Calculate microscopic nuclide nu-fission cross section
|
||||
micro_xs(index_nuclide) % nu_fission = &
|
||||
(ONE-f) * nuc % nu_fission(IE) + f * nuc % nu_fission(IE+1)
|
||||
end if
|
||||
|
||||
! If there is S(a,b) data for this nuclide, we need to do a few
|
||||
! things. Since the total cross section was based on non-S(a,b) data, we
|
||||
! need to correct it by subtracting the non-S(a,b) elastic cross section and
|
||||
! then add back in the calculated S(a,b) elastic+inelastic cross section.
|
||||
|
||||
if (index_sab > 0) call calculate_sab_xs(p, index_nuclide, index_sab)
|
||||
|
||||
! if the particle is in the unresolved resonance range and there are
|
||||
! probability tables, we need to determine cross sections from the table
|
||||
|
||||
if (urr_ptables_on .and. nuc % urr_present) then
|
||||
if (p % E > nuc % urr_data % energy(1) .and. &
|
||||
p % E < nuc % urr_data % energy(nuc % urr_data % n_energy)) then
|
||||
call calculate_urr_xs(p, index_nuclide)
|
||||
end if
|
||||
end if
|
||||
|
||||
! Set last evaluated energy
|
||||
micro_xs(index_nuclide) % last_E = p % E
|
||||
|
||||
end subroutine calculate_nuclide_xs
|
||||
|
||||
!===============================================================================
|
||||
! CALCULATE_SAB_XS determines the elastic and inelastic scattering
|
||||
! cross-sections in the thermal energy range. These cross sections replace
|
||||
! whatever data were taken from the normal Nuclide table.
|
||||
!===============================================================================
|
||||
|
||||
subroutine calculate_sab_xs(p, index_nuclide, index_sab)
|
||||
|
||||
type(Particle), pointer :: p
|
||||
integer, intent(in) :: index_nuclide ! index into nuclides array
|
||||
integer, intent(in) :: index_sab ! index into sab_tables array
|
||||
|
||||
integer :: IE_sab ! index on S(a,b) energy grid
|
||||
real(8) :: f_sab ! interp factor on S(a,b) energy grid
|
||||
real(8) :: inelastic ! S(a,b) inelastic cross section
|
||||
real(8) :: elastic ! S(a,b) elastic cross section
|
||||
type(SAB_Table), pointer :: sab => null()
|
||||
|
||||
! Set flag that S(a,b) treatment should be used for scattering
|
||||
micro_xs(index_nuclide) % use_sab = .true.
|
||||
|
||||
! Get pointer to S(a,b) table
|
||||
sab => sab_tables(index_sab)
|
||||
|
||||
! Get index and interpolation factor for inelastic grid
|
||||
if (p%E < sab % inelastic_e_in(1)) then
|
||||
IE_sab = 1
|
||||
f_sab = ZERO
|
||||
else
|
||||
IE_sab = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, p%E)
|
||||
f_sab = (p%E - sab%inelastic_e_in(IE_sab)) / &
|
||||
(sab%inelastic_e_in(IE_sab+1) - sab%inelastic_e_in(IE_sab))
|
||||
end if
|
||||
|
||||
! Calculate S(a,b) inelastic scattering cross section
|
||||
inelastic = (ONE-f_sab) * sab % inelastic_sigma(IE_sab) + f_sab * &
|
||||
sab % inelastic_sigma(IE_sab + 1)
|
||||
|
||||
! Check for elastic data
|
||||
if (p % E < sab % threshold_elastic) then
|
||||
! Determine whether elastic scattering is given in the coherent or
|
||||
! incoherent approximation. For coherent, the cross section is
|
||||
! represented as P/E whereas for incoherent, it is simply P
|
||||
|
||||
if (sab % elastic_mode == SAB_ELASTIC_EXACT) then
|
||||
if (p % E < sab % elastic_e_in(1)) then
|
||||
! If energy is below that of the lowest Bragg peak, the elastic
|
||||
! cross section will be zero
|
||||
elastic = ZERO
|
||||
else
|
||||
IE_sab = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E)
|
||||
elastic = sab % elastic_P(IE_sab) / p % E
|
||||
end if
|
||||
else
|
||||
! Determine index on elastic energy grid
|
||||
if (p % E < sab % elastic_e_in(1)) then
|
||||
IE_sab = 1
|
||||
else
|
||||
IE_sab = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E)
|
||||
end if
|
||||
|
||||
! Get interpolation factor for elastic grid
|
||||
f_sab = (p%E - sab%elastic_e_in(IE_sab))/(sab%elastic_e_in(IE_sab+1) - &
|
||||
sab%elastic_e_in(IE_sab))
|
||||
|
||||
! Calculate S(a,b) elastic scattering cross section
|
||||
elastic = (ONE-f_sab) * sab % elastic_P(IE_sab) + f_sab * &
|
||||
sab % elastic_P(IE_sab + 1)
|
||||
end if
|
||||
else
|
||||
! No elastic data
|
||||
elastic = ZERO
|
||||
end if
|
||||
|
||||
! Correct total and elastic cross sections
|
||||
micro_xs(index_nuclide) % total = micro_xs(index_nuclide) % total - &
|
||||
micro_xs(index_nuclide) % elastic + inelastic + elastic
|
||||
micro_xs(index_nuclide) % elastic = inelastic + elastic
|
||||
|
||||
! Store S(a,b) elastic cross section for sampling later
|
||||
micro_xs(index_nuclide) % elastic_sab = elastic
|
||||
|
||||
end subroutine calculate_sab_xs
|
||||
|
||||
!===============================================================================
|
||||
! CALCULATE_URR_XS determines cross sections in the unresolved resonance range
|
||||
! from probability tables
|
||||
!===============================================================================
|
||||
|
||||
subroutine calculate_urr_xs(p, index_nuclide)
|
||||
|
||||
type(Particle), pointer :: p
|
||||
integer, intent(in) :: index_nuclide ! index into nuclides array
|
||||
|
||||
integer :: i_energy ! index for energy
|
||||
integer :: i_table ! index for table
|
||||
real(8) :: f ! interpolation factor
|
||||
real(8) :: r ! pseudo-random number
|
||||
real(8) :: elastic ! smooth elastic cross section
|
||||
real(8) :: absorption ! smooth absorption cross section
|
||||
real(8) :: fission ! smooth fission cross section
|
||||
real(8) :: inelastic ! smooth inelastic cross section
|
||||
type(UrrData), pointer :: urr => null()
|
||||
type(Nuclide), pointer :: nuc => null()
|
||||
type(Reaction), pointer :: rxn => null()
|
||||
|
||||
! copy cross-sections already calculated
|
||||
elastic = micro_xs(index_nuclide) % elastic
|
||||
absorption = micro_xs(index_nuclide) % absorption
|
||||
fission = micro_xs(index_nuclide) % fission
|
||||
|
||||
! get pointer to probability table
|
||||
nuc => nuclides(index_nuclide)
|
||||
urr => nuc % urr_data
|
||||
|
||||
! determine energy table
|
||||
i_energy = 1
|
||||
do
|
||||
if (p % E > urr % energy(i_energy)) exit
|
||||
i_energy = i_energy + 1
|
||||
end do
|
||||
|
||||
! determine interpolation factor on table
|
||||
f = (p % E - urr % energy(i_energy)) / &
|
||||
(urr % energy(i_energy + 1) - urr % energy(i_energy))
|
||||
|
||||
! sample probability table using the cumulative distribution
|
||||
r = prn()
|
||||
i_table = 1
|
||||
do
|
||||
if (urr % prob(i_energy, URR_CUM_PROB, i_table) > r) exit
|
||||
i_table = i_table + 1
|
||||
end do
|
||||
|
||||
! determine elastic, fission, and capture cross sections from probability
|
||||
! table
|
||||
if (urr % interp == LINEAR_LINEAR) then
|
||||
micro_xs(index_nuclide) % elastic = &
|
||||
(ONE - f) * urr % prob(i_energy, URR_ELASTIC, i_table) + &
|
||||
f * urr % prob(i_energy + 1, URR_ELASTIC, i_table)
|
||||
micro_xs(index_nuclide) % fission = &
|
||||
(ONE - f) * urr % prob(i_energy, URR_FISSION, i_table) + &
|
||||
f * urr % prob(i_energy + 1, URR_FISSION, i_table)
|
||||
micro_xs(index_nuclide) % absorption = &
|
||||
micro_xs(index_nuclide) % fission + &
|
||||
(ONE - f) * urr % prob(i_energy, URR_N_GAMMA, i_table) + &
|
||||
f * urr % prob(i_energy + 1, URR_N_GAMMA, i_table)
|
||||
elseif (urr % interp == LOG_LOG) then
|
||||
message = "Log-log interpolation on probability table not yet supported."
|
||||
call fatal_error()
|
||||
end if
|
||||
|
||||
! Determine treatment of inelastic scattering
|
||||
inelastic = ZERO
|
||||
if (urr % inelastic_flag > 0) then
|
||||
! Get pointer to inelastic scattering reaction
|
||||
rxn => nuc % reactions(nuc % urr_inelastic)
|
||||
|
||||
! Get index on energy grid and interpolation factor
|
||||
i_energy = micro_xs(index_nuclide) % index_grid
|
||||
f = micro_xs(index_nuclide) % interp_factor
|
||||
|
||||
! Determine inelastic scattering cross section
|
||||
if (i_energy >= rxn % IE) then
|
||||
inelastic = (ONE - f) * rxn % sigma(i_energy - rxn%IE + 1) + &
|
||||
f * rxn % sigma(i_energy - rxn%IE + 2)
|
||||
end if
|
||||
end if
|
||||
|
||||
! Multiply by smooth cross-section if needed
|
||||
if (urr % multiply_smooth) then
|
||||
micro_xs(index_nuclide) % elastic = elastic * &
|
||||
micro_xs(index_nuclide) % elastic
|
||||
micro_xs(index_nuclide) % absorption = absorption * &
|
||||
micro_xs(index_nuclide) % absorption
|
||||
micro_xs(index_nuclide) % fission = fission * &
|
||||
micro_xs(index_nuclide) % fission
|
||||
end if
|
||||
|
||||
! Determine nu-fission cross section
|
||||
if (nuc % fissionable) then
|
||||
micro_xs(index_nuclide) % nu_fission = nu_total(nuc, p % E) * &
|
||||
micro_xs(index_nuclide) % fission
|
||||
end if
|
||||
|
||||
! Calculate total cross section as sum of partials
|
||||
micro_xs(index_nuclide) % total = micro_xs(index_nuclide) % elastic + &
|
||||
inelastic + micro_xs(index_nuclide) % absorption
|
||||
|
||||
end subroutine calculate_urr_xs
|
||||
|
||||
!===============================================================================
|
||||
! FIND_ENERGY_INDEX determines the index on the union energy grid and the
|
||||
! interpolation factor for a particle at a certain energy
|
||||
!===============================================================================
|
||||
|
||||
subroutine find_energy_index(p)
|
||||
|
||||
type(Particle), pointer :: p
|
||||
|
||||
integer :: IE ! index on union energy grid
|
||||
real(8) :: E ! energy of particle
|
||||
real(8) :: interp ! interpolation factor
|
||||
|
||||
! copy particle's energy
|
||||
E = p % E
|
||||
|
||||
! if particle's energy is outside of energy grid range, set to first or last
|
||||
! index. Otherwise, do a binary search through the union energy grid.
|
||||
if (E < e_grid(1)) then
|
||||
IE = 1
|
||||
elseif (E > e_grid(n_grid)) then
|
||||
IE = n_grid - 1
|
||||
else
|
||||
IE = binary_search(e_grid, n_grid, E)
|
||||
end if
|
||||
|
||||
! calculate the interpolation factor -- note this will be outside of [0,1)
|
||||
! for a particle outside the energy range of the union grid
|
||||
interp = (E - e_grid(IE))/(e_grid(IE+1) - e_grid(IE))
|
||||
|
||||
! set particle attributes
|
||||
p % IE = IE
|
||||
p % interp = interp
|
||||
|
||||
end subroutine find_energy_index
|
||||
|
||||
!===============================================================================
|
||||
! COLLISION samples a nuclide and reaction and then calls the appropriate
|
||||
! routine for that reaction
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue