crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-29 06:35:48 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Switch temperature tag from materials to cells

Browse source 
Note that this allows for distributed temperatures
This commit is contained in:
Sterling Harper 2016-02-03 20:32:20 -05:00
parent c27cb008d3
commit 4748f43f9a
7 changed files with 168 additions and 58 deletions
Download patch file Download diff file Expand all files Collapse all files

36
src/cross_section.F90
View file

@ -105,10 +105,13 @@ contains
i_nuclide = mat % nuclide(i)
! Calculate microscopic cross section for this nuclide
if (p % E /= micro_xs(i_nuclide) % last_E .or. mat % sqrtkT /= micro_xs(i_nuclide) % last_sqrtkT) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, i_grid, mat % sqrtkT)
if (p % E /= micro_xs(i_nuclide) % last_E &
.or. p % sqrtkT /= micro_xs(i_nuclide) % last_sqrtkT) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, &
i_grid, p % sqrtkT)
else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, i_grid, mat % sqrtkT)
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, &
i_grid, p % sqrtkT)
end if
! ========================================================================
@ -145,7 +148,8 @@ contains
! given index in the nuclides array at the energy of the given particle
!===============================================================================
subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_mat, i_nuc_mat, i_log_union, sqrtkT)
subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_mat, i_nuc_mat, &
i_log_union, sqrtkT)
integer, intent(in) :: i_nuclide ! index into nuclides array
integer, intent(in) :: i_sab ! index into sab_tables array
real(8), intent(in) :: E ! energy
@ -590,14 +594,20 @@ contains
! sections in the resolved resonance regions
!===============================================================================
subroutine multipole_eval(multipole, Emev, sqrtkT, sigT, sigA, sigF)
type(MultipoleArray), intent(in) :: multipole ! The windowed multipole object to process.
real(8), intent(in) :: Emev ! The energy at which to evaluate the cross section in MeV
real(8), intent(in) :: sqrtkT ! The temperature in the form sqrt(kT (in eV)), at which to evaluate the cross section.
real(8), intent(out) :: sigT ! Total cross section
real(8), intent(out) :: sigA ! Absorption cross section
real(8), intent(out) :: sigF ! Fission cross section
complex(8) :: psi_ki ! The value of the psi-ki function for the asymptotic form
subroutine multipole_eval(multipole, Emev, sqrtkT_, sigT, sigA, sigF)
type(MultipoleArray), intent(in) :: multipole ! The windowed multipole
! object to process.
real(8), intent(in) :: Emev ! The energy at which to
! evaluate the cross section
! in MeV
real(8), intent(in) :: sqrtkT_ ! The temperature in the form
! sqrt(kT (in MeV)), at which
! to evaluate the XS.
real(8), intent(out) :: sigT ! Total cross section
real(8), intent(out) :: sigA ! Absorption cross section
real(8), intent(out) :: sigF ! Fission cross section
complex(8) :: psi_ki ! The value of the psi-ki function for the asymptotic
! form
complex(8) :: c_temp ! complex temporary variable
complex(8) :: w_val ! The faddeeva function evaluated at Z
complex(8) :: Z ! sqrt(atomic weight ratio / kT) * (sqrt(E) - pole)
@ -607,6 +617,7 @@ contains
real(8) :: dopp_ecoef ! sqrt(atomic weight ratio * pi / kT) / E
real(8) :: temp ! real temporary value
real(8) :: E ! energy, eV
real(8) :: sqrtkT ! sqrt(kT (in eV))
integer :: iP ! index of pole
integer :: iC ! index of curvefit
integer :: iW ! index of window
@ -617,6 +628,7 @@ contains
! Convert to eV
E = Emev * 1.0e6_8
sqrtkT = sqrtkT_ * 1.0e3_8
sqrtE = sqrt(E)
invE = ONE/E

30
src/geometry.F90
View file

@ -279,6 +279,36 @@ contains
p % material = c % material(offset + 1)
end if
! Set the particle temperature
if (size(c % sqrtkT) == 1) then
! Only one temperature for this cell; assign that one to the particle.
p % sqrtkT = c % sqrtkT(1)
else
! Distributed instances of this cell have different temperatures.
! Determine which instance this is and assign the matching temp.
distribcell_index = c % distribcell_index
offset = 0
do k = 1, p % n_coord
if (cells(p % coord(k) % cell) % type == CELL_FILL) then
offset = offset + cells(p % coord(k) % cell) % &
offset(distribcell_index)
elseif (cells(p % coord(k) % cell) % type == CELL_LATTICE) then
if (lattices(p % coord(k + 1) % lattice) % obj &
% are_valid_indices([&
p % coord(k + 1) % lattice_x, &
p % coord(k + 1) % lattice_y, &
p % coord(k + 1) % lattice_z])) then
offset = offset + lattices(p % coord(k + 1) % lattice) % obj % &
offset(distribcell_index, &
p % coord(k + 1) % lattice_x, &
p % coord(k + 1) % lattice_y, &
p % coord(k + 1) % lattice_z)
end if
end if
end do
p % sqrtkT = c % sqrtkT(offset + 1)
end if
elseif (c % type == CELL_FILL) then CELL_TYPE
! ======================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL

