mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 06:05:58 -04:00
Added interpolation module with capabilities for all ENDF interpolation schemes. Changed laws 7,9,11.
This commit is contained in:
Paul Romano
parent
51d540273a
commit
d5bc1a7d80
4 changed files with 130 additions and 101 deletions
|
|
@ -84,6 +84,10 @@ input_xml.o: xml_materials.o
|
|||
input_xml.o: xml_settings.o
|
||||
input_xml.o: xml_tallies.o
|
||||
|
||||
interpolation.o: constants.o
|
||||
interpolation.o: error.o
|
||||
interpolation.o: search.o
|
||||
|
||||
logging.o: constants.o
|
||||
logging.o: global.o
|
||||
|
||||
|
|
@ -122,6 +126,7 @@ physics.o: fission.o
|
|||
physics.o: geometry.o
|
||||
physics.o: geometry_header.o
|
||||
physics.o: global.o
|
||||
physics.o: interpolation.o
|
||||
physics.o: mcnp_random.o
|
||||
physics.o: output.o
|
||||
physics.o: particle_header.o
|
||||
|
|
|
|||
|
|
@ -14,6 +14,7 @@ geometry.o \
|
|||
geometry_header.o \
|
||||
global.o \
|
||||
initialize.o \
|
||||
interpolation.o \
|
||||
input_xml.o \
|
||||
logging.o \
|
||||
main.o \
|
||||
|
|
|
|||
113
src/interpolation.f90
Normal file
113
src/interpolation.f90
Normal file
|
|
@ -0,0 +1,113 @@
|
|||
module interpolation
|
||||
|
||||
use constants
|
||||
use error, only: fatal_error
|
||||
use search, only: binary_search
|
||||
|
||||
implicit none
|
||||
|
||||
contains
|
||||
|
||||
!===============================================================================
|
||||
! INTERPOLATE_TAB1 interpolates a function between two points based on
|
||||
! particular interpolation scheme. The data needs to be organized as a ENDF TAB1
|
||||
! type function containing the interpolation regions, break points, and
|
||||
! tabulated x's and y's.
|
||||
!===============================================================================
|
||||
|
||||
function interpolate_tab1(data, x, loc_start) result(y)
|
||||
|
||||
real(8), intent(in) :: data(:) ! array of data whose
|
||||
real(8), intent(in) :: x ! x value to find y at
|
||||
integer, intent(in), optional :: loc_start ! starting location in data
|
||||
real(8) :: y ! y(x)
|
||||
|
||||
integer :: i ! bin in which to interpolate
|
||||
integer :: loc_0 ! starting location
|
||||
integer :: n_regions ! number of interpolation regions
|
||||
integer :: n_points ! number of tabulated values
|
||||
integer :: interp ! ENDF interpolation scheme
|
||||
integer :: loc_breakpoints ! location of breakpoints in data
|
||||
integer :: loc_interp ! location of interpolation schemes in data
|
||||
integer :: loc_x ! location of x's in data
|
||||
integer :: loc_y ! location of y's in data
|
||||
real(8) :: r ! interpolation factor
|
||||
real(8) :: x0, x1 ! bounding x values
|
||||
real(8) :: y0, y1 ! bounding y values
|
||||
character(MAX_LINE_LEN) :: msg
|
||||
|
||||
! determine starting location
|
||||
if (present(loc_start)) then
|
||||
loc_0 = loc_start - 1
|
||||
else
|
||||
loc_0 = 0
|
||||
end if
|
||||
|
||||
! determine number of interpolation regions
|
||||
n_regions = data(loc_0 + 1)
|
||||
|
||||
! set locations for breakpoints and interpolation schemes
|
||||
loc_breakpoints = loc_0 + 1
|
||||
loc_interp = loc_breakpoints + n_regions
|
||||
|
||||
! determine number of tabulated values
|
||||
n_points = data(loc_interp + n_regions + 1)
|
||||
|
||||
! set locations for x's and y's
|
||||
loc_x = loc_interp + n_regions + 1
|
||||
loc_y = loc_x + n_points
|
||||
|
||||
! find which bin the abscissa is in -- if the abscissa is outside the
|
||||
! tabulated range, the first or last point is chosen, i.e. no interpolation
|
||||
! is done outside the energy range
|
||||
if (x < data(loc_x + 1)) then
|
||||
y = data(loc_y + 1)
|
||||
return
|
||||
elseif (x > data(loc_x + n_points)) then
|
||||
y = data(loc_y + n_points)
|
||||
return
|
||||
else
|
||||
i = binary_search(data(loc_x + 1:loc_x + n_points), n_points, x)
|
||||
end if
|
||||
|
||||
! determine interpolation scheme
|
||||
if (n_regions == 0) then
|
||||
interp = LINEAR_LINEAR
|
||||
elseif (n_regions == 1) then
|
||||
interp = data(loc_interp + 1)
|
||||
elseif (n_regions > 1) then
|
||||
msg = "Multiple interpolation regions not yet supported."
|
||||
call fatal_error(msg)
|
||||
end if
|
||||
|
||||
! handle special case of histogram interpolation
|
||||
if (interp == HISTOGRAM) then
|
||||
y = data(loc_y + i)
|
||||
return
|
||||
end if
|
||||
|
||||
! determine bounding values
|
||||
x0 = data(loc_x + i)
|
||||
x1 = data(loc_x + i + 1)
|
||||
y0 = data(loc_y + i)
|
||||
y1 = data(loc_y + i + 1)
|
||||
|
||||
! determine interpolation factor and interpolated value
|
||||
select case (interp)
|
||||
case (LINEAR_LINEAR)
|
||||
r = (x - x0)/(x1 - x0)
|
||||
y = (1 - r)*y0 + r*y1
|
||||
case (LINEAR_LOG)
|
||||
r = (log(x) - log(x0))/(log(x1) - log(x0))
|
||||
y = (1 - r)*y0 + r*y1
|
||||
case (LOG_LINEAR)
|
||||
r = (x - x0)/(x1 - x0)
|
||||
y = exp((1-r)*log(y0) + r*log(y1))
|
||||
case (LOG_LOG)
|
||||
r = (log(x) - log(x0))/(log(x1) - log(x0))
|
||||
y = exp((1-r)*log(y0) + r*log(y1))
|
||||
end select
|
||||
|
||||
end function interpolate_tab1
|
||||
|
||||
end module interpolation
|
||||
112
src/physics.f90
112
src/physics.f90
|
|
@ -9,6 +9,7 @@ module physics
|
|||
cross_lattice
|
||||
use geometry_header, only: Universe, BASE_UNIVERSE
|
||||
use global
|
||||
use interpolation, only: interpolate_tab1
|
||||
use mcnp_random, only: rang
|
||||
use output, only: message, print_particle
|
||||
use particle_header, only: Particle
|
||||
|
|
@ -1115,36 +1116,13 @@ contains
|
|||
! =======================================================================
|
||||
! MAXWELL FISSION SPECTRUM
|
||||
|
||||
! read number of interpolation regions and incoming energies
|
||||
NR = edist % data(1)
|
||||
NE = edist % data(2 + 2*NR)
|
||||
if (NR > 0) then
|
||||
msg = "Multiple interpolation regions not supported while &
|
||||
&attempting to sample Maxwell fission spectrum."
|
||||
call fatal_error(msg)
|
||||
end if
|
||||
|
||||
! find incident energy bin and calculate interpolation factor
|
||||
loc = 2 + 2*NR
|
||||
if (E_in < edist % data(loc+1)) then
|
||||
i = 1
|
||||
r = ZERO
|
||||
elseif (E_in > edist % data(loc+NE)) then
|
||||
i = NE - 1
|
||||
r = ONE
|
||||
else
|
||||
i = binary_search(edist % data(loc+1), NE, E_in)
|
||||
r = (E_in - edist%data(loc+i)) / &
|
||||
& (edist%data(loc+i+1) - edist%data(loc+i))
|
||||
end if
|
||||
|
||||
! determine nuclear temperature from tabulated function
|
||||
loc = loc + NE
|
||||
T = edist%data(loc+i) + r * &
|
||||
& (edist%data(loc+i+1) - edist%data(loc+i))
|
||||
T = interpolate_tab1(edist % data, E_in)
|
||||
|
||||
! sample maxwell fission spectrum
|
||||
E_out = maxwell_spectrum(T)
|
||||
|
||||
! TODO: Add restriction energy constraint??
|
||||
|
||||
case (9)
|
||||
! =======================================================================
|
||||
|
|
@ -1153,35 +1131,12 @@ contains
|
|||
! read number of interpolation regions and incoming energies
|
||||
NR = edist % data(1)
|
||||
NE = edist % data(2 + 2*NR)
|
||||
if (NR > 0) then
|
||||
msg = "Multiple interpolation regions not supported while &
|
||||
&attempting to sample evaporation spectrum."
|
||||
call fatal_error(msg)
|
||||
end if
|
||||
|
||||
! find energy bin and calculate interpolation factor -- if the energy is
|
||||
! outside the range of the tabulated energies, choose the first or last
|
||||
! bins
|
||||
loc = 2 + 2*NR
|
||||
if (E_in < edist % data(loc+1)) then
|
||||
i = 1
|
||||
r = ZERO
|
||||
elseif (E_in > edist % data(loc+NE)) then
|
||||
i = NE - 1
|
||||
r = ONE
|
||||
else
|
||||
i = binary_search(edist % data(loc+1), NE, E_in)
|
||||
r = (E_in - edist%data(loc+i)) / &
|
||||
& (edist%data(loc+i+1) - edist%data(loc+i))
|
||||
end if
|
||||
|
||||
! determine nuclear temperature from tabulated function
|
||||
loc = loc + NE
|
||||
T = edist%data(loc+i) + r * &
|
||||
& (edist%data(loc+i+1) - edist%data(loc+i))
|
||||
T = interpolate_tab1(edist % data, E_in)
|
||||
|
||||
! determine restriction energy
|
||||
loc = loc + NE
|
||||
loc = 2 + 2*NR + 2*NE
|
||||
U = edist % data(loc + 1)
|
||||
|
||||
! sample outgoing energy based on evaporation spectrum probability
|
||||
|
|
@ -1201,64 +1156,19 @@ contains
|
|||
! parameter 'a'
|
||||
NR = edist % data(1)
|
||||
NE = edist % data(2 + 2*NR)
|
||||
if (NR > 0) then
|
||||
msg = "Multiple interpolation regions not supported while &
|
||||
&attempting to sample Watt fission spectrum."
|
||||
call fatal_error(msg)
|
||||
end if
|
||||
|
||||
! find incident energy bin and calculate interpolation factor
|
||||
loc = 2 + 2*NR
|
||||
if (E_in < edist % data(loc+1)) then
|
||||
i = 1
|
||||
r = ZERO
|
||||
elseif (E_in > edist % data(loc+NE)) then
|
||||
i = NE - 1
|
||||
r = ONE
|
||||
else
|
||||
i = binary_search(edist % data(loc+1), NE, E_in)
|
||||
r = (E_in - edist%data(loc+i)) / &
|
||||
& (edist%data(loc+i+1) - edist%data(loc+i))
|
||||
end if
|
||||
|
||||
! determine Watt parameter 'a' from tabulated function
|
||||
loc = loc + NE
|
||||
Watt_a = edist%data(loc+i) + r * &
|
||||
& (edist%data(loc+i+1) - edist%data(loc+i))
|
||||
|
||||
! read number of interpolation regions and incoming energies for
|
||||
! parameter 'b'
|
||||
loc = loc + NE
|
||||
NR = edist % data(loc + 1)
|
||||
NE = edist % data(loc + 2 + 2*NR)
|
||||
if (NR > 0) then
|
||||
msg = "Multiple interpolation regions not supported while &
|
||||
&attempting to sample Watt fission spectrum."
|
||||
call fatal_error(msg)
|
||||
end if
|
||||
|
||||
! find incident energy bin and calculate interpolation factor
|
||||
loc = loc + 2 + 2*NR
|
||||
if (E_in < edist % data(loc+1)) then
|
||||
i = 1
|
||||
r = ZERO
|
||||
elseif (E_in > edist % data(loc+NE)) then
|
||||
i = NE - 1
|
||||
r = ONE
|
||||
else
|
||||
i = binary_search(edist % data(loc+1), NE, E_in)
|
||||
r = (E_in - edist%data(loc+i)) / &
|
||||
& (edist%data(loc+i+1) - edist%data(loc+i))
|
||||
end if
|
||||
Watt_a = interpolate_tab1(edist % data, E_in)
|
||||
|
||||
! determine Watt parameter 'b' from tabulated function
|
||||
loc = loc + NE
|
||||
Watt_b = edist%data(loc+i) + r * &
|
||||
& (edist%data(loc+i+1) - edist%data(loc+i))
|
||||
loc = 3 + 2*(NR + NE)
|
||||
Watt_b = interpolate_tab1(edist % data, E_in, loc)
|
||||
|
||||
! Sample energy-dependent Watt fission spectrum
|
||||
E_out = watt_spectrum(Watt_a, Watt_b)
|
||||
|
||||
! TODO: Add restriction energy constraint??
|
||||
|
||||
case (44)
|
||||
! =======================================================================
|
||||
! KALBACH-MANN CORRELATED SCATTERING
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue