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

only use nuclide grids for resonance scattering

Browse source
This commit is contained in:
walshjon 2014-09-10 12:03:38 -04:00
parent 20cffb0c4a
commit 4f8634680e
5 changed files with 52 additions and 59 deletions
Download patch file Download diff file Expand all files Collapse all files

4
src/ace_header.F90
View file

@ -117,7 +117,6 @@ module ace_header
character(10) :: name_0K ! name of 0K nuclide, e.g. 92235.00c
character(16) :: scheme ! target velocity sampling scheme
integer :: n_grid_0K ! number of 0K energy grid points
integer, allocatable :: grid_index_0K(:) ! pointers to union grid
real(8), allocatable :: energy_0K(:) ! energy grid for 0K xs
real(8), allocatable :: elastic_0K(:) ! Microscopic elastic cross section
real(8), allocatable :: xs_cdf(:) ! CDF of v_rel times cross section
@ -372,9 +371,6 @@ module ace_header
if (allocated(this % grid_index)) &
deallocate(this % grid_index)
if (allocated(this % grid_index_0K)) &
deallocate(this % grid_index_0K)
if (allocated(this % energy)) &
deallocate(this % energy, this % total, this % elastic, &

39
src/cross_section.F90
View file

@ -527,37 +527,26 @@ contains
real(8) :: xs_out ! 0K xs at trial energy
! Determine index on nuclide energy grid
select case (grid_method)
case (GRID_UNION)
! If we're using the unionized grid with pointers, finding the index on
! the nuclide energy grid is as simple as looking up the pointer
call find_energy_index(E)
i_grid = nuc % grid_index_0K(union_grid_index)
case (GRID_NUCLIDE)
! If we're not using the unionized grid, we have to do a binary search on
! the nuclide energy grid in order to determine which points to
! interpolate between
if (E < nuc % energy(1)) then
i_grid = 1
elseif (E > nuc % energy(nuc % n_grid)) then
i_grid = nuc % n_grid - 1
else
i_grid = binary_search(nuc % energy, nuc % n_grid, E)
end if
end select
if (E < nuc % energy_0K(1)) then
i_grid = 1
elseif (E > nuc % energy_0K(nuc % n_grid_0K)) then
i_grid = nuc % n_grid_0K - 1
else
i_grid = binary_search(nuc % energy_0K, nuc % n_grid_0K, E)
end if
! check for rare case where two energy points are the same
if (nuc % energy_0K(i_grid) == nuc % energy_0K(i_grid+1)) i_grid = i_grid + 1
if (nuc % energy_0K(i_grid) == nuc % energy_0K(i_grid+1)) then
i_grid = i_grid + 1
end if
! calculate interpolation factor
f = (E - nuc%energy_0K(i_grid))/(nuc%energy_0K(i_grid+1) - nuc%energy_0K(i_grid))
f = (E - nuc % energy_0K(i_grid)) &
& / (nuc % energy_0K(i_grid + 1) - nuc % energy_0K(i_grid))
! Calculate microscopic nuclide elastic cross section
xs_out = (ONE - f) * nuc % elastic_0K(i_grid) + f * nuc % elastic_0K(i_grid+1)
xs_out = (ONE - f) * nuc % elastic_0K(i_grid) &
& + f * nuc % elastic_0K(i_grid + 1)
end function elastic_xs_0K

18
src/energy_grid.F90
View file

@ -147,24 +147,6 @@ contains
end if
nuc % grid_index(j) = index_e - 1
end do
! set pointers for 0K energy grid to the unionized grid
if (nuc % resonant) then
allocate(nuc % grid_index_0K(n_grid))
index_e = 1
energy = nuc % energy_0K(index_e)
do j = 1, n_grid
union_energy = e_grid(j)
if (union_energy >= energy .and. index_e < nuc % n_grid_0K) then
index_e = index_e + 1
energy = nuc % energy_0K(index_e)
end if
nuc % grid_index_0K(j) = index_e - 1
end do
end if
end do
end subroutine grid_pointers

4
src/input_xml.F90
View file

@ -842,7 +842,7 @@ contains
end if
! check that E_min is non-negative
if (E_min < ZERO) then
if (nuclides_0K(i) % E_min < ZERO) then
message = "Lower resonance scattering energy bound is negative"
call fatal_error()
end if
@ -853,7 +853,7 @@ contains
end if
! check that E_max is not less than E_min
if (E_max < E_min) then
if (nuclides_0K(i) % E_max < nuclides_0K(i) % E_min) then
message = "Lower resonance scattering energy bound exceeds upper"
call fatal_error()
end if

46
src/physics.F90
View file

@ -835,13 +835,26 @@ contains
E_red = sqrt((awr * E) / kT)
E_low = (((E_red - 4.0_8)**2) * kT) / awr
E_up = (((E_red + 4.0_8)**2) * kT) / awr
! find lower and upper energy bound indices
call find_energy_index(E_low)
i_E_low = nuc % grid_index_0K(union_grid_index)
call find_energy_index(E_up)
i_E_up = nuc % grid_index_0K(union_grid_index)
! lower index
if (E_low < nuc % energy_0K(1)) then
i_E_low = 1
elseif (E_low > nuc % energy_0K(nuc % n_grid_0K)) then
i_E_low = nuc % n_grid_0K - 1
else
i_E_low = binary_search(nuc % energy_0K, nuc % n_grid_0K, E_low)
end if
! upper index
if (E_up < nuc % energy_0K(1)) then
i_E_up = 1
elseif (E_up > nuc % energy_0K(nuc % n_grid_0K)) then
i_E_up = nuc % n_grid_0K - 1
else
i_E_up = binary_search(nuc % energy_0K, nuc % n_grid_0K, E_up)
end if
! interpolate xs since we're not exactly at the energy indices
xs_low = nuc % elastic_0K(i_E_low)
m = (nuc % elastic_0K(i_E_low + 1) - xs_low) &
@ -878,10 +891,23 @@ contains
E_up = (((E_red + 4.0_8)**2) * kT) / awr
! find lower and upper energy bound indices
call find_energy_index(E_low)
i_E_low = nuc % grid_index_0K(union_grid_index)
call find_energy_index(E_up)
i_E_up = nuc % grid_index_0K(union_grid_index)
! lower index
if (E_low < nuc % energy_0K(1)) then
i_E_low = 1
elseif (E_low > nuc % energy_0K(nuc % n_grid_0K)) then
i_E_low = nuc % n_grid_0K - 1
else
i_E_low = binary_search(nuc % energy_0K, nuc % n_grid_0K, E_low)
end if
! upper index
if (E_up < nuc % energy_0K(1)) then
i_E_up = 1
elseif (E_up > nuc % energy_0K(nuc % n_grid_0K)) then
i_E_up = nuc % n_grid_0K - 1
else
i_E_up = binary_search(nuc % energy_0K, nuc % n_grid_0K, E_up)
end if
! interpolate xs CDF since we're not exactly at the energy indices
! cdf value at lower bound attainable energy