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

Added cross section caches and changed physics module substantially. Moved calculation of nu into a new module named fission.

Browse source
This commit is contained in:
Paul Romano 2011-09-13 20:41:36 -04:00
parent 7392cc44b2
commit 3404aafaac
11 changed files with 598 additions and 316 deletions
Download patch file Download diff file Expand all files Collapse all files

6
src/DEPENDENCIES
View file

@ -29,6 +29,11 @@ fileio.o: constants.o
fileio.o: global.o
fileio.o: string.o
fission.o: constants.o
fission.o: cross_section_header.o
fission.o: error.o
fission.o: search.o
geometry.o: constants.o
geometry.o: datatypes.o
geometry.o: error.o
@ -113,6 +118,7 @@ physics.o: constants.o
physics.o: cross_section_header.o
physics.o: endf.o
physics.o: error.o
physics.o: fission.o
physics.o: geometry.o
physics.o: geometry_header.o
physics.o: global.o

1
src/OBJECTS
View file

@ -9,6 +9,7 @@ endf.o \
energy_grid.o \
error.o \
fileio.o \
fission.o \
geometry.o \
geometry_header.o \
global.o \

104
src/constants.f90
View file

@ -2,6 +2,11 @@ module constants
implicit none
! Versioning numbers
integer, parameter :: VERSION_MAJOR = 0
integer, parameter :: VERSION_MINOR = 3
integer, parameter :: VERSION_RELEASE = 2
! Physical constants
real(8), parameter :: &
PI = 3.1415926535898_8, & ! pi
@ -98,53 +103,53 @@ module constants
! Reaction types
integer, parameter :: &
TOTAL_XS = 1, &
ELASTIC = 2, &
N_LEVEL = 4, &
MISC = 5, &
N_2ND = 11, &
N_2N = 16, &
N_3N = 17, &
FISSION = 18, &
N_F = 19, &
N_NF = 20, &
N_2NF = 21, &
N_NA = 22, &
N_N3A = 23, &
N_2NA = 24, &
N_3NA = 25, &
N_NP = 28, &
N_N2A = 29, &
N_2N2A = 30, &
N_ND = 32, &
N_NT = 33, &
N_N3HE = 34, &
N_ND2A = 35, &
N_NT2A = 36, &
N_4N = 37, &
N_3NF = 38, &
N_2NP = 41, &
N_3NP = 42, &
N_N2P = 44, &
N_NPA = 45, &
N_N1 = 51, &
N_N40 = 90, &
N_NC = 91, &
N_GAMMA = 102, &
N_P = 103, &
N_D = 104, &
N_T = 105, &
N_3HE = 106, &
N_A = 107, &
N_2A = 108, &
N_3A = 109, &
N_2P = 111, &
N_PA = 112, &
N_T2A = 113, &
N_D2A = 114, &
N_PD = 115, &
N_PT = 116, &
N_DA = 117
TOTAL_XS = 1, &
ELASTIC = 2, &
N_LEVEL = 4, &
MISC = 5, &
N_2ND = 11, &
N_2N = 16, &
N_3N = 17, &
N_FISSION = 18, &
N_F = 19, &
N_NF = 20, &
N_2NF = 21, &
N_NA = 22, &
N_N3A = 23, &
N_2NA = 24, &
N_3NA = 25, &
N_NP = 28, &
N_N2A = 29, &
N_2N2A = 30, &
N_ND = 32, &
N_NT = 33, &
N_N3HE = 34, &
N_ND2A = 35, &
N_NT2A = 36, &
N_4N = 37, &
N_3NF = 38, &
N_2NP = 41, &
N_3NP = 42, &
N_N2P = 44, &
N_NPA = 45, &
N_N1 = 51, &
N_N40 = 90, &
N_NC = 91, &
N_GAMMA = 102, &
N_P = 103, &
N_D = 104, &
N_T = 105, &
N_3HE = 106, &
N_A = 107, &
N_2A = 108, &
N_3A = 109, &
N_2P = 111, &
N_PA = 112, &
N_T2A = 113, &
N_D2A = 114, &
N_PD = 115, &
N_PT = 116, &
N_DA = 117
! Tally distinguishers
integer, parameter :: &
@ -159,11 +164,6 @@ module constants
NU_POLYNOMIAL = 1, & ! Nu values given by polynomial
NU_TABULAR = 2 ! Nu values given by tabular distribution
! Versioning numbers
integer, parameter :: VERSION_MAJOR = 0
integer, parameter :: VERSION_MINOR = 3
integer, parameter :: VERSION_RELEASE = 1
! Maximum number of words in a single line, length of line, and length of
! single word
integer, parameter :: MAX_WORDS = 500

115
src/cross_section.f90
View file

