diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 615908dc46..d4167bcb60 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -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 diff --git a/src/OBJECTS b/src/OBJECTS index 8f15589b5e..3746d3fa3d 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -14,6 +14,7 @@ geometry.o \ geometry_header.o \ global.o \ initialize.o \ +interpolation.o \ input_xml.o \ logging.o \ main.o \ diff --git a/src/interpolation.f90 b/src/interpolation.f90 new file mode 100644 index 0000000000..c0e6f631b6 --- /dev/null +++ b/src/interpolation.f90 @@ -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 diff --git a/src/physics.f90 b/src/physics.f90 index 1dd2ef8543..98ef890ea2 100644 --- a/src/physics.f90 +++ b/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