crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Initial support for multiple S(a,b) tables. Still need sorting of i_sab_nuclides and i_sab_tables.

Browse source
This commit is contained in:
Paul Romano 2013-02-06 19:19:20 -05:00
parent b85ef9ae76
commit c4041512b0
8 changed files with 118 additions and 78 deletions
Download patch file Download diff file Expand all files Collapse all files

36
src/ace.F90
View file

@ -34,6 +34,7 @@ contains
integer :: i ! index in materials array
integer :: j ! index over nuclides in material
integer :: k ! index over S(a,b) tables in material
integer :: i_listing ! index in xs_listings array
integer :: i_nuclide ! index in nuclides
integer :: i_sab ! index in sab_tables
@ -79,8 +80,9 @@ contains
end if
end do
if (mat % has_sab_table) then
name = mat % sab_name
do k = 1, mat % n_sab
! Get name of S(a,b) table
name = mat % sab_names(k)
if (.not. already_read % contains(name)) then
i_listing = xs_listing_dict % get_key(name)
@ -93,7 +95,7 @@ contains
! Add name to dictionary
call already_read % add(name)
end if
end if
end do
end do
! ==========================================================================
@ -102,26 +104,26 @@ contains
do i = 1, n_materials
mat => materials(i)
if (mat % has_sab_table) then
! In order to know which nuclide the S(a,b) table applies to, we need
! to search through the list of nuclides for one which has a matching
! zaid
sab => sab_tables(mat % sab_table)
do k = 1, mat % n_sab
! In order to know which nuclide the S(a,b) table applies to, we need to
! search through the list of nuclides for one which has a matching zaid
sab => sab_tables(mat % i_sab_tables(k))
do j = 1, mat % n_nuclides
nuc => nuclides(mat % nuclide(j))
if (nuc % zaid == sab % zaid) then
mat % sab_nuclide = j
! Loop through nuclides and find match
FIND_NUCLIDE: do j = 1, mat % n_nuclides
if (nuclides(mat % nuclide(j)) % zaid == sab % zaid) then
mat % i_sab_nuclides(k) = j
exit FIND_NUCLIDE
end if
end do
end do FIND_NUCLIDE
! Check to make sure S(a,b) table matched a nuclide
if (mat % sab_nuclide == 0) then
message = "S(a,b) table " // trim(mat % sab_name) // " did not match &
&any nuclide on material " // trim(to_str(mat % id))
if (mat % i_sab_nuclides(k) == NONE) then
message = "S(a,b) table " // trim(mat % sab_names(k)) // " did not &
&match any nuclide on material " // trim(to_str(mat % id))
call fatal_error()
end if
end if
end do
end do
end subroutine read_xs

2
src/ace_header.F90
View file

@ -195,7 +195,7 @@ module ace_header
real(8) :: kappa_fission ! microscopic energy-released from fission
! Information for S(a,b) use
logical :: use_sab ! in S(a,b) energy range?
integer :: index_sab ! index in sab_tables (zero means no table)
real(8) :: elastic_sab ! microscopic elastic scattering on S(a,b) table
! Information for URR probability table use

54
src/cross_section.F90
View file