@ -84,6 +84,9 @@ contains
allocate(nuclides(n_nuclides_total))
allocate(sab_tables(n_sab_tables))
! allocate array for microscopic cross section cache
allocate(micro_xs(n_nuclides_total))
! loop over all nuclides in xsdata
call dict_create(temp_dict)
@ -241,25 +244,36 @@ contains
! determine number of energy points
NE = NXS(3)
nuc%n_grid = NE
nuc % n_grid = NE
! allocate storage for arrays
allocate(nuc%energy(NE))
allocate(nuc%sigma_t(NE))
allocate(nuc%sigma_a(NE))
allocate(nuc%sigma_el(NE))
allocate(nuc%heating(NE))
! allocate storage for energy grid and cross section arrays
allocate(nuc % energy(NE))
allocate(nuc % total(NE))
allocate(nuc % elastic(NE))
allocate(nuc % fission(NE))
allocate(nuc % absorption(NE))
allocate(nuc % heating(NE))
! initialize cross sections
nuc % total = ZERO
nuc % elastic = ZERO
nuc % fission = ZERO
nuc % absorption = ZERO
nuc % heating = ZERO
! Read data from XSS -- only the energy grid, elastic scattering and heating
! cross section values are actually read from here. The total and absorption
! cross sections are reconstructed from the partial reaction data.
! read data from XSS -- right now the total, absorption and elastic
! scattering are read in to these special arrays, but in reality, it should
! be necessary to only store elastic scattering and possibly total
! cross-section for total material xs generation.
XSS_index = 1
nuc%energy = get_real(NE)
nuc%sigma_t = get_real(NE)
nuc%sigma_a = get_real(NE)
nuc%sigma_el = get_real(NE)
nuc%heating = get_real(NE)
nuc % energy = get_real(NE)
! Skip total and absorption
XSS_index = XSS_index + 2*NE
! Continue reading elastic scattering and heating
nuc % elastic = get_real(NE)
nuc % heating = get_real(NE)
end subroutine read_esz
@ -515,6 +529,7 @@ contains
integer :: JXS7 ! index of reactions cross-sections in XSS
integer :: LXS ! location of cross-section locators
integer :: LOCA ! location of cross-section for given MT
integer :: IE ! reaction's starting index on energy grid
integer :: NE ! number of energies for reaction
type(Reaction), pointer :: rxn => null()
@ -527,34 +542,60 @@ contains
! allocate array of reactions. Add one since we need to include an elastic
! scattering channel
nuc%n_reaction = NMT + 1
allocate(nuc%reactions(NMT+1))
nuc % n_reaction = NMT + 1
allocate(nuc % reactions(NMT+1))
! Store elastic scattering cross-section on reaction one
rxn => nuc%reactions(1)
rxn%MT = 2
rxn%Q_value = ZERO
rxn%TY = 1
rxn%IE = 1
allocate(rxn%sigma(nuc%n_grid))
rxn%sigma = nuc%sigma_el
rxn => nuc % reactions(1)
rxn % MT = 2
rxn % Q_value = ZERO
rxn % TY = 1
rxn % IE = 1
allocate(rxn % sigma(nuc % n_grid))
rxn % sigma = nuc % elastic
! Add contribution of elastic scattering to total cross section
nuc % total = nuc % total + nuc % elastic
do i = 1, NMT
rxn => nuc%reactions(i+1)
rxn => nuc % reactions(i+1)
! read MT number, Q-value, and neutrons produced
rxn % MT = XSS(LMT+i-1)
rxn % Q_value = XSS(JXS4+i-1)
rxn % TY = XSS(JXS5+i-1)
rxn % MT = XSS(LMT + i - 1)
rxn % Q_value = XSS(JXS4 + i - 1)
rxn % TY = XSS(JXS5 + i - 1)
! read cross section values
LOCA = XSS(LXS+i-1)
rxn % IE = XSS(JXS7 + LOCA - 1)
! read starting energy index
LOCA = XSS(LXS + i - 1)
IE = XSS(JXS7 + LOCA - 1)
rxn % IE = IE
! read number of energies cross section values
NE = XSS(JXS7 + LOCA)
allocate(rxn%sigma(NE))
allocate(rxn % sigma(NE))
XSS_index = JXS7 + LOCA + 1
rxn % sigma = get_real(NE)
! Skip redundant reactions -- this includes total inelastic level
! scattering, gas production cross sections (MT=200+), and (n,p), (n,d),
! etc. reactions leaving the nucleus in an excited state
if (rxn % MT == N_LEVEL .or. rxn % MT > N_DA) cycle
! Add contribution to total cross section
nuc % total(IE:IE+NE-1) = nuc % total(IE:IE+NE-1) + rxn % sigma
! Add contribution to absorption cross section
if (rxn % MT >= N_GAMMA .and. rxn % MT <= N_DA) then
nuc % absorption(IE:IE+NE-1) = nuc % absorption(IE:IE+NE-1) + rxn % sigma
end if
! Add contribution to fission cross section
if (rxn % MT == N_FISSION .or. rxn % MT == N_F .or. rxn % MT == N_NF &
.or. rxn % MT == N_2NF .or. rxn % MT == N_3NF) then
nuc % fissionable = .true.
nuc % fission(IE:IE+NE-1) = nuc % fission(IE:IE+NE-1) + rxn % sigma
end if
! set defaults
rxn % has_angle_dist = .false.
rxn % has_energy_dist = .false.
@ -1197,8 +1238,8 @@ contains
! handle special case of total cross section
if (MT == 1) then
xs = xs + mat % density * (ONE-f) * nuc%sigma_t(IE) + &
& f * (nuc%sigma_t(IE+1))
xs = xs + mat % density * (ONE-f) * nuc%total(IE) + &
& f * (nuc%total(IE+1))
cycle
end if
@ -1376,12 +1417,12 @@ contains
! get nuclide grid energy and cross-section value
E_i = nuc % energy(IE)
sigma_i = nuc % sigma_t(IE)
sigma_i = nuc % total(IE)
! interpolate as necessary
if (E_i < e_grid(k)) then
E_i1 = nuc % energy(IE + 1)
sigma_i1 = nuc % sigma_t(IE + 1)
sigma_i1 = nuc % total(IE + 1)
r = (e_grid(k) - E_i)/(E_i1 - E_i)
val = (ONE - r)*sigma_i + r*sigma_i1