3
src/geometry_header.F90
View file

@ -139,6 +139,9 @@ module geometry_header
! only)
integer :: distribcell_index ! Index corresponding to this cell in
! distribcell arrays
real(8), allocatable :: sqrtkT(:) ! Square root of k_Boltzmann *
! temperature in MeV. Multiple for
! distribcell
! Rotation matrix and translation vector
real(8), allocatable :: translation(:)

78
src/initialize.F90
View file

@ -120,6 +120,9 @@ contains
! Create linked lists for multiple instances of the same nuclide
call same_nuclide_list()
! Set undefined cell temperatures to match the material data.
call lookup_material_temperatures()
! Construct unionized or log energy grid for cross-sections
select case (grid_method)
case (GRID_NUCLIDE)
@ -967,10 +970,11 @@ contains
end do
end do
! We also need distribcell if any distributed materials are present.
! We also need distribcell if any distributed materials or distributed
! temperatues are present.
if (.not. distribcell_active) then
do i = 1, n_cells
if (size(cells(i) % material) > 1) then
if (size(cells(i) % material) > 1 .or. size(cells(i) % sqrtkT) > 1) then
distribcell_active = .true.
exit
end if
@ -995,8 +999,8 @@ contains
end do
end do
! Make sure the number of materials matches the number of cell instances for
! distributed materials.
! Make sure the number of materials and temperatures matches the number of
! cell instances.
do i = 1, n_cells
associate (c => cells(i))
if (size(c % material) > 1) then
@ -1008,6 +1012,15 @@ contains
&equal one or the number of instances.")
end if
end if
if (size(c % sqrtkT) > 1) then
if (size(c % sqrtkT) /= c % instances) then
call fatal_error("Cell " // trim(to_str(c % id)) // " was &
&specified with " // trim(to_str(size(c % sqrtkT))) &
// " temperatures but has " // trim(to_str(c % instances)) &
// " distributed instances. The number of temperatures must &
&equal one or the number of instances.")
end if
end if
end associate
end do
@ -1049,9 +1062,9 @@ contains
end do
end do
! List all cells with multiple (distributed) materials.
! List all cells with multiple (distributed) materials or temperatures.
do i = 1, n_cells
if (size(cells(i) % material) > 1) then
if (size(cells(i) % material) > 1 .or. size(cells(i) % sqrtkT) > 1) then
call cell_list % add(i)
end if
end do
@ -1120,4 +1133,57 @@ contains
end subroutine allocate_offsets
!===============================================================================
! LOOKUP_MATERIAL_TEMPERATURES If any cells have undefined temperatures, try to
! find their temperatures from material data.
!===============================================================================
subroutine lookup_material_temperatures()
integer :: i, j, k
real(8) :: min_temp
logical :: warning_given
warning_given = .false.
do i = 1, n_cells
! Ignore non-normal cells and cells with defined temperature.
if (cells(i) % type /= CELL_NORMAL) cycle
if (cells(i) % sqrtkT(1) /= ERROR_REAL) cycle
! Set the number of temperatures equal to the number of materials.
deallocate(cells(i) % sqrtkT)
allocate(cells(i) % sqrtkT(size(cells(i) % material)))
! Check each of the cell materials for temperature data.
do j = 1, size(cells(i) % material)
! Arbitrarily set void regions to 0K.
if (cells(i) % material(j) == MATERIAL_VOID) then
cells(i) % sqrtkT(j) = ZERO
cycle
end if
associate (mat => materials(cells(i) % material(j)))
! Find the temperature of the coldest nuclide.
min_temp = nuclides(mat % nuclide(1)) % kT
do k = 2, mat % n_nuclides
! Warn the user if the nuclides don't have identical temperatues.
if (nuclides(mat % nuclide(k)) % kT /= min_temp &
.and. .not. warning_given) then
call warning("OpenMC cannot &
&identify the temperature of at least one cell. For the &
&purposes of multipole cross section evaluations, all cells &
&with unknown temperature will be set to the coldest &
&temperature found in the nuclear data for that cell's &
&material")
warning_given = .true.
end if
min_temp = min(min_temp, nuclides(mat % nuclide(k)) % kT)
end do
! Set the temperature for this cell instance.
cells(i) % sqrtkT(j) = sqrt(min_temp)
end associate
end do
end do
end subroutine lookup_material_temperatures
end module initialize

72
src/input_xml.F90
View file

@ -1297,6 +1297,40 @@ contains
call get_node_array(node_cell, "translation", c % translation)
end if
! Read cell temperatures. If the temperature is not specified, set it to
! ERROR_REAL for now. During initialization we'll replace ERROR_REAL with
! the temperature from the material data.
if (check_for_node(node_cell, "temperature")) then
n = get_arraysize_double(node_cell, "temperature")
if (n > 0) then
! Make sure this is a "normal" cell.
if (c % material(1) == NONE) call fatal_error("Cell " &
// trim(to_str(c % id)) // " was specified with a temperature &
&but no material. Temperature specification is only valid for &
&cells filled with a material.")
! Copy in temperatures
allocate(c % sqrtkT(n))
call get_node_array(node_cell, "temperature", c % sqrtkT)
! Make sure all temperatues are positive
do j = 1, size(c % sqrtkT)
if (c % sqrtkT(j) < ZERO) call fatal_error("Cell " &
// trim(to_str(c % id)) // " was specified with a negative &
&temperature. All cell temperatures must be non-negative.")
end do
! Convert to sqrt(kT)
c % sqrtkT(:) = sqrt(K_BOLTZMANN * c % sqrtkT(:))
else
allocate(c % sqrtkT(1))
c % sqrtkT(1) = ERROR_REAL
end if
else
allocate(c % sqrtkT(1))
c % sqrtkT = ERROR_REAL
end if
! Add cell to dictionary
call cell_dict % add_key(c % id, i)
@ -1854,7 +1888,6 @@ contains
real(8) :: temp_dble ! temporary double prec. real
logical :: file_exists ! does materials.xml exist?
logical :: sum_density ! density is taken to be sum of nuclide densities
logical :: temp_known ! Has the temperature yet been defined?
character(12) :: name ! name of isotope, e.g. 92235.03c
character(12) :: alias ! alias of nuclide, e.g. U-235.03c
character(MAX_WORD_LEN) :: units ! units on density
@ -1938,17 +1971,6 @@ contains
cycle
end if
! =======================================================================
! READ AND PARSE <temperature> TAG
if (check_for_node(node_mat, "temperature")) then
call get_node_ptr(node_mat, "temperature", node_temp)
call get_node_value(node_temp, "value", temp_dble)
mat % sqrtkT = sqrt(temp_dble * K_BOLTZMANN * 1.0D6)
temp_known = .true.
else
temp_known = .false.
end if
! =======================================================================
! READ AND PARSE <density> TAG
@ -2055,16 +2077,6 @@ contains
call get_node_value(node_nuc, "xs", name)
name = trim(temp_str) // "." // trim(name)
! If needed, look up temperature
if (.not. temp_known) then
! Find xs_listing and set the name/alias according to the listing
index_list = xs_listing_dict % get_key(to_lower(name))
if(xs_listings(index_list) % kT /= 0.0_8) then
mat % sqrtkT = sqrt(xs_listings(index_list) % kT * 1.0D6)
temp_known = .true.
end if
end if
! save name and density to list
call list_names % append(name)
@ -2153,16 +2165,6 @@ contains
temp_str = "data"
end if
! If still needed, look up temperature
if (.not. temp_known) then
! Find xs_listing and set kT
index_list = xs_listing_dict % get_key(to_lower(list_names % tail % data))
if(xs_listings(index_list) % kT /= 0.0_8) then
mat % sqrtkT = sqrt(xs_listings(index_list) % kT * 1.0D6)
temp_known = .true.
end if
end if
! Set ace or iso-in-lab scattering for each nuclide in element
do k = 1, n_nuc_ele
if (adjustl(to_lower(temp_str)) == "iso-in-lab") then
@ -2177,12 +2179,6 @@ contains
end do NATURAL_ELEMENTS
! If still undefined, set the temperature to zero
if (.not. temp_known) then
mat % sqrtkT = 0.0_8
temp_known = .true.
end if
! ========================================================================
! COPY NUCLIDES TO ARRAYS IN MATERIAL

1
src/material_header.F90
View file

@ -13,7 +13,6 @@ module material_header
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
real(8) :: sqrtkT ! sqrt(kT), kT in eV
! Energy grid information
integer :: n_grid ! # of union material grid points

6
src/particle_header.F90
View file

@ -2,7 +2,7 @@ module particle_header
use bank_header, only: Bank
use constants, only: NEUTRON, ONE, NONE, ZERO, MAX_SECONDARY, &
MAX_DELAYED_GROUPS
MAX_DELAYED_GROUPS, ERROR_REAL
use error, only: fatal_error
use geometry_header, only: BASE_UNIVERSE
@ -79,6 +79,9 @@ module particle_header
integer :: material ! index for current material
integer :: last_material ! index for last material
! Temperature of the current cell
real(8) :: sqrtkT ! sqrt(k_Boltzmann * temperature) in MeV
! Statistical data
integer :: n_collision ! # of collisions
@ -124,6 +127,7 @@ contains
this % absorb_wgt = ZERO
this % n_bank = 0
this % wgt_bank = ZERO
this % sqrtkT = ERROR_REAL
this % n_collision = 0
this % fission = .false.
this % delayed_group = 0