Merge pull request #282 from walshjon/ptable-resample

probability table resample

src/ace.F90

@@ -7,6 +7,7 @@ module ace
   use error,            only: fatal_error, warning
   use fission,          only: nu_total
   use global
+  use list_header,      only: ListElemInt
   use material_header,  only: Material
   use output,           only: write_message
   use set_header,       only: SetChar
@@ -1428,4 +1429,31 @@ contains
 
   end function get_real
 
+!===============================================================================
+! SAME_NUCLIDE_LIST creates a linked list for each nuclide containing the
+! indices in the nuclides array of all other instances of that nuclide.  For
+! example, the same nuclide may exist at multiple temperatures resulting
+! in multiple entries in the nuclides array for a single zaid number.
+!===============================================================================
+
+  subroutine same_nuclide_list()
+
+    integer :: i ! index in nuclides array
+    integer :: j ! index in nuclides array
+    type(ListElemInt), pointer :: nuc_list => null() ! pointer to nuclide list
+
+    do i = 1, n_nuclides_total
+      allocate(nuclides(i) % nuc_list)
+      nuc_list => nuclides(i) % nuc_list
+      do j = 1, n_nuclides_total
+        if (nuclides(i) % zaid == nuclides(j) % zaid) then
+          nuc_list % data = j
+          allocate(nuc_list % next)
+          nuc_list => nuc_list % next
+        end if
+      end do
+    end do
+
+  end subroutine same_nuclide_list
+
 end module ace

src/ace_header.F90

@@ -2,6 +2,7 @@ module ace_header
 
   use constants,   only: MAX_FILE_LEN
   use endf_header, only: Tab1
+  use list_header, only: ListElemInt
 
   implicit none
 
@@ -95,6 +96,9 @@ module ace_header
     real(8)       :: awr     ! weight of nucleus in neutron masses
     real(8)       :: kT      ! temperature in MeV (k*T)
 
+    ! Linked list of indices in nuclides array of instances of this same nuclide
+    type(ListElemInt), pointer :: nuc_list => null()
+
     ! Energy grid information
     integer :: n_grid                     ! # of nuclide grid points
     integer, allocatable :: grid_index(:) ! pointers to union grid
@@ -241,6 +245,7 @@ module ace_header
 
     ! Information for URR probability table use
     logical :: use_ptable  ! in URR range with probability tables?
+    real(8) :: last_prn
   end type NuclideMicroXS
 
 !===============================================================================

src/cross_section.F90

@@ -5,6 +5,7 @@ module cross_section
   use error,           only: fatal_error
   use fission,         only: nu_total
   use global
+  use list_header,     only: ListElemInt
   use material_header, only: Material
   use particle_header, only: Particle
   use random_lcg,      only: prn
@@ -184,6 +185,8 @@ contains
     ! Initialize sab treatment to false
     micro_xs(i_nuclide) % index_sab   = NONE
     micro_xs(i_nuclide) % elastic_sab = ZERO
+
+    ! Initialize URR probability table treatment to false
     micro_xs(i_nuclide) % use_ptable  = .false.
 
     ! Initialize nuclide cross-sections to zero
@@ -348,10 +351,12 @@ contains
     real(8) :: capture    ! (n,gamma) cross section
     real(8) :: fission    ! fission cross section
     real(8) :: inelastic  ! inelastic cross section
-    type(UrrData),  pointer, save :: urr => null()
-    type(Nuclide),  pointer, save :: nuc => null()
-    type(Reaction), pointer, save :: rxn => null()
-!$omp threadprivate(urr, nuc, rxn)
+    logical :: same_nuc   ! do we know the xs for this nuclide at this energy?
+    type(UrrData),  pointer, save :: urr      => null()
+    type(Nuclide),  pointer, save :: nuc      => null()
+    type(Reaction), pointer, save :: rxn      => null()
+    type(ListElemInt), pointer    :: nuc_list => null()
+!$omp threadprivate(urr, nuc, rxn, nuc_list)
 
     micro_xs(i_nuclide) % use_ptable = .true.
 
@@ -371,7 +376,28 @@ contains
          (urr % energy(i_energy + 1) - urr % energy(i_energy))
 
     ! sample probability table using the cumulative distribution
-    r = prn()
+
+    ! if we're dealing with a nuclide that we've previously encountered at
+    ! this energy but a different temperature, use the original random number to
+    ! preserve correlation of temperature in probability tables
+    same_nuc = .false.
+    nuc_list => nuc % nuc_list
+    do
+      if (E /= ZERO .and. E == micro_xs(nuc_list % data) % last_E) then
+        same_nuc = .true.
+        exit
+      end if
+      nuc_list => nuc_list % next
+      if (.not. associated(nuc_list % next)) exit
+    end do
+
+    if (same_nuc) then
+      r = micro_xs(nuc_list % data) % last_prn
+    else
+      r = prn()
+      micro_xs(i_nuclide) % last_prn = r
+    end if
+
     i_low = 1
     do
       if (urr % prob(i_energy, URR_CUM_PROB, i_low) > r) exit

src/initialize.F90

@@ -1,6 +1,6 @@
 module initialize
 
-  use ace,              only: read_xs
+  use ace,              only: read_xs, same_nuclide_list
   use bank_header,      only: Bank
   use constants
   use dict_header,      only: DictIntInt, ElemKeyValueII
@@ -104,6 +104,9 @@ contains
       call read_xs()
       call time_read_xs % stop()
 
+      ! Create linked lists for multiple instances of the same nuclide
+      call same_nuclide_list()
+
       ! Construct unionized energy grid from cross-sections
       if (grid_method == GRID_UNION) then
         call time_unionize % start()