48
src/cross_section_header.f90
View file

@ -56,14 +56,20 @@ module cross_section_header
type Nuclide
character(20) :: name
real(8) :: awr
real(8) :: temp
real(8) :: awr
real(8) :: temp
logical :: fissionable
! Energy grid information
integer :: n_grid
integer, allocatable :: grid_index(:)
real(8), allocatable :: energy(:)
real(8), allocatable :: sigma_t(:)
real(8), allocatable :: sigma_a(:)
real(8), allocatable :: sigma_el(:)
! Cross sections
real(8), allocatable :: total(:)
real(8), allocatable :: elastic(:)
real(8), allocatable :: fission(:)
real(8), allocatable :: absorption(:)
real(8), allocatable :: heating(:)
! Total fission neutron emission
@ -133,4 +139,36 @@ module cross_section_header
character(150) :: path
end type xsData
!===============================================================================
! NUCLIDEMICROXS contains cached microscopic cross sections for a
! particular nuclide at the current energy
!===============================================================================
type NuclideMicroXS
integer :: index_grid
integer :: index_temp
integer :: last_index_grid
integer :: last_index_temp
real(8) :: interp_factor
real(8) :: total
real(8) :: elastic
real(8) :: absorption
real(8) :: fission
real(8) :: nu_fission
end type NuclideMicroXS
!===============================================================================
! MATERIALMACROXS contains cached macroscopic cross sections for the material a
! particle is traveling through
!===============================================================================
type MaterialMacroXS
real(8) :: total
real(8) :: scatter
real(8) :: elastic
real(8) :: absorption
real(8) :: fission
real(8) :: nu_fission
end type MaterialMacroXS
end module cross_section_header

2
src/endf.f90
View file

@ -27,7 +27,7 @@ contains
string = '(n,2n)'
case (N_3N)
string = '(n,3n)'
case (FISSION)
case (N_FISSION)
string = '(n,fission)'
case (N_F)
string = '(n,f)'

208
src/fission.f90 Normal file
View file

