diff --git a/src/ace.F90 b/src/ace.F90 index adff30dbeb..58f7fd5707 100644 --- a/src/ace.F90 +++ b/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 diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 27cf647155..2067030c1e 100644 --- a/src/ace_header.F90 +++ b/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 !=============================================================================== diff --git a/src/cross_section.F90 b/src/cross_section.F90 index cdb4b1e2cc..18a23a6b0f 100644 --- a/src/cross_section.F90 +++ b/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 diff --git a/src/initialize.F90 b/src/initialize.F90 index 8d465ef7d7..3e67ae8662 100644 --- a/src/initialize.F90 +++ b/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()