@ -23,8 +23,9 @@ contains
integer :: i ! loop index over nuclides
integer :: i_nuclide ! index into nuclides array
integer :: i_sab ! index into sab_tables array
integer :: j ! index in mat % i_sab_nuclides
real(8) :: atom_density ! atom density of a nuclide
real(8) :: sab_threshold ! threshold for S(a,b) table
logical :: check_sab ! should we check for S(a,b) table?
type(Material), pointer :: mat => null() ! current material
! Set all material macroscopic cross sections to zero
@ -43,31 +44,50 @@ contains
! Find energy index on unionized grid
if (grid_method == GRID_UNION) call find_energy_index()
! 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
! Determine if this material has S(a,b) tables
check_sab = (mat % n_sab > 0)
! Initialize position in i_sab_nuclides
j = 1
! Add contribution from each nuclide in material
do i = 1, mat % n_nuclides
! ========================================================================
! CHECK FOR S(A,B) TABLE
i_sab = 0
! Check if this nuclide matches one of the S(a,b) tables specified -- this
! relies on i_sab_nuclides being in sorted order
if (check_sab .and. i == mat % i_sab_nuclides(j)) then
! Get index in sab_tables
i_sab = mat % i_sab_tables(j)
! If particle energy is greater than the highest energy for the S(a,b)
! table, don't use the S(a,b) table
if (p % E > sab_tables(i_sab) % threshold_inelastic) i_sab = 0
! Increment position in i_sab_nuclides
j = j + 1
! Don't check for S(a,b) tables if there are no more left
if (j > mat % n_sab) check_sab = .false.
end if
! ========================================================================
! CALCULATE MICROSCOPIC CROSS SECTION
! Determine microscopic cross sections for this nuclide
i_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
i_sab = mat % sab_table
else
i_sab = 0
end if
! Calculate microscopic cross section for this nuclide
if (p % E /= micro_xs(i_nuclide) % last_E) then
call calculate_nuclide_xs(i_nuclide, i_sab)
end if
! ========================================================================
! ADD TO MACROSCOPIC CROSS SECTION
! Copy atom density of nuclide in material
atom_density = mat % atom_density(i)
@ -148,7 +168,7 @@ contains
micro_xs(i_nuclide) % interp_factor = f
! Initialize sab treatment to false
micro_xs(i_nuclide) % use_sab = .false.
micro_xs(i_nuclide) % index_sab = NONE
micro_xs(i_nuclide) % elastic_sab = ZERO
micro_xs(i_nuclide) % use_ptable = .false.
@ -231,7 +251,7 @@ contains
type(SAlphaBeta), pointer :: sab => null()
! Set flag that S(a,b) treatment should be used for scattering
micro_xs(i_nuclide) % use_sab = .true.
micro_xs(i_nuclide) % index_sab = i_sab
! Get pointer to S(a,b) table
sab => sab_tables(i_sab)

8
src/hdf5_interface.F90
View file

@ -464,9 +464,11 @@ contains
dims, m % atom_density, hdf5_err)
! Write S(a,b) information if present
if (m % has_sab_table) then
call h5ltmake_dataset_string_f(temp_group, "sab_table", &
m % sab_name, hdf5_err)
if (m % n_sab > 0) then
call hdf5_write_integer(temp_group, "i_sab_nuclides", &
m % i_sab_nuclides)
call hdf5_write_integer(temp_group, "i_sab_tables", &
m % i_sab_tables)
end if
! Close group for i-th material

72
src/input_xml.F90
View file