@ -0,0 +1,208 @@
module fission
use constants
use cross_section_header, only: Nuclide
use error, only: fatal_error
use search, only: binary_search
contains
!===============================================================================
! NU_TOTAL calculates the total number of neutrons emitted per fission for a
! given nuclide and incoming neutron energy
!===============================================================================
function nu_total(nuc, E) result(nu)
type(Nuclide), pointer :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: nu ! number of total neutrons emitted per fission
integer :: i ! loop index
integer :: j ! index on nu energy grid / precursor group
integer :: loc ! index before start of energies/nu values
integer :: NC ! number of polynomial coefficients
integer :: NR ! number of interpolation regions
integer :: NE ! number of energies tabulated
real(8) :: c ! polynomial coefficient
real(8) :: f ! interpolation factor
character(MAX_LINE_LEN) :: msg ! error message
if (nuc % nu_t_type == NU_NONE) then
msg = "No neutron emission data for table: " // nuc % name
call fatal_error(msg)
elseif (nuc % nu_t_type == NU_POLYNOMIAL) then
! determine number of coefficients
NC = int(nuc % nu_t_data(1))
! sum up polynomial in energy
nu = ZERO
do i = 0, NC - 1
c = nuc % nu_t_data(i+2)
nu = nu + c * E**i
end do
elseif (nuc % nu_t_type == NU_TABULAR) then
! determine number of interpolation regions -- as far as I can tell, no
! nu data has multiple interpolation regions. Furthermore, it seems all
! are lin-lin.
NR = int(nuc % nu_t_data(1))
if (NR /= 0) then
msg = "Multiple interpolation regions not supported while &
&attempting to determine total nu."
call fatal_error(msg)
end if
! determine number of energies
loc = 2 + 2*NR
NE = int(nuc % nu_t_data(loc))
! do binary search over tabuled energies to determine appropriate index
! and interpolation factor
if (E < nuc % nu_t_data(loc+1)) then
j = 1
f = ZERO
elseif (E > nuc % nu_t_data(loc+NE)) then
j = NE - 1
f = ONE
else
j = binary_search(nuc % nu_t_data(loc+1), NE, E)
f = (E - nuc % nu_t_data(loc+j)) / &
& (nuc % nu_t_data(loc+j+1) - nuc % nu_t_data(loc+j))
end if
! determine nu total
loc = loc + NE
nu = nuc % nu_t_data(loc+j) + f * &
& (nuc % nu_t_data(loc+j+1) - nuc % nu_t_data(loc+j))
end if
end function nu_total
!===============================================================================
! NU_PROMPT calculates the total number of prompt neutrons emitted per fission
! for a given nuclide and incoming neutron energy
!===============================================================================
function nu_prompt(nuc, E) result(nu)
type(Nuclide), pointer :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: nu ! number of prompt neutrons emitted per fission
integer :: i ! loop index
integer :: j ! index on nu energy grid / precursor group
integer :: loc ! index before start of energies/nu values
integer :: NC ! number of polynomial coefficients
integer :: NR ! number of interpolation regions
integer :: NE ! number of energies tabulated
real(8) :: c ! polynomial coefficient
real(8) :: f ! interpolation factor
character(MAX_LINE_LEN) :: msg ! error message
if (nuc % nu_p_type == NU_NONE) then
! since no prompt or delayed data is present, this means all neutron
! emission is prompt -- WARNING: This currently returns zero. The calling
! routine needs to know this situation is occurring since we don't want
! to call nu_total unnecessarily if it's already been called
nu = ZERO
elseif (nuc % nu_p_type == NU_POLYNOMIAL) then
! determine number of coefficients
NC = int(nuc % nu_p_data(1))
! sum up polynomial in energy
nu = ZERO
do i = 0, NC - 1
c = nuc % nu_p_data(i+2)
nu = nu + c * E**i
end do
elseif (nuc % nu_p_type == NU_TABULAR) then
! determine number of interpolation regions
NR = int(nuc % nu_p_data(1))
if (NR /= 0) then
msg = "Multiple interpolation regions not supported while &
&attempting to determine prompt nu."
call fatal_error(msg)
end if
! determine number of energies
loc = 2 + 2*NR
NE = int(nuc % nu_p_data(loc))
! do binary search over tabuled energies to determine appropriate index
! and interpolation factor
if (E < nuc % nu_p_data(loc+1)) then
j = 1
f = ZERO
elseif (E > nuc % nu_p_data(loc+NE)) then
j = NE - 1
f = ONE
else
j = binary_search(nuc % nu_p_data(loc+1), NE, E)
f = (E - nuc % nu_p_data(loc+j)) / &
& (nuc % nu_p_data(loc+j+1) - nuc % nu_p_data(loc+j))
end if
! determine nu total
loc = loc + NE
nu = nuc % nu_p_data(loc+j) + f * &
& (nuc % nu_p_data(loc+j+1) - nuc % nu_p_data(loc+j))
end if
end function nu_prompt
!===============================================================================
! NU_DELAYED calculates the total number of delayed neutrons emitted per fission
! for a given nuclide and incoming neutron energy
!===============================================================================
function nu_delayed(nuc, E) result(nu)
type(Nuclide), pointer :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: nu ! number of delayed neutrons emitted per fission
integer :: j ! index on nu energy grid / precursor group
integer :: loc ! index before start of energies/nu values
integer :: NR ! number of interpolation regions
integer :: NE ! number of energies tabulated
real(8) :: f ! interpolation factor
character(MAX_LINE_LEN) :: msg ! error message
if (nuc % nu_d_type == NU_NONE) then
nu = ZERO
elseif (nuc % nu_d_type == NU_TABULAR) then
! determine number of interpolation regions
NR = int(nuc % nu_d_data(1))
if (NR /= 0) then
msg = "Multiple interpolation regions not supported while &
&attempting to determine delayed nu."
call fatal_error(msg)
end if
! determine number of energies
loc = 2 + 2*NR
NE = int(nuc % nu_d_data(loc))
! do binary search over tabuled energies to determine appropriate index
! and interpolation factor
if (E < nuc % nu_d_data(loc+1)) then
j = 1
f = ZERO
elseif (E > nuc % nu_d_data(loc+NE)) then
j = NE - 1
f = ONE
else
j = binary_search(nuc % nu_d_data(loc+1), NE, E)
f = (E - nuc % nu_d_data(loc+j)) / &
& (nuc % nu_d_data(loc+j+1) - nuc % nu_d_data(loc+j))
end if
! determine nu total
loc = loc + NE
nu = nuc % nu_d_data(loc+j) + f * &
& (nuc % nu_d_data(loc+j+1) - nuc % nu_d_data(loc+j))
end if
end function nu_delayed
end module fission

