diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 27cf647155..11cb109a93 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -97,7 +97,7 @@ module ace_header ! Energy grid information integer :: n_grid ! # of nuclide grid points - integer, allocatable :: grid_index(:) ! pointers to union grid + integer, allocatable :: grid_index(:) ! union grid pointers / log grid mapping real(8), allocatable :: energy(:) ! energy values corresponding to xs ! Microscopic cross sections diff --git a/src/constants.F90 b/src/constants.F90 index 94edd9d2a2..7bee73b1af 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -382,9 +382,10 @@ module constants ! Energy grid methods integer, parameter :: & - GRID_NUCLIDE = 1, & ! non-unionized energy grid - GRID_UNION = 2, & ! union grid with pointers - GRID_LETHARGY = 3 ! lethargy mapping + GRID_NUCLIDE = 1, & ! non-unionized energy grid + GRID_UNION = 2, & ! union grid with pointers + GRID_LOGARITHM = 3 ! logarithmic mapping + integer, parameter :: N_LOG_BINS = 8000 ! Running modes integer, parameter :: & diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 7426519cda..a150670724 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -141,8 +141,10 @@ contains integer, intent(in) :: i_sab ! index into sab_tables array real(8), intent(in) :: E ! energy - integer :: i_grid ! index on nuclide energy grid - real(8) :: f ! interp factor on nuclide energy grid + integer :: i_grid ! index on nuclide energy grid + integer :: i_low, i_high ! bounding indices from logarithmic mapping + integer :: u ! index into logarithmic mapping array + real(8) :: f ! interp factor on nuclide energy grid type(Nuclide), pointer, save :: nuc => null() !$omp threadprivate(nuc) @@ -157,10 +159,29 @@ contains i_grid = nuc % grid_index(union_grid_index) + case (GRID_LOGARITHM) + ! Determine the energy grid index using a logarithmic mapping to reduce + ! the energy range over which a binary search needs to be performed + + if (E < nuc % energy(1)) then + i_grid = 1 + elseif (E > nuc % energy(nuc % n_grid)) then + i_grid = nuc % n_grid - 1 + else + ! Determine bounding indices based on which equal log-spaced interval + ! the energy is in + u = int(log(E/1.0e-11_8)/log_spacing) + i_low = nuc % grid_index(u) + i_high = nuc % grid_index(u + 1) + 1 + + ! Perform binary search over reduced range + i_grid = binary_search(nuc % energy(i_low:i_high), & + i_high - i_low + 1, E) + i_low - 1 + end if + 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 + ! Perform 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 diff --git a/src/energy_grid.F90 b/src/energy_grid.F90 index 66b95be0de..934cf27534 100644 --- a/src/energy_grid.F90 +++ b/src/energy_grid.F90 @@ -1,6 +1,6 @@ module energy_grid - use constants, only: MAX_LINE_LEN + use constants, only: MAX_LINE_LEN, N_LOG_BINS use global use list_header, only: ListReal use output, only: write_message @@ -148,4 +148,53 @@ contains end subroutine grid_pointers +!=============================================================================== +! LOGARITHMIC_GRID determines a logarithmic mapping for energies to bounding +! indices on a nuclide energy grid +!=============================================================================== + + subroutine logarithmic_grid() + + integer :: i, j, k ! Loop indices + integer :: M ! Number of equally log-spaced bins + real(8) :: E_max ! Maximum energy in MeV + real(8) :: E_min ! Minimum energy in MeV + real(8), allocatable :: umesh(:) ! Equally log-spaced energy grid + type(Nuclide), pointer :: nuc => null() + + ! Set minimum/maximum energies + E_max = 20.0_8 + E_min = 1.0e-11_8 + + ! Determine equal-logarithmic energy spacing + M = N_LOG_BINS + log_spacing = log(E_max/E_min)/M + + ! Create equally log-spaced energy grid + allocate(umesh(0:M)) + umesh(:) = [(i*log_spacing, i=0, M)] + + do i = 1, n_nuclides_total + ! Allocate logarithmic mapping for nuclide + nuc => nuclides(i) + allocate(nuc % grid_index(0:M)) + + ! Determine corresponding indices in nuclide grid to energies on + ! equal-logarithmic grid + j = 1 + do k = 0, M - 1 + do while (log(nuc%energy(j + 1)/E_min) <= umesh(k)) + j = j + 1 + end do + nuc % grid_index(k) = j + end do + + ! Set the last point explicitly so that we don't have out-of-bounds issues + nuc % grid_index(M) = size(nuc % energy) - 1 + end do + + deallocate(umesh) + + end subroutine logarithmic_grid + end module energy_grid diff --git a/src/global.F90 b/src/global.F90 index 7d728a25d6..33674eb2b1 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -80,6 +80,7 @@ module global ! Unionized energy grid integer :: grid_method ! how to treat the energy grid integer :: n_grid ! number of points on unionized grid + real(8) :: log_spacing ! spacing on logarithmic grid real(8), allocatable :: e_grid(:) ! energies on unionized grid ! Unreoslved resonance probablity tables diff --git a/src/initialize.F90 b/src/initialize.F90 index c14a8a02a7..19586a5e5e 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -4,7 +4,7 @@ module initialize use bank_header, only: Bank use constants use dict_header, only: DictIntInt, ElemKeyValueII - use energy_grid, only: unionized_grid + use energy_grid, only: unionized_grid, logarithmic_grid use error, only: fatal_error, warning use geometry, only: neighbor_lists use geometry_header, only: Cell, Universe, Lattice, BASE_UNIVERSE @@ -109,6 +109,8 @@ contains call time_unionize % start() call unionized_grid() call time_unionize % stop() + elseif (grid_method == GRID_LOGARITHM) then + call logarithmic_grid() end if ! Allocate and setup tally stride, matching_bins, and tally maps diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 46461b78da..45b6067cb0 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -220,9 +220,8 @@ contains grid_method = GRID_NUCLIDE case ('union') grid_method = GRID_UNION - case ('lethargy') - message = "Lethargy mapped energy grid not yet supported." - call fatal_error() + case ('logarithm', 'logarithmic', 'log') + grid_method = GRID_LOGARITHM case default message = "Unknown energy grid method: " // trim(temp_str) call fatal_error()