From fc2681d65fcb79849a0f4cb38e49d10ce0dd30c8 Mon Sep 17 00:00:00 2001 From: "Jonathan A. Walsh" Date: Thu, 12 Jun 2014 15:30:38 -0600 Subject: [PATCH 1/7] don't resample ptable xs if diff temperature but same E --- src/ace_header.F90 | 1 + src/cross_section.F90 | 28 ++++++++++++++++++++++------ src/geometry.F90 | 1 + src/material_header.F90 | 1 + 4 files changed, 25 insertions(+), 6 deletions(-) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index e5bf1bb3a4..532a355d3f 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -241,6 +241,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 7426519cda..80995da861 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -92,9 +92,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) + call calculate_nuclide_xs(i_nuclide, i_sab, p % E, mat % first_E) else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then - call calculate_nuclide_xs(i_nuclide, i_sab, p % E) + call calculate_nuclide_xs(i_nuclide, i_sab, p % E, mat % first_E) end if ! ======================================================================== @@ -135,11 +135,12 @@ contains ! given index in the nuclides array at the energy of the given particle !=============================================================================== - subroutine calculate_nuclide_xs(i_nuclide, i_sab, E) + subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, first_E) integer, intent(in) :: i_nuclide ! index into nuclides array integer, intent(in) :: i_sab ! index into sab_tables array real(8), intent(in) :: E ! energy + real(8), intent(in) :: first_E ! particle energy when it entered this mat integer :: i_grid ! index on nuclide energy grid real(8) :: f ! interp factor on nuclide energy grid @@ -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 @@ -233,7 +236,7 @@ contains if (urr_ptables_on .and. nuc % urr_present) then if (E > nuc % urr_data % energy(1) .and. & E < nuc % urr_data % energy(nuc % urr_data % n_energy)) then - call calculate_urr_xs(i_nuclide, E) + call calculate_urr_xs(i_nuclide, E, first_E) end if end if @@ -334,10 +337,11 @@ contains ! from probability tables !=============================================================================== - subroutine calculate_urr_xs(i_nuclide, E) + subroutine calculate_urr_xs(i_nuclide, E, first_E) integer, intent(in) :: i_nuclide ! index into nuclides array real(8), intent(in) :: E ! energy + real(8), intent(in) :: first_E ! particle energy when it entered this mat integer :: i_energy ! index for energy integer :: i_table ! index for table @@ -370,8 +374,20 @@ contains (urr % energy(i_energy + 1) - urr % energy(i_energy)) ! sample probability table using the cumulative distribution - r = prn() + ! (if we're dealing with an isotope 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) + if (E /= ZERO .and. E == first_E .and. & + & minval(dble(abs(nuclides(:) % zaid - nuclides(i_nuclide) % zaid)) & + & + abs(micro_xs(:) % last_E - E)) == ZERO) then + r = micro_xs(minloc(dble(abs(nuclides(:) % zaid & + & - nuclides(i_nuclide) % zaid)) + abs(micro_xs(:) % last_E & + & - E),1)) % last_prn + else + r = prn() + end if + i_table = 1 + do if (urr % prob(i_energy, URR_CUM_PROB, i_table) > r) exit i_table = i_table + 1 diff --git a/src/geometry.F90 b/src/geometry.F90 index 24d3ef8b06..1a2544d28d 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -190,6 +190,7 @@ contains ! set material p % last_material = p % material p % material = c % material + materials(p % material) % first_E = p % E elseif (c % type == CELL_FILL) then ! ==================================================================== diff --git a/src/material_header.F90 b/src/material_header.F90 index cbfd917498..d4c7e81290 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -12,6 +12,7 @@ module material_header integer, allocatable :: nuclide(:) ! index in nuclides array real(8) :: density ! total atom density in atom/b-cm real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm + real(8) :: first_E ! particle energy upon entering mat ! S(a,b) data references integer :: n_sab = 0 ! number of S(a,b) tables From ad3519a0efd8a49679e60486f0b2528196130fe9 Mon Sep 17 00:00:00 2001 From: "Jonathan A. Walsh" Date: Thu, 12 Jun 2014 16:48:24 -0600 Subject: [PATCH 2/7] set last prn of a nuclide whenever a prn is needed --- src/cross_section.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 80995da861..2f7bcf9e3a 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -384,6 +384,7 @@ contains & - E),1)) % last_prn else r = prn() + micro_xs(i_nuclide) % last_prn = r end if i_table = 1 From 3d38155179383ab2ab284dc7aa3a393e7e07874a Mon Sep 17 00:00:00 2001 From: Jon Walsh Date: Wed, 30 Jul 2014 19:06:31 -0600 Subject: [PATCH 3/7] fixed comment that was over 80 characters --- src/cross_section.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 2f7bcf9e3a..815e6e90d7 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -374,8 +374,10 @@ contains (urr % energy(i_energy + 1) - urr % energy(i_energy)) ! sample probability table using the cumulative distribution + ! (if we're dealing with an isotope 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) + ! this energy but a different temperature, use the original random number to + ! preserve correlation of temperature in probability tables) if (E /= ZERO .and. E == first_E .and. & & minval(dble(abs(nuclides(:) % zaid - nuclides(i_nuclide) % zaid)) & & + abs(micro_xs(:) % last_E - E)) == ZERO) then From 560e979836afb8d0cb135af565da09f235c16d0f Mon Sep 17 00:00:00 2001 From: Jon Walsh Date: Wed, 30 Jul 2014 22:10:49 -0600 Subject: [PATCH 4/7] check for void material to avoid out of bounds --- src/geometry.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 1a2544d28d..03eddb5ec4 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -190,7 +190,8 @@ contains ! set material p % last_material = p % material p % material = c % material - materials(p % material) % first_E = p % E + if(p % material /= MATERIAL_VOID) & + & materials(p % material) % first_E = p % E elseif (c % type == CELL_FILL) then ! ==================================================================== From 60127faddca151c328f397012c118ae2afd13760 Mon Sep 17 00:00:00 2001 From: Jon Walsh Date: Thu, 31 Jul 2014 21:22:14 -0600 Subject: [PATCH 5/7] removed vestigial check, variable --- src/cross_section.F90 | 14 ++++++-------- src/geometry.F90 | 2 -- src/material_header.F90 | 1 - 3 files changed, 6 insertions(+), 11 deletions(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 815e6e90d7..cb97e5512a 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -92,9 +92,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, mat % first_E) + call calculate_nuclide_xs(i_nuclide, i_sab, p % E) else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then - call calculate_nuclide_xs(i_nuclide, i_sab, p % E, mat % first_E) + call calculate_nuclide_xs(i_nuclide, i_sab, p % E) end if ! ======================================================================== @@ -135,12 +135,11 @@ contains ! given index in the nuclides array at the energy of the given particle !=============================================================================== - subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, first_E) + subroutine calculate_nuclide_xs(i_nuclide, i_sab, E) integer, intent(in) :: i_nuclide ! index into nuclides array integer, intent(in) :: i_sab ! index into sab_tables array real(8), intent(in) :: E ! energy - real(8), intent(in) :: first_E ! particle energy when it entered this mat integer :: i_grid ! index on nuclide energy grid real(8) :: f ! interp factor on nuclide energy grid @@ -236,7 +235,7 @@ contains if (urr_ptables_on .and. nuc % urr_present) then if (E > nuc % urr_data % energy(1) .and. & E < nuc % urr_data % energy(nuc % urr_data % n_energy)) then - call calculate_urr_xs(i_nuclide, E, first_E) + call calculate_urr_xs(i_nuclide, E) end if end if @@ -337,11 +336,10 @@ contains ! from probability tables !=============================================================================== - subroutine calculate_urr_xs(i_nuclide, E, first_E) + subroutine calculate_urr_xs(i_nuclide, E) integer, intent(in) :: i_nuclide ! index into nuclides array real(8), intent(in) :: E ! energy - real(8), intent(in) :: first_E ! particle energy when it entered this mat integer :: i_energy ! index for energy integer :: i_table ! index for table @@ -378,7 +376,7 @@ contains ! (if we're dealing with an isotope 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) - if (E /= ZERO .and. E == first_E .and. & + if (E /= ZERO .and. & & minval(dble(abs(nuclides(:) % zaid - nuclides(i_nuclide) % zaid)) & & + abs(micro_xs(:) % last_E - E)) == ZERO) then r = micro_xs(minloc(dble(abs(nuclides(:) % zaid & diff --git a/src/geometry.F90 b/src/geometry.F90 index 03eddb5ec4..24d3ef8b06 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -190,8 +190,6 @@ contains ! set material p % last_material = p % material p % material = c % material - if(p % material /= MATERIAL_VOID) & - & materials(p % material) % first_E = p % E elseif (c % type == CELL_FILL) then ! ==================================================================== diff --git a/src/material_header.F90 b/src/material_header.F90 index d4c7e81290..cbfd917498 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -12,7 +12,6 @@ module material_header integer, allocatable :: nuclide(:) ! index in nuclides array real(8) :: density ! total atom density in atom/b-cm real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm - real(8) :: first_E ! particle energy upon entering mat ! S(a,b) data references integer :: n_sab = 0 ! number of S(a,b) tables From e6e4898045adc982d4a7e4b89d0818bb8eb7107b Mon Sep 17 00:00:00 2001 From: Jon Walsh Date: Sat, 9 Aug 2014 20:18:38 -0600 Subject: [PATCH 6/7] use linked list to find multiple instances of a single nuclide --- src/ace.F90 | 28 ++++++++++++++++++++++++++++ src/ace_header.F90 | 4 ++++ src/cross_section.F90 | 32 +++++++++++++++++++++----------- src/initialize.F90 | 5 ++++- 4 files changed, 57 insertions(+), 12 deletions(-) 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 fe645683be..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 diff --git a/src/cross_section.F90 b/src/cross_section.F90 index cb97e5512a..cc810dc5e6 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 @@ -349,9 +350,11 @@ 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() + 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) micro_xs(i_nuclide) % use_ptable = .true. @@ -373,15 +376,22 @@ contains ! sample probability table using the cumulative distribution - ! (if we're dealing with an isotope that we've previously encountered at + ! 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) - if (E /= ZERO .and. & - & minval(dble(abs(nuclides(:) % zaid - nuclides(i_nuclide) % zaid)) & - & + abs(micro_xs(:) % last_E - E)) == ZERO) then - r = micro_xs(minloc(dble(abs(nuclides(:) % zaid & - & - nuclides(i_nuclide) % zaid)) + abs(micro_xs(:) % last_E & - & - E),1)) % last_prn + ! 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 diff --git a/src/initialize.F90 b/src/initialize.F90 index c14a8a02a7..a5c64101f9 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() From 251f71deeb7b7947802345a5b8c53e646cad96f8 Mon Sep 17 00:00:00 2001 From: Jon Walsh Date: Sat, 9 Aug 2014 21:26:22 -0600 Subject: [PATCH 7/7] make the linked list threadprivate --- src/cross_section.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 23df55acd5..18a23a6b0f 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -356,7 +356,7 @@ contains type(Nuclide), pointer, save :: nuc => null() type(Reaction), pointer, save :: rxn => null() type(ListElemInt), pointer :: nuc_list => null() -!$omp threadprivate(urr, nuc, rxn) +!$omp threadprivate(urr, nuc, rxn, nuc_list) micro_xs(i_nuclide) % use_ptable = .true.