14
src/geometry.f90
View file

@ -121,8 +121,9 @@ contains
cUniverse => univ
! set particle attributes
p%cell = univ % cells(i)
p%universe = dict_get_key(universe_dict, univ % uid)
p % cell = univ % cells(i)
p % universe = dict_get_key(universe_dict, univ % uid)
p % material = c % material
exit
elseif (c % type == CELL_FILL) then
lower_univ => universes(c % fill)
@ -337,7 +338,8 @@ contains
end if
else
! set current pointers
p%cell = index_cell
p % cell = index_cell
p % material = c % material
cCell => c
cMaterial => materials(cCell%material)
end if
@ -360,7 +362,8 @@ contains
end if
else
! set current pointers
p%cell = index_cell
p % cell = index_cell
p % material = c % material
cCell => c
cMaterial => materials(cCell%material)
end if
@ -375,7 +378,8 @@ contains
do i = 1, size(cells)
c => cells(i)
if (cell_contains(c, p)) then
p%cell = i
p % cell = i
p % material = c % material
cCell => c
cMaterial => materials(cCell%material)
return

7
src/global.f90
View file

@ -2,7 +2,8 @@ module global
use bank_header, only: Bank
use constants
use cross_section_header, only: Nuclide, SAB_Table, xsData
use cross_section_header, only: Nuclide, SAB_Table, xsData, NuclideMicroXS, &
MaterialMacroXS
use datatypes_header, only: DictionaryII, DictionaryCI
use geometry_header, only: Cell, Universe, Lattice, Surface
use material_header, only: Material
@ -51,6 +52,10 @@ module global
integer :: n_nuclides_total
integer :: n_sab_tables
! Cross section caches
type(NuclideMicroXS), allocatable :: micro_xs(:)
type(MaterialMacroXS) :: material_xs
! Current cell, surface, material
type(Cell), pointer :: cCell
type(Universe), pointer :: cUniverse

36
src/particle_header.f90
View file

@ -8,23 +8,25 @@ module particle_header
!===============================================================================
type Particle
integer(8) :: uid ! Unique ID
integer :: type ! Particle type (n, p, e, etc)
real(8) :: xyz(3) ! location
real(8) :: xyz_local(3) ! local location (after transformations)
real(8) :: uvw(3) ! directional cosines
real(8) :: wgt ! particle weight
real(8) :: E ! energy
integer :: IE ! index on energy grid
real(8) :: interp ! interpolation factor for energy grid
integer :: cell ! current cell
integer :: universe ! current universe
integer :: lattice ! current lattice
integer :: surface ! current surface
integer :: index_x ! lattice index for x direction
integer :: index_y ! lattice index for y direction
logical :: alive ! is particle alive?
integer :: n_coll ! # of collisions
integer(8) :: uid ! Unique ID
integer :: type ! Particle type (n, p, e, etc)
real(8) :: xyz(3) ! location
real(8) :: xyz_local(3) ! local location (after transformations)
real(8) :: uvw(3) ! directional cosines
real(8) :: wgt ! particle weight
real(8) :: E ! energy
integer :: IE ! index on energy grid
real(8) :: interp ! interpolation factor for energy grid
integer :: cell ! current cell
integer :: universe ! current universe
integer :: lattice ! current lattice
integer :: surface ! current surface
integer :: material ! current material
integer :: last_material ! last material
integer :: index_x ! lattice index for x direction
integer :: index_y ! lattice index for y direction
logical :: alive ! is particle alive?
integer :: n_coll ! # of collisions
end type Particle
end module particle_header

373
src/physics.f90
View file

