Add new energy grid search based on an equal log-spaced mapping technique as

outlined in LA-UR-14-24530.

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

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 ::        &

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

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

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

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

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()