Merge pull request #941 from brbass/energy_bins

Fix energy bin sampling in source

src/input_xml.F90

@@ -4183,7 +4183,6 @@ contains
     integer :: n_libraries
     logical :: file_exists ! does mgxs.h5 exist?
     integer(HID_T) :: file_id
-    real(8), allocatable :: rev_energy_bins(:)
     character(len=MAX_WORD_LEN), allocatable :: names(:)
 
     ! Check if MGXS Library exists
@@ -4224,13 +4223,19 @@ contains
     end if
 
     ! First reverse the order of energy_groups
+    rev_energy_bins = energy_bins
     energy_bins = energy_bins(num_energy_groups + 1:1:-1)
 
+    ! Get the midpoint of the energy groups
     allocate(energy_bin_avg(num_energy_groups))
     do i = 1, num_energy_groups
       energy_bin_avg(i) = HALF * (energy_bins(i) + energy_bins(i + 1))
     end do
 
+    ! Get the minimum and maximum energies
+    energy_min_neutron = energy_bins(num_energy_groups + 1)
+    energy_max_neutron = energy_bins(1)
+
     ! Get the datasets present in the library
     call get_groups(file_id, names)
     n_libraries = size(names)

src/mgxs_header.F90

@@ -221,12 +221,15 @@ module mgxs_header
   ! Number of delayed groups
   integer :: num_delayed_groups
 
-  ! Energy group structure
+  ! Energy group structure with decreasing energy
   real(8), allocatable :: energy_bins(:)
 
   ! Midpoint of the energy group structure
   real(8), allocatable :: energy_bin_avg(:)
 
+  ! Energy group structure with increasing energy
+  real(8), allocatable :: rev_energy_bins(:)
+
 contains
 
 !===============================================================================

src/source.F90

@@ -15,7 +15,7 @@ module source
   use hdf5_interface,   only: file_create, file_open, file_close, read_dataset
   use math
   use message_passing,  only: rank
-  use mgxs_header,      only: energy_bins, num_energy_groups
+  use mgxs_header,      only: rev_energy_bins, num_energy_groups
   use output,           only: write_message
   use particle_header,  only: Particle
   use random_lcg,       only: prn, set_particle_seed, prn_set_stream
@@ -129,14 +129,9 @@ contains
 
     ! If running in MG, convert site % E to group
     if (.not. run_CE) then
-      if (site % E <= energy_bins(1)) then
-        site % E = real(1, 8)
-      else if (site % E > energy_bins(num_energy_groups + 1)) then
-        site % E = real(num_energy_groups, 8)
-      else
-        site % E = real(binary_search(energy_bins, num_energy_groups + 1, &
-             site % E), 8)
-      end if
+      site % E = real(binary_search(rev_energy_bins, num_energy_groups + 1, &
+           site % E), 8)
+      site % E = num_energy_groups + 1 - site % E
     end if
 
     ! Set the random number generator back to the tracking stream.

src/source_header.F90

@@ -9,7 +9,7 @@ module source_header
   use error
   use geometry, only: find_cell
   use material_header, only: materials
-  use nuclide_header, only: energy_max_neutron
+  use nuclide_header, only: energy_min_neutron, energy_max_neutron
   use particle_header, only: Particle
   use string, only: to_lower
   use xml_interface
@@ -275,9 +275,12 @@ contains
     ! Check for monoenergetic source above maximum neutron energy
     select type (energy => this % energy)
     type is (Discrete)
-      if (any(energy % x >= energy_max_neutron)) then
+      if (any(energy % x > energy_max_neutron)) then
         call fatal_error("Source energy above range of energies of at least &
              &one cross section table")
+      else if (any(energy % x < energy_min_neutron)) then
+        call fatal_error("Source energy below range of energies of at least &
+             &one cross section table")
       end if
     end select
 
@@ -285,8 +288,8 @@ contains
       ! Sample energy spectrum
       site % E = this % energy % sample()
 
-      ! resample if energy is greater than maximum neutron energy
-      if (site % E < energy_max_neutron) exit
+      ! Resample if energy falls outside minimum or maximum neutron energy
+      if (site % E < energy_max_neutron .and. site % E > energy_min_neutron) exit
     end do
 
     ! Set delayed group

src/tallies/tally_filter_energy.F90

@@ -8,7 +8,7 @@ module tally_filter_energy
   use constants
   use error
   use hdf5_interface
-  use mgxs_header,         only: num_energy_groups, energy_bins
+  use mgxs_header,         only: num_energy_groups, rev_energy_bins
   use particle_header,     only: Particle
   use settings,            only: run_CE
   use string,              only: to_str
@@ -75,7 +75,7 @@ contains
     ! ordering of the library and tallying systems).
     if (.not. run_CE) then
       if (n == num_energy_groups + 1) then
-        if (all(this % bins == energy_bins(num_energy_groups + 1:1:-1))) &
+        if (all(this % bins == rev_energy_bins)) &
              then
           this % matches_transport_groups = .true.
         end if