@ -4,6 +4,7 @@ module physics
use cross_section_header, only: Nuclide, Reaction, DistEnergy
use endf, only: reaction_name
use error, only: fatal_error, warning
use fission, only: nu_total, nu_prompt, nu_delayed
use geometry, only: find_cell, dist_to_boundary, cross_surface, &
cross_lattice
use geometry_header, only: Universe, BASE_UNIVERSE
@ -29,12 +30,9 @@ contains
integer :: surf ! surface which particle is on
integer :: last_cell ! most recent cell particle was in
integer :: IE ! index on energy grid
real(8) :: d_to_boundary ! distance to nearest boundary
real(8) :: d_to_collision ! sampled distance to collision
real(8) :: distance ! distance particle travels
real(8) :: Sigma ! total cross-section
real(8) :: f ! interpolation factor
logical :: found_cell ! found cell which particle is in?
logical :: in_lattice ! is surface crossing in lattice?
character(MAX_LINE_LEN) :: msg ! output/error message
@ -46,8 +44,8 @@ contains
! if particle couldn't be located, print error
if (.not. found_cell) then
write(msg, 100) "Could not locate cell for particle at: ", p % xyz
100 format (A,3ES11.3)
write(msg, '(A,3ES11.3)') &
"Could not locate cell for particle at: ", p % xyz
call fatal_error(msg)
end if
end if
@ -63,22 +61,16 @@ contains
end if
! find energy index, interpolation factor
call find_energy_index(p)
do while (p % alive)
! Copy energy index and interpolation factor
IE = p % IE
f = p % interp
! Determine material total cross-section
Sigma = (1-f)*cMaterial%total_xs(IE) + f*cMaterial%total_xs(IE+1)
! Calculate microscopic and macroscopic cross sections
call calculate_xs(p)
! Find the distance to the nearest boundary
call dist_to_boundary(p, d_to_boundary, surf, in_lattice)
! Sample a distance to collision
d_to_collision = -log(rang()) / Sigma
d_to_collision = -log(rang()) / material_xs % total
! Select smaller of the two distances
distance = min(d_to_boundary, d_to_collision)
@ -107,6 +99,131 @@ 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
real(8) :: atom_density ! atom density of a nuclide
type(Material), pointer :: mat ! current material
! If the material is the same as the last material and the energy of the
! particle hasn't changed, we don't need to lookup cross sections again.
if (p % material == p % last_material) return
! Set all material macroscopic cross sections to zero
material_xs % total = ZERO
material_xs % scatter = 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)
! 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)
call calculate_nuclide_xs(p, index_nuclide)
! 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)
type(Particle), pointer :: p
integer, intent(in) :: index_nuclide ! index into nuclides array
integer :: i ! index into nuclides array
integer :: IE ! index on nuclide energy grid
real(8) :: f ! interpolation factor on nuclide energy grid
real(8) :: nu_t ! total number of neutrons emitted per fission
type(Nuclide), pointer :: nuc ! pointer to nuclide cross section table
! Copy index of nuclide
i = index_nuclide
! Set pointer to nuclide
nuc => nuclides(i)
! TODO: Check if last energy/temp combination is same as current. If so, we
! can return.
! 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
! search nuclide energy grid
IE = nuc % grid_index(p % IE)
f = (p%E - nuc%energy(IE))/(nuc%energy(IE+1) - nuc%energy(IE))
micro_xs(i) % index_grid = IE
micro_xs(i) % interp_factor = f
! Initialize nuclide cross-sections to zero
micro_xs(i) % nu_fission = ZERO
! Calculate microscopic nuclide total cross section
micro_xs(i) % total = &
(ONE-f) * nuc % total(IE) + f * nuc % total(IE+1)
! Calculate microscopic nuclide total cross section
micro_xs(i) % elastic = &
(ONE-f) * nuc % elastic(IE) + f * nuc % elastic(IE+1)
! Calculate microscopic nuclide absorption cross section
micro_xs(i) % absorption = &
(ONE-f) * nuc % absorption(IE) + f * nuc % absorption(IE+1)
if (nuc % fissionable) then
! Calculate microscopic nuclide total cross section
micro_xs(i) % fission = &
(ONE-f) * nuc % fission(IE) + f * nuc % fission(IE+1)
! Calculate microscopic nuclide nu-fission cross section
nu_t = nu_total(nuc, p % E)
micro_xs(i) % nu_fission = nu_t * micro_xs(i) % fission
end if
end subroutine calculate_nuclide_xs
!===============================================================================
! FIND_ENERGY_INDEX determines the index on the union energy grid and the
! interpolation factor for a particle at a certain energy
@ -153,63 +270,47 @@ contains
type(Particle), pointer :: p
integer :: i ! index over nuclides in a material
integer :: n_nuclides ! number of nuclides in material
integer :: index_nuclide
integer :: IE ! index on nuclide energy grid
real(8) :: f ! interpolation factor
real(8) :: sigma ! microscopic total xs for nuclide
real(8) :: total ! total macroscopic xs for material
real(8) :: density_i ! atom density of nuclide
real(8) :: prob ! cumulative probability
real(8) :: r1 ! random number
real(8) :: cutoff ! random number
real(8) :: flux ! collision estimator of flux
real(8), allocatable :: Sigma_rxn(:) ! macroscopic xs for each nuclide
character(MAX_LINE_LEN) :: msg ! output/error message
real(8) :: atom_density ! atom density of nuclide in atom/b-cm
character(MAX_LINE_LEN) :: msg ! output/error message
type(Material), pointer :: mat
type(Nuclide), pointer :: nuc
type(Reaction), pointer :: rxn
mat => cMaterial
! calculate total cross-section for each nuclide at current energy in order
! to create discrete pdf for sampling nuclide
n_nuclides = mat % n_nuclides
allocate(Sigma_rxn(n_nuclides))
do i = 1, n_nuclides
nuc => nuclides(mat % nuclide(i))
! determine nuclide atom density
density_i = mat % atom_density(i)
! search nuclide energy grid
IE = nuc%grid_index(p % IE)
f = (p%E - nuc%energy(IE))/(nuc%energy(IE+1) - nuc%energy(IE))
! calculate nuclide macroscopic cross-section
Sigma_rxn(i) = density_i*((ONE-f)*nuc%Sigma_t(IE) + &
& f*(nuc%Sigma_t(IE+1)))
end do
total = sum(Sigma_rxn)
mat => materials(p % material)
! sample nuclide
r1 = rang()
i = 0
total = material_xs % total
cutoff = rang() * total
prob = ZERO
do i = 1, n_nuclides
prob = prob + Sigma_rxn(i) / total
if (r1 < prob) exit
do while (prob < cutoff)
i = i + 1
if (i > mat % n_nuclides) then
msg = "Did not sample any nuclide during collision."
call fatal_error(msg)
end if
index_nuclide = mat % nuclide(i)
atom_density = mat % atom_density(i)
sigma = atom_density * micro_xs(index_nuclide) % total
prob = prob + sigma
end do
! Get table, total xs, interpolation factor
density_i = mat % atom_density(i)
nuc => nuclides(mat%nuclide(i))
sigma = Sigma_rxn(i) / density_i
nuc => nuclides(index_nuclide)
IE = nuc%grid_index(p % IE)
f = (p%E - nuc%energy(IE))/(nuc%energy(IE+1) - nuc%energy(IE))
! free memory
deallocate(Sigma_rxn)
! sample reaction channel
r1 = rang()*sigma
cutoff = rang() * sigma / atom_density
prob = ZERO
do i = 1, nuc % n_reaction
rxn => nuc % reactions(i)
@ -225,7 +326,7 @@ contains
! add to cumulative probability
prob = prob + ((ONE-f)*rxn%sigma(IE-rxn%IE+1) &
& + f*(rxn%sigma(IE-rxn%IE+2)))
if (r1 < prob) exit
if (cutoff < prob) exit
end do
if (verbosity >= 10) then
@ -243,8 +344,8 @@ contains
case (N_NA, N_N3A, N_NP, N_N2A, N_ND, N_NT, N_N3HE, &
& N_NT2A, N_N2P, N_NPA, N_N1 : N_NC)
call inelastic_scatter(p, nuc, rxn)
case (FISSION, N_F, N_NF, N_2NF, N_3NF)
call n_fission(p, nuc, rxn)
case (N_FISSION, N_F, N_NF, N_2NF, N_3NF)
call n_fission_event(p, nuc, rxn)
case (N_GAMMA : N_DA)
call n_absorption(p)
case default
@ -350,7 +451,7 @@ contains
! with implicit absorption, namely sampling of the number of neutrons!
!===============================================================================
subroutine n_fission(p, nuc, rxn)
subroutine n_fission_event(p, nuc, rxn)
type(Particle), pointer :: p
type(Nuclide), pointer :: nuc
@ -360,18 +461,16 @@ contains
integer :: j ! index on nu energy grid / precursor group
integer :: k ! index on precursor yield grid
integer :: loc ! index before start of energies/nu values
integer :: NC ! number of polynomial coefficients
integer :: NR ! number of interpolation regions
integer :: NE ! number of energies tabulated
integer :: nu ! actual number of neutrons produced
integer :: law ! energy distribution law
real(8) :: E ! incoming energy of neutron
real(8) :: E_out ! outgoing energy of fission neutron
real(8) :: c ! polynomial coefficient
real(8) :: f ! interpolation factor
real(8) :: nu_total ! total nu
real(8) :: nu_prompt ! prompt nu
real(8) :: nu_delay ! delayed nu
real(8) :: nu_t ! total nu
real(8) :: nu_p ! prompt nu
real(8) :: nu_d ! delayed nu
real(8) :: mu ! fission neutron angular cosine
real(8) :: phi ! fission neutron azimuthal angle
real(8) :: beta ! delayed neutron fraction
@ -383,152 +482,30 @@ contains
! copy energy of neutron
E = p % E
! ==========================================================================
! DETERMINE TOTAL NU
if (nuc % nu_t_type == NU_NONE) then
msg = "No neutron emission data for table: " // nuc % name
call fatal_error(msg)
elseif (nuc % nu_t_type == NU_POLYNOMIAL) then
! determine number of coefficients
NC = int(nuc % nu_t_data(1))
! Determine total nu
nu_t = nu_total(nuc, E)
! sum up polynomial in energy
nu_total = ZERO
do i = 0, NC - 1
c = nuc % nu_t_data(i+2)
nu_total = nu_total + c * E**i
end do
elseif (nuc % nu_t_type == NU_TABULAR) then
! determine number of interpolation regions -- as far as I can tell, no
! nu data has multiple interpolation regions. Furthermore, it seems all
! are lin-lin.
NR = int(nuc % nu_t_data(1))
if (NR /= 0) then
msg = "Multiple interpolation regions not supported while &
&attempting to determine total nu."
call fatal_error(msg)
end if
! determine number of energies
loc = 2 + 2*NR
NE = int(nuc % nu_t_data(loc))
! do binary search over tabuled energies to determine appropriate index
! and interpolation factor
if (E < nuc % nu_t_data(loc+1)) then
j = 1
f = ZERO
elseif (E > nuc % nu_t_data(loc+NE)) then
j = NE - 1
f = ONE
else
j = binary_search(nuc % nu_t_data(loc+1), NE, E)
f = (E - nuc % nu_t_data(loc+j)) / &
& (nuc % nu_t_data(loc+j+1) - nuc % nu_t_data(loc+j))
end if
! determine nu total
loc = loc + NE
nu_total = nuc % nu_t_data(loc+j) + f * &
& (nuc % nu_t_data(loc+j+1) - nuc % nu_t_data(loc+j))
! Determine prompt nu
if (nuc % nu_p_type == NU_NONE) then
nu_p = nu_t
else
nu_p = nu_prompt(nuc, E)
end if
! ==========================================================================
! DETERMINE PROMPT NU
if (nuc % nu_p_type == NU_NONE) then
! since no prompt or delayed data is present, this means all neutron
! emission is prompt
nu_prompt = nu_total
elseif (nuc % nu_p_type == NU_POLYNOMIAL) then
! determine number of coefficients
NC = int(nuc % nu_p_data(1))
! Determine delayed nu
nu_d = nu_delayed(nuc, E)
! sum up polynomial in energy
nu_prompt = ZERO
do i = 0, NC - 1
c = nuc % nu_p_data(i+2)
nu_prompt = nu_prompt + c * E**i
end do
elseif (nuc % nu_p_type == NU_TABULAR) then
! determine number of interpolation regions
NR = int(nuc % nu_p_data(1))
if (NR /= 0) then
msg = "Multiple interpolation regions not supported while &
&attempting to determine prompt nu."
call fatal_error(msg)
end if
! determine number of energies
loc = 2 + 2*NR
NE = int(nuc % nu_p_data(loc))
! do binary search over tabuled energies to determine appropriate index
! and interpolation factor
if (E < nuc % nu_p_data(loc+1)) then
j = 1
f = ZERO
elseif (E > nuc % nu_p_data(loc+NE)) then
j = NE - 1
f = ONE
else
j = binary_search(nuc % nu_p_data(loc+1), NE, E)
f = (E - nuc % nu_p_data(loc+j)) / &
& (nuc % nu_p_data(loc+j+1) - nuc % nu_p_data(loc+j))
end if
! determine nu total
loc = loc + NE
nu_prompt = nuc % nu_p_data(loc+j) + f * &
& (nuc % nu_p_data(loc+j+1) - nuc % nu_p_data(loc+j))
end if
! ==========================================================================
! DETERMINE DELAYED NU
if (nuc % nu_d_type == NU_NONE) then
nu_delay = ZERO
elseif (nuc % nu_d_type == NU_TABULAR) then
! determine number of interpolation regions
NR = int(nuc % nu_d_data(1))
if (NR /= 0) then
msg = "Multiple interpolation regions not supported while &
&attempting to determine delayed nu."
call fatal_error(msg)
end if
! determine number of energies
loc = 2 + 2*NR
NE = int(nuc % nu_d_data(loc))
! do binary search over tabuled energies to determine appropriate index
! and interpolation factor
if (E < nuc % nu_d_data(loc+1)) then
j = 1
f = ZERO
elseif (E > nuc % nu_d_data(loc+NE)) then
j = NE - 1
f = ONE
else
j = binary_search(nuc % nu_d_data(loc+1), NE, E)
f = (E - nuc % nu_d_data(loc+j)) / &
& (nuc % nu_d_data(loc+j+1) - nuc % nu_d_data(loc+j))
end if
! determine nu total
loc = loc + NE
nu_delay = nuc % nu_d_data(loc+j) + f * &
& (nuc % nu_d_data(loc+j+1) - nuc % nu_d_data(loc+j))
end if
beta = nu_delay / nu_total
! Determine delayed neutron fraction
beta = nu_d / nu_t
! TODO: Heat generation from fission
! Sample number of neutrons produced
nu_total = p%wgt / keff * nu_total
if (rang() > nu_total - int(nu_total)) then
nu = int(nu_total)
nu_t = p%wgt / keff * nu_t
if (rang() > nu_t - int(nu_t)) then
nu = int(nu_t)
else
nu = int(nu_total) + 1
nu = int(nu_t) + 1
end if
! Bank source neutrons
@ -633,7 +610,7 @@ contains
! kill original neutron
p % alive = .false.
end subroutine n_fission
end subroutine n_fission_event
!===============================================================================
! INELASTIC_SCATTER handles all reactions with a single secondary neutron (other