Merge pull request #466 from cjosey/improve-xs-perf

Improvements to XS calculation performance

src/cross_section.F90

@@ -33,6 +33,7 @@ contains
     integer :: i_nuclide     ! index into nuclides array
     integer :: i_sab         ! index into sab_tables array
     integer :: j             ! index in mat % i_sab_nuclides
+    integer :: u             ! index into logarithmic mapping array
     real(8) :: atom_density  ! atom density of a nuclide
     logical :: check_sab     ! should we check for S(a,b) table?
     type(Material), pointer :: mat ! current material
@@ -50,9 +51,13 @@ contains
 
     mat => materials(p % material)
 
-    ! Find energy index on global or material unionized grid
-    if (grid_method == GRID_MAT_UNION) &
-         call find_energy_index(p % E, p % material)
+    ! Find energy index on energy grid
+    u = 0
+    if (grid_method == GRID_MAT_UNION) then
+      call find_energy_index(p % E, p % material)
+    else if (grid_method == GRID_LOGARITHM) then
+      u = int(log(p % E/1.0e-11_8)/log_spacing)
+    end if
 
     ! Determine if this material has S(a,b) tables
     check_sab = (mat % n_sab > 0)
@@ -94,9 +99,9 @@ contains
 
       ! Calculate microscopic cross section for this nuclide
       if (p % E /= micro_xs(i_nuclide) % last_E) then
-        call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i)
+        call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, u)
       else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then
-        call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i)
+        call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, u)
       end if
 
       ! ========================================================================
@@ -137,16 +142,16 @@ contains
 ! given index in the nuclides array at the energy of the given particle
 !===============================================================================
 
-  subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_mat, i_nuc_mat)
+  subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_mat, i_nuc_mat, u)
 
     integer, intent(in) :: i_nuclide ! index into nuclides array
     integer, intent(in) :: i_sab     ! index into sab_tables array
     integer, intent(in) :: i_mat     ! index into materials array
     integer, intent(in) :: i_nuc_mat ! index into nuclides array for a material
+    integer, intent(in) :: u         ! index into logarithmic mapping array
     integer :: i_grid ! index on nuclide energy grid
     integer :: i_low  ! lower logarithmic mapping index
     integer :: i_high ! upper logarithmic mapping index
-    integer :: u      ! index into logarithmic mapping array
     real(8), intent(in) :: E ! energy
     real(8) :: f             ! interp factor on nuclide energy grid
     type(Nuclide),  pointer :: nuc
@@ -173,7 +178,6 @@ contains
       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
 

src/search.F90

@@ -28,7 +28,6 @@ contains
     integer :: L
     integer :: R
     integer :: n_iteration
-    real(8) :: testval
 
     L = 1
     R = n
@@ -39,22 +38,11 @@ contains
 
     n_iteration = 0
     do while (R - L > 1)
-
-      ! Check boundaries
-      if (val > array(L) .and. val < array(L+1)) then
-        array_index = L
-        return
-      elseif (val > array(R-1) .and. val < array(R)) then
-        array_index = R - 1
-        return
-      end if
-
       ! Find values at midpoint
       array_index = L + (R - L)/2
-      testval = array(array_index)
-      if (val >= testval) then
+      if (val >= array(array_index)) then
         L = array_index
-      elseif (val < testval) then
+      else
         R = array_index
       end if
 
@@ -80,7 +68,6 @@ contains
     integer :: L
     integer :: R
     integer :: n_iteration
-    real(8) :: testval
 
     L = 1
     R = n
@@ -91,22 +78,11 @@ contains
 
     n_iteration = 0
     do while (R - L > 1)
-
-      ! Check boundaries
-      if (val > array(L) .and. val < array(L+1)) then
-        array_index = L
-        return
-      elseif (val > array(R-1) .and. val < array(R)) then
-        array_index = R - 1
-        return
-      end if
-
       ! Find values at midpoint
       array_index = L + (R - L)/2
-      testval = array(array_index)
-      if (val >= testval) then
+      if (val >= array(array_index)) then
         L = array_index
-      elseif (val < testval) then
+      else
         R = array_index
       end if
 
@@ -132,7 +108,6 @@ contains
     integer :: L
     integer :: R
     integer :: n_iteration
-    real(8) :: testval
 
     L = 1
     R = n
@@ -143,22 +118,11 @@ contains
 
     n_iteration = 0
     do while (R - L > 1)
-
-      ! Check boundaries
-      if (val > array(L) .and. val < array(L+1)) then
-        array_index = L
-        return
-      elseif (val > array(R-1) .and. val < array(R)) then
-        array_index = R - 1
-        return
-      end if
-
       ! Find values at midpoint
       array_index = L + (R - L)/2
-      testval = array(array_index)
-      if (val >= testval) then
+      if (val >= array(array_index)) then
         L = array_index
-      elseif (val < testval) then
+      else
         R = array_index
       end if