@ -952,6 +952,7 @@ contains
integer :: i ! loop index for materials
integer :: j ! loop index for nuclides
integer :: n ! number of nuclides
integer :: n_sab ! number of sab tables for a material
integer :: index_list ! index in xs_listings array
integer :: index_nuclide ! index in nuclides
integer :: index_sab ! index in sab_tables
@ -1152,43 +1153,52 @@ contains
! =======================================================================
! READ AND PARSE <sab> TAG FOR S(a,b) DATA
if (size(material_(i) % sab) == 1) then
! Get pointer to S(a,b) table
sab => material_(i) % sab(1)
n_sab = size(material_(i) % sab)
if (n_sab > 0) then
! Set number of S(a,b) tables
mat % n_sab = n_sab
! Determine name of S(a,b) table
name = trim(sab % name) // "." // trim(sab % xs)
mat % sab_name = name
! Allocate names and indices for nuclides and tables
allocate(mat % sab_names(n_sab))
allocate(mat % i_sab_nuclides(n_sab))
allocate(mat % i_sab_tables(n_sab))
! Check that this nuclide is listed in the cross_sections.xml file
if (.not. xs_listing_dict % has_key(name)) then
message = "Could not find S(a,b) table " // trim(name) // &
" in cross_sections.xml file!"
call fatal_error()
end if
mat % has_sab_table = .true.
! Initialize i_sab_nuclides
mat % i_sab_nuclides = NONE
! Find index in xs_listing and set the name and alias according to the
! listing
index_list = xs_listing_dict % get_key(name)
name = xs_listings(index_list) % name
alias = xs_listings(index_list) % alias
do j = 1, n_sab
! Get pointer to S(a,b) table
sab => material_(i) % sab(j)
! If this S(a,b) table hasn't been encountered yet, we need to add its
! name and alias to the sab_dict
if (.not. sab_dict % has_key(name)) then
index_sab = index_sab + 1
mat % sab_table = index_sab
! Determine name of S(a,b) table
name = trim(sab % name) // "." // trim(sab % xs)
mat % sab_names(j) = name
call sab_dict % add_key(name, index_sab)
call sab_dict % add_key(alias, index_sab)
else
mat % sab_table = sab_dict % get_key(name)
end if
! Check that this nuclide is listed in the cross_sections.xml file
if (.not. xs_listing_dict % has_key(name)) then
message = "Could not find S(a,b) table " // trim(name) // &
" in cross_sections.xml file!"
call fatal_error()
end if
elseif (size(material_(i) % sab) > 1) then
message = "Cannot have multiple S(a,b) tables on a single material."
call fatal_error()
! Find index in xs_listing and set the name and alias according to the
! listing
index_list = xs_listing_dict % get_key(name)
name = xs_listings(index_list) % name
alias = xs_listings(index_list) % alias
! If this S(a,b) table hasn't been encountered yet, we need to add its
! name and alias to the sab_dict
if (.not. sab_dict % has_key(name)) then
index_sab = index_sab + 1
mat % i_sab_tables(j) = index_sab
call sab_dict % add_key(name, index_sab)
call sab_dict % add_key(alias, index_sab)
else
mat % i_sab_tables(j) = sab_dict % get_key(name)
end if
end do
end if
! Add material to dictionary

12
src/material_header.F90
View file

@ -9,16 +9,18 @@ module material_header
type Material
integer :: id ! unique identifier
integer :: n_nuclides ! number of nuclides
character(12), allocatable :: names(:) ! isotope names
integer, allocatable :: nuclide(:) ! index in nuclides array
real(8) :: density ! total atom density in atom/b-cm
real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm
! S(a,b) data references
logical :: has_sab_table = .false.
character(12) :: sab_name ! name of S(a,b) table
integer :: sab_table = 0 ! index in sab_tables
integer :: sab_nuclide = 0 ! index of nuclide which has S(a,b) table
integer :: n_sab = 0 ! number of S(a,b) tables
integer, allocatable :: i_sab_nuclides(:) ! index of corresponding nuclide
integer, allocatable :: i_sab_tables(:) ! index in sab_tables
! Temporary names read during initialization
character(12), allocatable :: names(:) ! isotope names
character(12), allocatable :: sab_names(:) ! name of S(a,b) table
end type Material
end module material_header

8
src/output.F90
View file

@ -624,8 +624,12 @@ contains
end do
! Write information on S(a,b) table
if (mat % has_sab_table) then
write(unit_,*) ' S(a,b) table = ' // trim(mat % sab_name)
if (mat % n_sab > 0) then
write(unit_,*) ' S(a,b) tables:'
do i = 1, mat % n_sab
write(unit_,*) ' ' // trim(&
sab_tables(mat % i_sab_tables(i)) % name)
end do
end if
write(unit_,*)

4
src/physics.F90
View file

@ -422,9 +422,9 @@ contains
! =======================================================================
! ELASTIC SCATTERING
if (micro_xs(i_nuclide) % use_sab) then
if (micro_xs(i_nuclide) % index_sab /= NONE) then
! S(a,b) scattering
call sab_scatter(i_nuclide, mat % sab_table)
call sab_scatter(i_nuclide, micro_xs(i_nuclide) % index_sab)
else
! get pointer to elastic scattering reaction