From 7f83000812b04149903fdb22cac1b00b8b549257 Mon Sep 17 00:00:00 2001 From: Jose Salcedo Perez Date: Mon, 17 Jul 2017 17:56:27 +0000 Subject: [PATCH 01/14] Adding direct-mapping implementation. --- src/input_xml.F90 | 10 ++++++++++ src/material_header.F90 | 1 + src/tallies/tally.F90 | 25 ++++++++----------------- 3 files changed, 19 insertions(+), 17 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 91cd57ae29..8dcd9335fb 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2450,6 +2450,16 @@ contains n_nuclides = index_nuclide n_sab_tables = index_sab + + do i=1, n_materials + mat => materials(i) + allocate(mat % mat_nuclide_list(n_nuclides_total)) + mat % mat_nuclide_list(:) = 0 + do j=1, mat % n_nuclides + mat % mat_nuclide_list(mat % nuclide(j))=j + end do + end do + ! Close materials XML file call doc % clear() diff --git a/src/material_header.F90 b/src/material_header.F90 index 9d93131acf..4f4b66e780 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -32,6 +32,7 @@ module material_header character(len=104) :: name = "" ! User-defined name integer :: n_nuclides = 0 ! number of nuclides integer, allocatable :: nuclide(:) ! index in nuclides array + integer, allocatable :: mat_nuclide_list(:) real(8) :: density ! total atom density in atom/b-cm real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm real(8) :: density_gpcc ! total density in g/cm^3 diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 8b1e7ca9a2..fbbccfde11 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -2730,7 +2730,7 @@ contains type(Particle), intent(in) :: p real(8), intent(in) :: distance - + integer :: check_value integer :: i integer :: i_tally integer :: i_filt @@ -2744,7 +2744,7 @@ contains real(8) :: filter_weight ! combined weight of all filters logical :: finished ! found all valid bin combinations type(Material), pointer :: mat - + ! Determine track-length estimate of flux flux = p % wgt * distance @@ -2811,26 +2811,17 @@ contains if (p % material /= MATERIAL_VOID) then ! Get pointer to current material mat => materials(p % material) + check_value= mat % mat_nuclide_list(i_nuclide) + + if (check_value == 0) then + cycle NUCLIDE_BIN_LOOP + end if - ! Determine if nuclide is actually in material - NUCLIDE_MAT_LOOP: do j = 1, mat % n_nuclides - ! If index of nuclide matches the j-th nuclide listed in the - ! material, break out of the loop - if (i_nuclide == mat % nuclide(j)) exit - - ! If we've reached the last nuclide in the material, it means - ! the specified nuclide to be tallied is not in this material - if (j == mat % n_nuclides) then - cycle NUCLIDE_BIN_LOOP - end if - end do NUCLIDE_MAT_LOOP - - atom_density = mat % atom_density(j) + atom_density = mat % atom_density(check_value) else atom_density = ZERO end if end if - ! Determine score for each bin call score_general(p, t, (k-1)*t % n_score_bins, filter_index, & i_nuclide, atom_density, flux * filter_weight) From 4e40788128aca5bbd1e108443081a039ad9782e7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 21 Dec 2017 17:16:26 +0700 Subject: [PATCH 02/14] Fix up material nuclide mapping --- src/input_xml.F90 | 13 ++++---- src/material_header.F90 | 7 +++- src/tallies/tally.F90 | 72 ++++++++++++++++------------------------- 3 files changed, 41 insertions(+), 51 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 8dcd9335fb..363c088672 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2450,13 +2450,14 @@ contains n_nuclides = index_nuclide n_sab_tables = index_sab - - do i=1, n_materials + ! Create direct address tables for tally optimization purposes + do i = 1, n_materials mat => materials(i) - allocate(mat % mat_nuclide_list(n_nuclides_total)) - mat % mat_nuclide_list(:) = 0 - do j=1, mat % n_nuclides - mat % mat_nuclide_list(mat % nuclide(j))=j + allocate(mat % mat_nuclide_index(n_nuclides)) + mat % mat_nuclide_index(:) = 0 + + do j = 1, mat % n_nuclides + mat % mat_nuclide_index(mat % nuclide(j)) = j end do end do diff --git a/src/material_header.F90 b/src/material_header.F90 index 4f4b66e780..3287e796f2 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -32,11 +32,16 @@ module material_header character(len=104) :: name = "" ! User-defined name integer :: n_nuclides = 0 ! number of nuclides integer, allocatable :: nuclide(:) ! index in nuclides array - integer, allocatable :: mat_nuclide_list(:) real(8) :: density ! total atom density in atom/b-cm real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm real(8) :: density_gpcc ! total density in g/cm^3 + ! To improve performance of tallying, we store an array (direct address + ! table) that indicates for each nuclide in the global nuclides(:) array the + ! index of the corresponding nuclide in the Material % nuclide(:) array. If + ! it is not present in the material, the entry is set to zero. + integer, allocatable :: mat_nuclide_index(:) + ! Energy grid information integer :: n_grid ! # of union material grid points real(8), allocatable :: e_grid(:) ! union material grid energies diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index fbbccfde11..eaee96fa7e 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -2377,10 +2377,6 @@ contains i_tally = active_analog_tallies % data(i) associate (t => tallies(i_tally) % obj) - ! Get pointer to current material. We need this in order to determine what - ! nuclides are in the material - mat => materials(p % material) - ! Find all valid bins in each filter if they have not already been found ! for a previous tally. do j = 1, size(t % filter) @@ -2423,33 +2419,28 @@ contains ! Nuclide logic ! Check for nuclide bins - k = 0 - NUCLIDE_LOOP: do while (k < t % n_nuclide_bins) - - ! Increment the index in the list of nuclide bins - k = k + 1 - + NUCLIDE_LOOP: do k = 1, t % n_nuclide_bins + ! Get index of nuclide in nuclides array i_nuclide = t % nuclide_bins(k) if (i_nuclide > 0) then - atom_density = -ONE - ! Check to see if this nuclide was in the material of our collision - do m = 1, mat % n_nuclides - if (mat % nuclide(m) == i_nuclide) then - atom_density = mat % atom_density(m) - exit - end if - end do + if (p % material /= MATERIAL_VOID) then + ! Get pointer to current material + mat => materials(p % material) + + ! Determine index of nuclide in Material % atom_density array + j = mat % mat_nuclide_index(i_nuclide) + if (j == 0) cycle NUCLIDE_LOOP + + ! Copy corresponding atom density + atom_density = mat % atom_density(j) + end if else atom_density = ZERO end if - ! If we found the nuclide, determine the score for each bin - if (atom_density >= ZERO) then - call score_general(p, t, (k-1)*t % n_score_bins, filter_index, & - i_nuclide, atom_density, filter_weight) - end if - + call score_general(p, t, (k-1)*t % n_score_bins, filter_index, & + i_nuclide, atom_density, filter_weight) end do NUCLIDE_LOOP ! ====================================================================== @@ -2730,7 +2721,7 @@ contains type(Particle), intent(in) :: p real(8), intent(in) :: distance - integer :: check_value + integer :: i integer :: i_tally integer :: i_filt @@ -2744,7 +2735,7 @@ contains real(8) :: filter_weight ! combined weight of all filters logical :: finished ! found all valid bin combinations type(Material), pointer :: mat - + ! Determine track-length estimate of flux flux = p % wgt * distance @@ -2811,17 +2802,18 @@ contains if (p % material /= MATERIAL_VOID) then ! Get pointer to current material mat => materials(p % material) - check_value= mat % mat_nuclide_list(i_nuclide) - - if (check_value == 0) then - cycle NUCLIDE_BIN_LOOP - end if - atom_density = mat % atom_density(check_value) + ! Determine index of nuclide in Material % atom_density array + j = mat % mat_nuclide_index(i_nuclide) + if (j == 0) cycle NUCLIDE_BIN_LOOP + + ! Copy corresponding atom density + atom_density = mat % atom_density(j) else atom_density = ZERO end if end if + ! Determine score for each bin call score_general(p, t, (k-1)*t % n_score_bins, filter_index, & i_nuclide, atom_density, flux * filter_weight) @@ -2968,19 +2960,11 @@ contains ! Get pointer to current material mat => materials(p % material) - ! Determine if nuclide is actually in material - NUCLIDE_MAT_LOOP: do j = 1, mat % n_nuclides - ! If index of nuclide matches the j-th nuclide listed in the - ! material, break out of the loop - if (i_nuclide == mat % nuclide(j)) exit - - ! If we've reached the last nuclide in the material, it means - ! the specified nuclide to be tallied is not in this material - if (j == mat % n_nuclides) then - cycle NUCLIDE_BIN_LOOP - end if - end do NUCLIDE_MAT_LOOP + ! Determine index of nuclide in Material % atom_density array + j = mat % mat_nuclide_index(i_nuclide) + if (j == 0) cycle NUCLIDE_BIN_LOOP + ! Copy corresponding atom density atom_density = mat % atom_density(j) else atom_density = ZERO From 29a8737b25e78e5f2cdb3f9a87c85bff2f4a941b Mon Sep 17 00:00:00 2001 From: Jose Salcedo Perez Date: Tue, 15 Aug 2017 14:46:04 +0000 Subject: [PATCH 03/14] Use direct address table for reactions in nuclides --- src/nuclide_header.F90 | 10 ++++++---- src/tallies/tally.F90 | 14 +++++++------- 2 files changed, 13 insertions(+), 11 deletions(-) diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 676e1b500d..1831f0a892 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -87,9 +87,9 @@ module nuclide_header ! Reactions type(Reaction), allocatable :: reactions(:) - type(DictIntInt) :: reaction_index ! map MT values to index in reactions + !type(DictIntInt) :: reaction_index ! map MT values to index in reactions ! array; used at tally-time - + integer, allocatable :: rxn_index_MT(:) ! Fission energy release class(Function1D), allocatable :: fission_q_prompt ! prompt neutrons, gammas class(Function1D), allocatable :: fission_q_recov ! neutrons, gammas, betas @@ -561,7 +561,8 @@ contains n_temperature = size(this % kTs) allocate(this % sum_xs(n_temperature)) - + allocate(this % rxn_index_MT(1000)) + this % rxn_index_MT(:)=0 do i = 1, n_temperature ! Allocate and initialize derived cross sections n_grid = size(this % grid(i) % energy) @@ -580,8 +581,9 @@ contains i_fission = 0 do i = 1, size(this % reactions) + call MTs % push_back(this % reactions(i) % MT) - call this % reaction_index % set(this % reactions(i) % MT, i) + this % rxn_index_MT(this % reactions(i) % MT) = i associate (rx => this % reactions(i)) ! Skip total inelastic level scattering, gas production cross sections diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index eaee96fa7e..c47f94aa76 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -247,7 +247,7 @@ contains ! multiplicities of one. score = p % last_wgt * flux else - m = nuclides(p % event_nuclide) % reaction_index % get(p % event_MT) + m = nuclides(p % event_nuclide) % rxn_index_MT(p % event_MT) ! Get yield and apply to score associate (rxn => nuclides(p % event_nuclide) % reactions(m)) @@ -273,7 +273,7 @@ contains ! multiplicities of one. score = p % last_wgt * flux else - m = nuclides(p % event_nuclide) % reaction_index % get(p % event_MT) + m = nuclides(p % event_nuclide) % rxn_index_MT(p % event_MT) ! Get yield and apply to score associate (rxn => nuclides(p % event_nuclide) % reactions(m)) @@ -299,7 +299,7 @@ contains ! multiplicities of one. score = p % last_wgt * flux else - m = nuclides(p % event_nuclide) % reaction_index % get(p % event_MT) + m = nuclides(p % event_nuclide) % rxn_index_MT(p % event_MT) ! Get yield and apply to score associate (rxn => nuclides(p%event_nuclide)%reactions(m)) @@ -1148,8 +1148,8 @@ contains score = ZERO if (i_nuclide > 0) then - m = nuclides(i_nuclide) % reaction_index % get(score_bin) - if (m /= EMPTY) then + m = nuclides(i_nuclide) % rxn_index_MT(score_bin) + if (m /= 0) then ! Retrieve temperature and energy grid index and interpolation ! factor i_temp = micro_xs(i_nuclide) % index_temp @@ -1187,8 +1187,8 @@ contains ! Get index in nuclides array i_nuc = materials(p % material) % nuclide(l) - m = nuclides(i_nuc) % reaction_index % get(score_bin) - if (m /= EMPTY) then + m = nuclides(i_nuc) % rxn_index_MT(score_bin) + if (m /= 0) then ! Retrieve temperature and energy grid index and ! interpolation factor i_temp = micro_xs(i_nuc) % index_temp From fe19b9b6dd7fec5ce5ebae806419d5dca0f18cd3 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 22 Dec 2017 13:11:17 +0700 Subject: [PATCH 04/14] Make rxn_index_MT allocated as part of Nuclide type --- src/nuclide_header.F90 | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 1831f0a892..2c5de1b140 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -87,9 +87,11 @@ module nuclide_header ! Reactions type(Reaction), allocatable :: reactions(:) - !type(DictIntInt) :: reaction_index ! map MT values to index in reactions - ! array; used at tally-time - integer, allocatable :: rxn_index_MT(:) + + ! Array that maps MT values to index in reactions; used at tally-time. Note + ! that ENDF-102 does not assign any MT values above 891. + integer :: rxn_index_MT(891) + ! Fission energy release class(Function1D), allocatable :: fission_q_prompt ! prompt neutrons, gammas class(Function1D), allocatable :: fission_q_recov ! neutrons, gammas, betas @@ -561,8 +563,7 @@ contains n_temperature = size(this % kTs) allocate(this % sum_xs(n_temperature)) - allocate(this % rxn_index_MT(1000)) - this % rxn_index_MT(:)=0 + this % rxn_index_MT(:) = 0 do i = 1, n_temperature ! Allocate and initialize derived cross sections n_grid = size(this % grid(i) % energy) From b10b7204b2cf17ae1a3eb58f1cdc58732f6c488d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 22 Dec 2017 13:37:15 +0700 Subject: [PATCH 05/14] Create mat_nuclide_index as part of simulation_init --- src/input_xml.F90 | 11 ---------- src/material_header.F90 | 47 ++++++++++++++++++++++++++++++----------- src/simulation.F90 | 19 ++++++++++++----- 3 files changed, 49 insertions(+), 28 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 363c088672..91cd57ae29 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2450,17 +2450,6 @@ contains n_nuclides = index_nuclide n_sab_tables = index_sab - ! Create direct address tables for tally optimization purposes - do i = 1, n_materials - mat => materials(i) - allocate(mat % mat_nuclide_index(n_nuclides)) - mat % mat_nuclide_index(:) = 0 - - do j = 1, mat % n_nuclides - mat % mat_nuclide_index(mat % nuclide(j)) = j - end do - end do - ! Close materials XML file call doc % clear() diff --git a/src/material_header.F90 b/src/material_header.F90 index 3287e796f2..d772204156 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -69,6 +69,7 @@ module material_header contains procedure :: set_density => material_set_density + procedure :: init_nuclide_index => material_init_nuclide_index procedure :: assign_sab_tables => material_assign_sab_tables end type Material @@ -85,8 +86,8 @@ contains ! MATERIAL_SET_DENSITY sets the total density of a material in atom/b-cm. !=============================================================================== - function material_set_density(m, density) result(err) - class(Material), intent(inout) :: m + function material_set_density(this, density) result(err) + class(Material), intent(inout) :: this real(8), intent(in) :: density integer :: err @@ -94,23 +95,23 @@ contains real(8) :: sum_percent real(8) :: awr - if (allocated(m % atom_density)) then + if (allocated(this % atom_density)) then ! Set total density based on value provided - m % density = density + this % density = density ! Determine normalized atom percents - sum_percent = sum(m % atom_density) - m % atom_density(:) = m % atom_density / sum_percent + sum_percent = sum(this % atom_density) + this % atom_density(:) = this % atom_density / sum_percent ! Recalculate nuclide atom densities based on given density - m % atom_density(:) = density * m % atom_density + this % atom_density(:) = density * this % atom_density ! Calculate density in g/cm^3. - m % density_gpcc = ZERO - do i = 1, m % n_nuclides - awr = nuclides(m % nuclide(i)) % awr - m % density_gpcc = m % density_gpcc & - + m % atom_density(i) * awr * MASS_NEUTRON / N_AVOGADRO + this % density_gpcc = ZERO + do i = 1, this % n_nuclides + awr = nuclides(this % nuclide(i)) % awr + this % density_gpcc = this % density_gpcc & + + this % atom_density(i) * awr * MASS_NEUTRON / N_AVOGADRO end do err = 0 else @@ -119,6 +120,28 @@ contains end if end function material_set_density +!=============================================================================== +! INIT_NUCLIDE_INDEX creates a mapping from indices in the global nuclides(:) +! array to the Material % nuclides array +!=============================================================================== + + subroutine material_init_nuclide_index(this) + class(Material), intent(inout) :: this + + integer :: i + + ! Allocate nuclide index array and set to zeros + if (allocated(this % mat_nuclide_index)) & + deallocate(this % mat_nuclide_index) + allocate(this % mat_nuclide_index(n_nuclides)) + this % mat_nuclide_index(:) = 0 + + ! Assign entries in the index array + do i = 1, this % n_nuclides + this % mat_nuclide_index(this % nuclide(i)) = i + end do + end subroutine material_init_nuclide_index + !=============================================================================== ! ASSIGN_SAB_TABLES assigns S(alpha,beta) tables to specific nuclides within ! materials so the code knows when to apply bound thermal scattering data diff --git a/src/simulation.F90 b/src/simulation.F90 index 5eacbd40b8..ade209e361 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -18,6 +18,7 @@ module simulation #endif use error, only: fatal_error use geometry_header, only: n_cells + use material_header, only: n_materials, materials use message_passing use mgxs_header, only: energy_bins, energy_bin_avg use nuclide_header, only: micro_xs, n_nuclides @@ -401,6 +402,7 @@ contains !=============================================================================== subroutine openmc_simulation_init() bind(C) + integer :: i ! Skip if simulation has already been initialized if (simulation_initialized) return @@ -418,6 +420,11 @@ contains ! Allocate tally results arrays if they're not allocated yet call configure_tallies() + ! Set up material nuclide index mapping + do i = 1, n_materials + call materials(i) % init_nuclide_index() + end do + !$omp parallel ! Allocate array for microscopic cross section cache allocate(micro_xs(n_nuclides)) @@ -459,8 +466,8 @@ contains subroutine openmc_simulation_finalize() bind(C) + integer :: i ! loop index #ifdef MPI - integer :: i ! loop index for tallies integer :: n ! size of arrays integer :: mpi_err ! MPI error code integer(8) :: temp @@ -470,9 +477,14 @@ contains ! Skip if simulation was never run if (.not. simulation_initialized) return - ! Stop active batch timer + ! Stop active batch timer and start finalization timer call time_active % stop() + call time_finalize % start() + ! Free up simulation-specific memory + do i = 1, n_materials + deallocate(materials(i) % mat_nuclide_index) + end do !$omp parallel deallocate(micro_xs, filter_matches) !$omp end parallel @@ -480,9 +492,6 @@ contains ! Increment total number of generations total_gen = total_gen + current_batch*gen_per_batch - ! Start finalization timer - call time_finalize % start() - #ifdef MPI ! Broadcast tally results so that each process has access to results if (allocated(tallies)) then From aa5d093a592f641720671b9e832f799fd78a4262 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 22 Dec 2017 13:40:24 +0700 Subject: [PATCH 06/14] Remove unused components of the Material derived type --- src/cross_section.F90 | 22 ---------------------- src/material_header.F90 | 7 ------- 2 files changed, 29 deletions(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 0bdadd87ee..46cbb93abc 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -573,28 +573,6 @@ contains end subroutine calculate_urr_xs -!=============================================================================== -! FIND_ENERGY_INDEX determines the index on the union energy grid at a certain -! energy -!=============================================================================== - - pure function find_energy_index(mat, E) result(i) - type(Material), intent(in) :: mat ! pointer to current material - real(8), intent(in) :: E ! energy of particle - integer :: i ! energy grid index - - ! if the energy is outside of energy grid range, set to first or last - ! index. Otherwise, do a binary search through the union energy grid. - if (E <= mat % e_grid(1)) then - i = 1 - elseif (E > mat % e_grid(mat % n_grid)) then - i = mat % n_grid - 1 - else - i = binary_search(mat % e_grid, mat % n_grid, E) - end if - - end function find_energy_index - !=============================================================================== ! MULTIPOLE_EVAL evaluates the windowed multipole equations for cross ! sections in the resolved resonance regions diff --git a/src/material_header.F90 b/src/material_header.F90 index d772204156..a3ce32db3d 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -42,13 +42,6 @@ module material_header ! it is not present in the material, the entry is set to zero. integer, allocatable :: mat_nuclide_index(:) - ! Energy grid information - integer :: n_grid ! # of union material grid points - real(8), allocatable :: e_grid(:) ! union material grid energies - - ! Unionized energy grid information - integer, allocatable :: nuclide_grid_index(:,:) ! nuclide e_grid pointers - ! S(a,b) data integer :: n_sab = 0 ! number of S(a,b) tables integer, allocatable :: i_sab_nuclides(:) ! index of corresponding nuclide From f774c29254313b026120ba34b31bcdba5ee629ce Mon Sep 17 00:00:00 2001 From: Jose Salcedo Perez Date: Sat, 2 Sep 2017 19:15:19 +0000 Subject: [PATCH 07/14] Cache cross sections for depletion-related reactions (ngamma, n2n, etc.) --- src/constants.F90 | 11 +- src/cross_section.F90 | 111 ++++++++++++++++- src/endf.F90 | 36 ++++-- src/input_xml.F90 | 33 +++-- src/nuclide_header.F90 | 21 ++++ src/tallies/tally.F90 | 231 ++++++++++++++++++++++++++++++++++- src/tallies/tally_header.F90 | 1 + 7 files changed, 414 insertions(+), 30 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index 853a40c8d4..936c84499a 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -308,7 +308,7 @@ module constants ! Tally score type -- if you change these, make sure you also update the ! _SCORES dictionary in openmc/capi/tally.py - integer, parameter :: N_SCORE_TYPES = 24 + integer, parameter :: N_SCORE_TYPES = 30 integer, parameter :: & SCORE_FLUX = -1, & ! flux SCORE_TOTAL = -2, & ! total reaction rate @@ -333,8 +333,13 @@ module constants SCORE_INVERSE_VELOCITY = -21, & ! flux-weighted inverse velocity SCORE_FISS_Q_PROMPT = -22, & ! prompt fission Q-value SCORE_FISS_Q_RECOV = -23, & ! recoverable fission Q-value - SCORE_DECAY_RATE = -24 ! delayed neutron precursor decay rate - + SCORE_DECAY_RATE = -24, & ! delayed neutron precursor decay rate + SCORE_N2N = -25, & + SCORE_NGAMMA = -26, & + SCORE_N3N = -27, & + SCORE_N4N = -28, & + SCORE_NP = -29, & + SCORE_NALPHA = -30 ! Maximum scattering order supported integer, parameter :: MAX_ANG_ORDER = 10 diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 46cbb93abc..2484be51b9 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -2,6 +2,7 @@ module cross_section use algorithm, only: binary_search use constants + use endf use error, only: fatal_error use list_header, only: ListElemInt use material_header, only: Material, materials @@ -45,7 +46,12 @@ contains material_xs % absorption = ZERO material_xs % fission = ZERO material_xs % nu_fission = ZERO - + material_xs % n2n = ZERO + material_xs % n3n = ZERO + material_xs % n4n = ZERO + material_xs % ngamma = ZERO + material_xs % np = ZERO + material_xs % nalpha = ZERO ! Exit subroutine if material is void if (p % material == MATERIAL_VOID) return @@ -129,6 +135,35 @@ contains ! Add contributions to material macroscopic nu-fission cross section material_xs % nu_fission = material_xs % nu_fission + & atom_density * micro_xs(i_nuclide) % nu_fission + + + ! Add contributions to material macroscopic n2n cross section + material_xs % n2n = material_xs % n2n + & + atom_density * micro_xs(i_nuclide) % n2n + + + ! Add contributions to material macroscopic ngamma cross section + material_xs % ngamma = material_xs % ngamma + & + atom_density * micro_xs(i_nuclide) % ngamma + + ! Add contributions to material macroscopic n3n cross section + material_xs % n3n = material_xs % n3n + & + atom_density * micro_xs(i_nuclide) % n3n + + ! Add contributions to material macroscopic n4n cross section + material_xs % n4n = material_xs % n4n + & + atom_density * micro_xs(i_nuclide) % n4n + + ! Add contributions to material macroscopic np cross section + material_xs % np = material_xs % np + & + atom_density * micro_xs(i_nuclide) % np + + + ! Add contributions to material macroscopic nalpha cross section + material_xs % nalpha = material_xs % nalpha + & + atom_density * micro_xs(i_nuclide) % nalpha + + end do end associate @@ -157,7 +192,9 @@ contains real(8) :: f ! interp factor on nuclide energy grid real(8) :: kT ! temperature in eV real(8) :: sigT, sigA, sigF ! Intermediate multipole variables - + integer :: rxn + type(TallyObject) :: t + character(15) :: threshold_rxn_type associate (nuc => nuclides(i_nuclide)) ! Check to see if there is multipole data present at this energy use_mp = .false. @@ -202,6 +239,8 @@ contains else ! Find the appropriate temperature index. kT = sqrtkT**2 + + !rxn = nuc % rxn_index_MT(p % event_MT) select case (temperature_method) case (TEMPERATURE_NEAREST) i_temp = minloc(abs(nuclides(i_nuclide) % kTs - kT), dim=1) @@ -247,12 +286,19 @@ contains (grid % energy(i_grid + 1) - grid % energy(i_grid)) micro_xs(i_nuclide) % index_temp = i_temp + micro_xs(i_nuclide) % index_grid = i_grid + micro_xs(i_nuclide) % interp_factor = f ! Initialize nuclide cross-sections to zero micro_xs(i_nuclide) % fission = ZERO micro_xs(i_nuclide) % nu_fission = ZERO + micro_xs(i_nuclide) % n2n = ZERO + micro_xs(i_nuclide) % n3n = ZERO + micro_xs(i_nuclide) % n4n = ZERO + micro_xs(i_nuclide) % np = ZERO + micro_xs(i_nuclide) % nalpha = ZERO micro_xs(i_nuclide) % thermal = ZERO micro_xs(i_nuclide) % thermal_elastic = ZERO @@ -274,11 +320,66 @@ contains + f * xs % fission(i_grid + 1) ! Calculate microscopic nuclide nu-fission cross section - micro_xs(i_nuclide) % nu_fission = (ONE - f) * xs % nu_fission( & + micro_xs(i_nuclide) % nu_fission = (ONE - f) * xs % nu_fission(& i_grid) + f * xs % nu_fission(i_grid + 1) + end if - end associate - end if + + micro_xs(i_nuclide) % ngamma = (micro_xs(i_nuclide) % absorption - & + micro_xs(i_nuclide) % fission) + if (t % has_threshold_rxn) then + !if (nuclides(i_nuclide)%reaction_index%has_key(p % event_MT)) the + rxn = nuclides(i_nuclide) % rxn_index_MT(p % event_MT) + associate(xs_threshold => nuc % reactions(rxn) % xs(i_temp)) + !end if + threshold_rxn_type = reaction_name(p % event_MT) + if (threshold_rxn_type == 'n2n') then + if (i_grid >= xs_threshold % threshold) then + micro_xs(i_nuclide) % n2n = (ONE - f) * xs_threshold % value(i_grid - xs_threshold & + % threshold + 1) + f * xs_threshold % value(i_grid - xs_threshold & + % threshold + 2) + end if + end if + if (threshold_rxn_type == 'n3n') then + if (i_grid >= xs_threshold % threshold) then + micro_xs(i_nuclide) % n3n = (ONE - f) * xs_threshold % value(i_grid - xs_threshold & + % threshold + 1) + f * xs_threshold % value(i_grid - xs_threshold & + % threshold + 2) + end if + end if + + + if (threshold_rxn_type == 'n4n') then + if (i_grid >= xs_threshold % threshold) then + micro_xs(i_nuclide) % n4n = (ONE - f) * xs_threshold % value(i_grid - xs_threshold & + % threshold + 1) + f * xs_threshold % value(i_grid - xs_threshold & + % threshold + 2) + end if + end if + + if (threshold_rxn_type == 'np') then + if (i_grid >= xs_threshold % threshold) then + micro_xs(i_nuclide) % np = (ONE - f) * xs_threshold % value(i_grid - xs_threshold & + % threshold + 1) + f * xs_threshold % value(i_grid - xs_threshold & + % threshold + 2) + end if + end if + + if (threshold_rxn_type == 'nalpha') then + if (i_grid >= xs_threshold % threshold) then + micro_xs(i_nuclide) % nalpha = (ONE - f) * xs_threshold % value(i_grid - xs_threshold & + % threshold + 1) + f * xs_threshold % value(i_grid - xs_threshold & + % threshold + 2) + end if + end if + + + + end associate + + end if + end associate + end if ! Initialize sab treatment to false micro_xs(i_nuclide) % index_sab = NONE diff --git a/src/endf.F90 b/src/endf.F90 index 0ba2db388b..9aace42d7b 100644 --- a/src/endf.F90 +++ b/src/endf.F90 @@ -66,6 +66,18 @@ contains string = "fission-q-prompt" case (SCORE_FISS_Q_RECOV) string = "fission-q-recoverable" + case (SCORE_N2N) + string = "n2n" + case (SCORE_N3N) + string = "n3n" + case (SCORE_N4N) + string = "n4n" + case (SCORE_NGAMMA) + string = "ngamma" + case (SCORE_NP) + string = "np" + case (SCORE_NALPHA) + string = "nalpha" ! Normal ENDF-based reactions case (TOTAL_XS) @@ -76,12 +88,12 @@ contains string = '(n,level)' case (N_2ND) string = '(n,2nd)' - case (N_2N) - string = '(n,2n)' - case (N_3N) - string = '(n,3n)' - case (N_FISSION) - string = '(n,fission)' + !case (N_2N) + !string = '(n,2n)' + !case (N_3N) + !string = '(n,3n)' + !case (N_FISSION) + ! string = '(n,fission)' case (N_F) string = '(n,f)' case (N_NF) @@ -130,18 +142,18 @@ contains string = '(n,nc)' case (N_DISAPPEAR) string = '(n,disappear)' - case (N_GAMMA) - string = '(n,gamma)' - case (N_P) - string = '(n,p)' + !case (N_GAMMA) + ! string = '(n,gamma)' + !case (N_P) + !string = '(n,p)' case (N_D) string = '(n,d)' case (N_T) string = '(n,t)' case (N_3HE) string = '(n,3He)' - case (N_A) - string = '(n,a)' + !case (N_A) + !string = '(n,a)' case (N_2A) string = '(n,2a)' case (N_3A) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 91cd57ae29..4587072f5d 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2450,6 +2450,17 @@ contains n_nuclides = index_nuclide n_sab_tables = index_sab + + do i=1, n_materials + mat => materials(i) + allocate(mat % mat_nuclide_list(n_nuclides_total)) + mat % mat_nuclide_list(:) = 0 + do j=1, mat % n_nuclides + mat % mat_nuclide_list(mat % nuclide(j))=j + end do + end do + + ! Close materials XML file call doc % clear() @@ -3047,14 +3058,14 @@ contains call fatal_error("n1n score no longer supported for tallies, & &please remove") case ('n2n', '(n,2n)') - t % score_bins(j) = N_2N - + t % score_bins(j) = SCORE_N2N + t % has_threshold_rxn = .true. case ('n3n', '(n,3n)') - t % score_bins(j) = N_3N - + t % score_bins(j) = SCORE_N3N + t % has_threshold_rxn = .true. case ('n4n', '(n,4n)') - t % score_bins(j) = N_4N - + t % score_bins(j) = SCORE_N4N + t % has_threshold_rxn = .true. case ('absorption') t % score_bins(j) = SCORE_ABSORPTION if (t % find_filter(FILTER_ENERGYOUT) > 0) then @@ -3238,9 +3249,12 @@ contains case ('(n,nc)') t % score_bins(j) = N_NC case ('(n,gamma)') - t % score_bins(j) = N_GAMMA + t % score_bins(j) = SCORE_NGAMMA case ('(n,p)') - t % score_bins(j) = N_P + t % score_bins(j) = SCORE_NP + + t % has_threshold_rxn = .true. + case ('(n,d)') t % score_bins(j) = N_D case ('(n,t)') @@ -3248,7 +3262,8 @@ contains case ('(n,3He)') t % score_bins(j) = N_3HE case ('(n,a)') - t % score_bins(j) = N_A + t % score_bins(j) = SCORE_NALPHA + t % has_threshold_rxn = .true. case ('(n,2a)') t % score_bins(j) = N_2A case ('(n,3a)') diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 2c5de1b140..1f6d45165d 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -43,6 +43,14 @@ module nuclide_header real(8), allocatable :: nu_fission(:) ! neutron production real(8), allocatable :: absorption(:) ! absorption (MT > 100) real(8), allocatable :: heating(:) ! heating + real(8), allocatable :: ngamma(:) + real(8), allocatable :: n2n(:) + real(8), allocatable :: n3n(:) + real(8), allocatable :: n4n(:) + + real(8), allocatable :: np(:) + real(8), allocatable :: nalpha(:) + end type SumXS type :: Nuclide @@ -118,6 +126,12 @@ module nuclide_header real(8) :: absorption real(8) :: fission real(8) :: nu_fission + real(8) :: n2n + real(8) :: ngamma + real(8) :: n3n + real(8) :: n4n + real(8) :: np + real(8) :: nalpha real(8) :: thermal ! Bound thermal elastic & inelastic scattering real(8) :: thermal_elastic ! Bound thermal elastic scattering @@ -148,6 +162,13 @@ module nuclide_header real(8) :: absorption ! macroscopic absorption xs real(8) :: fission ! macroscopic fission xs real(8) :: nu_fission ! macroscopic production xs + real(8) :: n2n + real(8) :: n3n + real(8) :: n4n + real(8) :: ngamma + real(8) :: np + real(8) :: nalpha + end type MaterialMacroXS !=============================================================================== diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index c47f94aa76..811ead4f8b 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -1131,6 +1131,234 @@ contains end if end if + + case (SCORE_N2N) + + score = ZERO + + if (i_nuclide > 0) then + !if (nuclides(i_nuclide)%reaction_index%has_key(score_bin)) then + !m = nuclides(i_nuclide)%reaction_index%get_key(score_bin) + + ! Retrieve temperature and energy grid index and interpolation + ! factor + !i_temp = micro_xs(i_nuclide) % index_temp + !if (i_temp > 0) then + !i_energy = micro_xs(i_nuclide) % index_grid + !f = micro_xs(i_nuclide) % interp_factor + + !associate (xs => nuclides(i_nuclide) % reactions(m) % xs(i_temp)) + !if (i_energy >= xs % threshold) then + score = micro_xs(i_nuclide) % n2n * atom_density * flux + + + !else + ! This block is reached if multipole is turned on and we're in + ! the resolved range. For (n,gamma), use absorption - + ! fission. For everything else, assume it's zero. + + ! score = ZERO + + else + do l = 1, materials(p % material) % n_nuclides + ! Get atom density + atom_density_ = materials(p % material) % atom_density(l) + + ! Get index in nuclides array + i_nuc = materials(p % material) % nuclide(l) + + score = score + micro_xs(i_nuc) % n2n * atom_density_ * flux + end do + end if + + + case (SCORE_N3N) + + score = ZERO + + if (i_nuclide > 0) then + !if (nuclides(i_nuclide)%reaction_index%has_key(score_bin)) then + !m = nuclides(i_nuclide)%reaction_index%get_key(score_bin) + + ! Retrieve temperature and energy grid index and interpolation + ! factor + !i_temp = micro_xs(i_nuclide) % index_temp + !if (i_temp > 0) then + !i_energy = micro_xs(i_nuclide) % index_grid + !f = micro_xs(i_nuclide) % interp_factor + + !associate (xs => nuclides(i_nuclide) % reactions(m) % + !xs(i_temp)) + !if (i_energy >= xs % threshold) then + score = micro_xs(i_nuclide) % n3n * atom_density * flux + + + !else + ! This block is reached if multipole is turned on and we're in + ! the resolved range. For (n,gamma), use absorption - + ! fission. For everything else, assume it's zero. + + ! score = ZERO + + else + do l = 1, materials(p % material) % n_nuclides + ! Get atom density + atom_density_ = materials(p % material) % atom_density(l) + + ! Get index in nuclides array + i_nuc = materials(p % material) % nuclide(l) + + score = score + micro_xs(i_nuc) % n3n * atom_density_ * flux + end do + end if + + + case (SCORE_N4N) + + score = ZERO + + if (i_nuclide > 0) then + !if (nuclides(i_nuclide)%reaction_index%has_key(score_bin)) then + !m = nuclides(i_nuclide)%reaction_index%get_key(score_bin) + + ! Retrieve temperature and energy grid index and interpolation + ! factor + !i_temp = micro_xs(i_nuclide) % index_temp + !if (i_temp > 0) then + !i_energy = micro_xs(i_nuclide) % index_grid + !f = micro_xs(i_nuclide) % interp_factor + + !associate (xs => nuclides(i_nuclide) % reactions(m) % + !xs(i_temp)) + !if (i_energy >= xs % threshold) then + score = micro_xs(i_nuclide) % n4n * atom_density * flux + + + !else + ! This block is reached if multipole is turned on and we're in + ! the resolved range. For (n,gamma), use absorption - + ! fission. For everything else, assume it's zero. + + ! score = ZERO + + else + do l = 1, materials(p % material) % n_nuclides + ! Get atom density + atom_density_ = materials(p % material) % atom_density(l) + + ! Get index in nuclides array + i_nuc = materials(p % material) % nuclide(l) + + score = score + micro_xs(i_nuc) % n4n * atom_density_ * flux + end do + end if + + case (SCORE_NP) + + score = ZERO + + if (i_nuclide > 0) then + !if (nuclides(i_nuclide)%reaction_index%has_key(score_bin)) then + !m = nuclides(i_nuclide)%reaction_index%get_key(score_bin) + + ! Retrieve temperature and energy grid index and interpolation + ! factor + !i_temp = micro_xs(i_nuclide) % index_temp + !if (i_temp > 0) then + !i_energy = micro_xs(i_nuclide) % index_grid + !f = micro_xs(i_nuclide) % interp_factor + + !associate (xs => nuclides(i_nuclide) % reactions(m) % + !xs(i_temp)) + !if (i_energy >= xs % threshold) then + score = micro_xs(i_nuclide) % np * atom_density * flux + + + !else + ! This block is reached if multipole is turned on and we're in + ! the resolved range. For (n,gamma), use absorption - + ! fission. For everything else, assume it's zero. + + ! score = ZERO + + else + do l = 1, materials(p % material) % n_nuclides + ! Get atom density + atom_density_ = materials(p % material) % atom_density(l) + + ! Get index in nuclides array + i_nuc = materials(p % material) % nuclide(l) + + score = score + micro_xs(i_nuc) % np * atom_density_ * flux + end do + end if + + case (SCORE_NALPHA) + + score = ZERO + + if (i_nuclide > 0) then + !if (nuclides(i_nuclide)%reaction_index%has_key(score_bin)) then + !m = nuclides(i_nuclide)%reaction_index%get_key(score_bin) + + ! Retrieve temperature and energy grid index and interpolation + ! factor + !i_temp = micro_xs(i_nuclide) % index_temp + !if (i_temp > 0) then + !i_energy = micro_xs(i_nuclide) % index_grid + !f = micro_xs(i_nuclide) % interp_factor + + !associate (xs => nuclides(i_nuclide) % reactions(m) % + !xs(i_temp)) + !if (i_energy >= xs % threshold) then + score = micro_xs(i_nuclide) % nalpha * atom_density * flux + + + !else + ! This block is reached if multipole is turned on and we're in + ! the resolved range. For (n,gamma), use absorption - + ! fission. For everything else, assume it's zero. + + ! score = ZERO + + else + do l = 1, materials(p % material) % n_nuclides + ! Get atom density + atom_density_ = materials(p % material) % atom_density(l) + + ! Get index in nuclides array + i_nuc = materials(p % material) % nuclide(l) + + score = score + micro_xs(i_nuc) % nalpha * atom_density_ * flux + end do + end if + + + + case (SCORE_NGAMMA) + if (t % estimator == ESTIMATOR_ANALOG) then + if (survival_biasing) then + ! No absorption events actually occur if survival biasing is on -- + ! just use weight absorbed in survival biasing + score = p % absorb_wgt * flux + else + ! Skip any event where the particle wasn't absorbed + if (p % event == EVENT_SCATTER) cycle SCORE_LOOP + ! All fission and absorption events will contribute here, so we + ! can just use the particle's weight entering the collision + score = p % last_wgt * flux + end if + + else + if (i_nuclide > 0) then + score = micro_xs(i_nuclide) % ngamma * atom_density * flux + else + score = material_xs % ngamma * flux + end if + end if + + + case default if (t % estimator == ESTIMATOR_ANALOG) then ! Any other score is assumed to be a MT number. Thus, we just need @@ -2721,7 +2949,7 @@ contains type(Particle), intent(in) :: p real(8), intent(in) :: distance - + integer :: check_value integer :: i integer :: i_tally integer :: i_filt @@ -2802,6 +3030,7 @@ contains if (p % material /= MATERIAL_VOID) then ! Get pointer to current material mat => materials(p % material) + check_value= mat % mat_nuclide_list(i_nuclide) ! Determine index of nuclide in Material % atom_density array j = mat % mat_nuclide_index(i_nuclide) diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index 1778474e04..c283c3f16c 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -88,6 +88,7 @@ module tally_header ! reset property - allows a tally to be reset after every batch logical :: reset = .false. + logical :: has_threshold_rxn = .false. ! Number of realizations of tally random variables integer :: n_realizations = 0 From 9ff4e620b6786e8ad0d3cf0afab15d5a3e873ddc Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 22 Dec 2017 14:47:42 +0700 Subject: [PATCH 08/14] Improvements to depletion reaction caching --- src/constants.F90 | 11 +- src/cross_section.F90 | 141 ++++++---------- src/endf.F90 | 36 ++-- src/input_xml.F90 | 33 +--- src/nuclide_header.F90 | 23 +-- src/tallies/tally.F90 | 313 +++++++++-------------------------- src/tallies/tally_header.F90 | 1 - 7 files changed, 162 insertions(+), 396 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index 936c84499a..853a40c8d4 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -308,7 +308,7 @@ module constants ! Tally score type -- if you change these, make sure you also update the ! _SCORES dictionary in openmc/capi/tally.py - integer, parameter :: N_SCORE_TYPES = 30 + integer, parameter :: N_SCORE_TYPES = 24 integer, parameter :: & SCORE_FLUX = -1, & ! flux SCORE_TOTAL = -2, & ! total reaction rate @@ -333,13 +333,8 @@ module constants SCORE_INVERSE_VELOCITY = -21, & ! flux-weighted inverse velocity SCORE_FISS_Q_PROMPT = -22, & ! prompt fission Q-value SCORE_FISS_Q_RECOV = -23, & ! recoverable fission Q-value - SCORE_DECAY_RATE = -24, & ! delayed neutron precursor decay rate - SCORE_N2N = -25, & - SCORE_NGAMMA = -26, & - SCORE_N3N = -27, & - SCORE_N4N = -28, & - SCORE_NP = -29, & - SCORE_NALPHA = -30 + SCORE_DECAY_RATE = -24 ! delayed neutron precursor decay rate + ! Maximum scattering order supported integer, parameter :: MAX_ANG_ORDER = 10 diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 2484be51b9..7066d96425 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -2,7 +2,6 @@ module cross_section use algorithm, only: binary_search use constants - use endf use error, only: fatal_error use list_header, only: ListElemInt use material_header, only: Material, materials @@ -52,6 +51,7 @@ contains material_xs % ngamma = ZERO material_xs % np = ZERO material_xs % nalpha = ZERO + ! Exit subroutine if material is void if (p % material == MATERIAL_VOID) return @@ -136,17 +136,11 @@ contains material_xs % nu_fission = material_xs % nu_fission + & atom_density * micro_xs(i_nuclide) % nu_fission - - ! Add contributions to material macroscopic n2n cross section + ! Add contributions to material macroscopic n2n cross section material_xs % n2n = material_xs % n2n + & atom_density * micro_xs(i_nuclide) % n2n - - ! Add contributions to material macroscopic ngamma cross section - material_xs % ngamma = material_xs % ngamma + & - atom_density * micro_xs(i_nuclide) % ngamma - - ! Add contributions to material macroscopic n3n cross section + ! Add contributions to material macroscopic n3n cross section material_xs % n3n = material_xs % n3n + & atom_density * micro_xs(i_nuclide) % n3n @@ -154,16 +148,17 @@ contains material_xs % n4n = material_xs % n4n + & atom_density * micro_xs(i_nuclide) % n4n - ! Add contributions to material macroscopic np cross section + ! Add contributions to material macroscopic ngamma cross section + material_xs % ngamma = material_xs % ngamma + & + atom_density * micro_xs(i_nuclide) % ngamma + + ! Add contributions to material macroscopic np cross section material_xs % np = material_xs % np + & atom_density * micro_xs(i_nuclide) % np - - ! Add contributions to material macroscopic nalpha cross section + ! Add contributions to material macroscopic nalpha cross section material_xs % nalpha = material_xs % nalpha + & atom_density * micro_xs(i_nuclide) % nalpha - - end do end associate @@ -189,12 +184,23 @@ contains integer :: i_grid ! index on nuclide energy grid integer :: i_low ! lower logarithmic mapping index integer :: i_high ! upper logarithmic mapping index + integer :: i_rxn ! reaction index + integer :: j + real(8) :: val real(8) :: f ! interp factor on nuclide energy grid real(8) :: kT ! temperature in eV real(8) :: sigT, sigA, sigF ! Intermediate multipole variables - integer :: rxn - type(TallyObject) :: t - character(15) :: threshold_rxn_type + integer, parameter :: DEPLETION_RX(6) = [N_2N, N_3N, N_4N, N_GAMMA, N_P, N_A] + + ! Initialize cached cross sections to zero + micro_xs(i_nuclide) % n2n = ZERO + micro_xs(i_nuclide) % n3n = ZERO + micro_xs(i_nuclide) % n4n = ZERO + micro_xs(i_nuclide) % np = ZERO + micro_xs(i_nuclide) % nalpha = ZERO + micro_xs(i_nuclide) % thermal = ZERO + micro_xs(i_nuclide) % thermal_elastic = ZERO + associate (nuc => nuclides(i_nuclide)) ! Check to see if there is multipole data present at this energy use_mp = .false. @@ -221,6 +227,7 @@ contains micro_xs(i_nuclide) % fission = ZERO micro_xs(i_nuclide) % nu_fission = ZERO end if + micro_xs(i_nuclide) % ngamma = sigF - sigA ! Ensure these values are set ! Note, the only time either is used is in one of 4 places: @@ -239,8 +246,6 @@ contains else ! Find the appropriate temperature index. kT = sqrtkT**2 - - !rxn = nuc % rxn_index_MT(p % event_MT) select case (temperature_method) case (TEMPERATURE_NEAREST) i_temp = minloc(abs(nuclides(i_nuclide) % kTs - kT), dim=1) @@ -286,22 +291,9 @@ contains (grid % energy(i_grid + 1) - grid % energy(i_grid)) micro_xs(i_nuclide) % index_temp = i_temp - micro_xs(i_nuclide) % index_grid = i_grid - micro_xs(i_nuclide) % interp_factor = f - ! Initialize nuclide cross-sections to zero - micro_xs(i_nuclide) % fission = ZERO - micro_xs(i_nuclide) % nu_fission = ZERO - micro_xs(i_nuclide) % n2n = ZERO - micro_xs(i_nuclide) % n3n = ZERO - micro_xs(i_nuclide) % n4n = ZERO - micro_xs(i_nuclide) % np = ZERO - micro_xs(i_nuclide) % nalpha = ZERO - micro_xs(i_nuclide) % thermal = ZERO - micro_xs(i_nuclide) % thermal_elastic = ZERO - ! Calculate microscopic nuclide total cross section micro_xs(i_nuclide) % total = (ONE - f) * xs % total(i_grid) & + f * xs % total(i_grid + 1) @@ -320,66 +312,41 @@ contains + f * xs % fission(i_grid + 1) ! Calculate microscopic nuclide nu-fission cross section - micro_xs(i_nuclide) % nu_fission = (ONE - f) * xs % nu_fission(& + micro_xs(i_nuclide) % nu_fission = (ONE - f) * xs % nu_fission( & i_grid) + f * xs % nu_fission(i_grid + 1) - + else + micro_xs(i_nuclide) % fission = ZERO + micro_xs(i_nuclide) % nu_fission = ZERO end if + end associate - micro_xs(i_nuclide) % ngamma = (micro_xs(i_nuclide) % absorption - & - micro_xs(i_nuclide) % fission) - if (t % has_threshold_rxn) then - !if (nuclides(i_nuclide)%reaction_index%has_key(p % event_MT)) the - rxn = nuclides(i_nuclide) % rxn_index_MT(p % event_MT) - associate(xs_threshold => nuc % reactions(rxn) % xs(i_temp)) - !end if - threshold_rxn_type = reaction_name(p % event_MT) - if (threshold_rxn_type == 'n2n') then - if (i_grid >= xs_threshold % threshold) then - micro_xs(i_nuclide) % n2n = (ONE - f) * xs_threshold % value(i_grid - xs_threshold & - % threshold + 1) + f * xs_threshold % value(i_grid - xs_threshold & - % threshold + 2) - end if - end if - if (threshold_rxn_type == 'n3n') then - if (i_grid >= xs_threshold % threshold) then - micro_xs(i_nuclide) % n3n = (ONE - f) * xs_threshold % value(i_grid - xs_threshold & - % threshold + 1) + f * xs_threshold % value(i_grid - xs_threshold & - % threshold + 2) - end if - end if - - - if (threshold_rxn_type == 'n4n') then - if (i_grid >= xs_threshold % threshold) then - micro_xs(i_nuclide) % n4n = (ONE - f) * xs_threshold % value(i_grid - xs_threshold & - % threshold + 1) + f * xs_threshold % value(i_grid - xs_threshold & - % threshold + 2) - end if - end if - - if (threshold_rxn_type == 'np') then - if (i_grid >= xs_threshold % threshold) then - micro_xs(i_nuclide) % np = (ONE - f) * xs_threshold % value(i_grid - xs_threshold & - % threshold + 1) + f * xs_threshold % value(i_grid - xs_threshold & - % threshold + 2) - end if - end if - - if (threshold_rxn_type == 'nalpha') then - if (i_grid >= xs_threshold % threshold) then - micro_xs(i_nuclide) % nalpha = (ONE - f) * xs_threshold % value(i_grid - xs_threshold & - % threshold + 1) + f * xs_threshold % value(i_grid - xs_threshold & - % threshold + 2) - end if - end if - - - + ! Depletion-related reactions + do j = 1, 6 + i_rxn = nuc % rxn_index_MT(DEPLETION_RX(j)) + if (i_rxn > 0) then + associate (xs => nuc % reactions(i_rxn) % xs(i_temp)) + if (i_grid >= xs % threshold) then + val = (ONE - f) * xs % value(i_grid - xs % threshold + 1) + & + f * xs % value(i_grid - xs % threshold + 2) + select case (DEPLETION_RX(j)) + case (N_2N) + micro_xs(i_nuclide) % n2n = val + case (N_3N) + micro_xs(i_nuclide) % n3n = val + case (N_4N) + micro_xs(i_nuclide) % n4n = val + case (N_GAMMA) + micro_xs(i_nuclide) % ngamma = val + case (N_P) + micro_xs(i_nuclide) % np = val + case (N_A) + micro_xs(i_nuclide) % nalpha = val + end select + end if end associate - end if - end associate - end if + end do + end if ! Initialize sab treatment to false micro_xs(i_nuclide) % index_sab = NONE diff --git a/src/endf.F90 b/src/endf.F90 index 9aace42d7b..0ba2db388b 100644 --- a/src/endf.F90 +++ b/src/endf.F90 @@ -66,18 +66,6 @@ contains string = "fission-q-prompt" case (SCORE_FISS_Q_RECOV) string = "fission-q-recoverable" - case (SCORE_N2N) - string = "n2n" - case (SCORE_N3N) - string = "n3n" - case (SCORE_N4N) - string = "n4n" - case (SCORE_NGAMMA) - string = "ngamma" - case (SCORE_NP) - string = "np" - case (SCORE_NALPHA) - string = "nalpha" ! Normal ENDF-based reactions case (TOTAL_XS) @@ -88,12 +76,12 @@ contains string = '(n,level)' case (N_2ND) string = '(n,2nd)' - !case (N_2N) - !string = '(n,2n)' - !case (N_3N) - !string = '(n,3n)' - !case (N_FISSION) - ! string = '(n,fission)' + case (N_2N) + string = '(n,2n)' + case (N_3N) + string = '(n,3n)' + case (N_FISSION) + string = '(n,fission)' case (N_F) string = '(n,f)' case (N_NF) @@ -142,18 +130,18 @@ contains string = '(n,nc)' case (N_DISAPPEAR) string = '(n,disappear)' - !case (N_GAMMA) - ! string = '(n,gamma)' - !case (N_P) - !string = '(n,p)' + case (N_GAMMA) + string = '(n,gamma)' + case (N_P) + string = '(n,p)' case (N_D) string = '(n,d)' case (N_T) string = '(n,t)' case (N_3HE) string = '(n,3He)' - !case (N_A) - !string = '(n,a)' + case (N_A) + string = '(n,a)' case (N_2A) string = '(n,2a)' case (N_3A) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 4587072f5d..91cd57ae29 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2450,17 +2450,6 @@ contains n_nuclides = index_nuclide n_sab_tables = index_sab - - do i=1, n_materials - mat => materials(i) - allocate(mat % mat_nuclide_list(n_nuclides_total)) - mat % mat_nuclide_list(:) = 0 - do j=1, mat % n_nuclides - mat % mat_nuclide_list(mat % nuclide(j))=j - end do - end do - - ! Close materials XML file call doc % clear() @@ -3058,14 +3047,14 @@ contains call fatal_error("n1n score no longer supported for tallies, & &please remove") case ('n2n', '(n,2n)') - t % score_bins(j) = SCORE_N2N - t % has_threshold_rxn = .true. + t % score_bins(j) = N_2N + case ('n3n', '(n,3n)') - t % score_bins(j) = SCORE_N3N - t % has_threshold_rxn = .true. + t % score_bins(j) = N_3N + case ('n4n', '(n,4n)') - t % score_bins(j) = SCORE_N4N - t % has_threshold_rxn = .true. + t % score_bins(j) = N_4N + case ('absorption') t % score_bins(j) = SCORE_ABSORPTION if (t % find_filter(FILTER_ENERGYOUT) > 0) then @@ -3249,12 +3238,9 @@ contains case ('(n,nc)') t % score_bins(j) = N_NC case ('(n,gamma)') - t % score_bins(j) = SCORE_NGAMMA + t % score_bins(j) = N_GAMMA case ('(n,p)') - t % score_bins(j) = SCORE_NP - - t % has_threshold_rxn = .true. - + t % score_bins(j) = N_P case ('(n,d)') t % score_bins(j) = N_D case ('(n,t)') @@ -3262,8 +3248,7 @@ contains case ('(n,3He)') t % score_bins(j) = N_3HE case ('(n,a)') - t % score_bins(j) = SCORE_NALPHA - t % has_threshold_rxn = .true. + t % score_bins(j) = N_A case ('(n,2a)') t % score_bins(j) = N_2A case ('(n,3a)') diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 1f6d45165d..14e16b268c 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -43,14 +43,6 @@ module nuclide_header real(8), allocatable :: nu_fission(:) ! neutron production real(8), allocatable :: absorption(:) ! absorption (MT > 100) real(8), allocatable :: heating(:) ! heating - real(8), allocatable :: ngamma(:) - real(8), allocatable :: n2n(:) - real(8), allocatable :: n3n(:) - real(8), allocatable :: n4n(:) - - real(8), allocatable :: np(:) - real(8), allocatable :: nalpha(:) - end type SumXS type :: Nuclide @@ -127,9 +119,9 @@ module nuclide_header real(8) :: fission real(8) :: nu_fission real(8) :: n2n - real(8) :: ngamma real(8) :: n3n real(8) :: n4n + real(8) :: ngamma real(8) :: np real(8) :: nalpha real(8) :: thermal ! Bound thermal elastic & inelastic scattering @@ -162,13 +154,12 @@ module nuclide_header real(8) :: absorption ! macroscopic absorption xs real(8) :: fission ! macroscopic fission xs real(8) :: nu_fission ! macroscopic production xs - real(8) :: n2n - real(8) :: n3n - real(8) :: n4n - real(8) :: ngamma - real(8) :: np - real(8) :: nalpha - + real(8) :: n2n ! macroscopic (n,2n) xs + real(8) :: n3n ! macroscopic (n,3n) xs + real(8) :: n4n ! macroscopic (n,4n) xs + real(8) :: ngamma ! macroscopic (n,gamma) xs + real(8) :: np ! macroscopic (n,p) xs + real(8) :: nalpha ! macroscopic (n,alpha) xs end type MaterialMacroXS !=============================================================================== diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 811ead4f8b..373c695109 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -1131,223 +1131,81 @@ contains end if end if - - case (SCORE_N2N) - - score = ZERO - - if (i_nuclide > 0) then - !if (nuclides(i_nuclide)%reaction_index%has_key(score_bin)) then - !m = nuclides(i_nuclide)%reaction_index%get_key(score_bin) - - ! Retrieve temperature and energy grid index and interpolation - ! factor - !i_temp = micro_xs(i_nuclide) % index_temp - !if (i_temp > 0) then - !i_energy = micro_xs(i_nuclide) % index_grid - !f = micro_xs(i_nuclide) % interp_factor - - !associate (xs => nuclides(i_nuclide) % reactions(m) % xs(i_temp)) - !if (i_energy >= xs % threshold) then - score = micro_xs(i_nuclide) % n2n * atom_density * flux - - - !else - ! This block is reached if multipole is turned on and we're in - ! the resolved range. For (n,gamma), use absorption - - ! fission. For everything else, assume it's zero. - - ! score = ZERO - - else - do l = 1, materials(p % material) % n_nuclides - ! Get atom density - atom_density_ = materials(p % material) % atom_density(l) - - ! Get index in nuclides array - i_nuc = materials(p % material) % nuclide(l) - - score = score + micro_xs(i_nuc) % n2n * atom_density_ * flux - end do - end if - - - case (SCORE_N3N) - - score = ZERO - - if (i_nuclide > 0) then - !if (nuclides(i_nuclide)%reaction_index%has_key(score_bin)) then - !m = nuclides(i_nuclide)%reaction_index%get_key(score_bin) - - ! Retrieve temperature and energy grid index and interpolation - ! factor - !i_temp = micro_xs(i_nuclide) % index_temp - !if (i_temp > 0) then - !i_energy = micro_xs(i_nuclide) % index_grid - !f = micro_xs(i_nuclide) % interp_factor - - !associate (xs => nuclides(i_nuclide) % reactions(m) % - !xs(i_temp)) - !if (i_energy >= xs % threshold) then - score = micro_xs(i_nuclide) % n3n * atom_density * flux - - - !else - ! This block is reached if multipole is turned on and we're in - ! the resolved range. For (n,gamma), use absorption - - ! fission. For everything else, assume it's zero. - - ! score = ZERO - - else - do l = 1, materials(p % material) % n_nuclides - ! Get atom density - atom_density_ = materials(p % material) % atom_density(l) - - ! Get index in nuclides array - i_nuc = materials(p % material) % nuclide(l) - - score = score + micro_xs(i_nuc) % n3n * atom_density_ * flux - end do - end if - - - case (SCORE_N4N) - - score = ZERO - - if (i_nuclide > 0) then - !if (nuclides(i_nuclide)%reaction_index%has_key(score_bin)) then - !m = nuclides(i_nuclide)%reaction_index%get_key(score_bin) - - ! Retrieve temperature and energy grid index and interpolation - ! factor - !i_temp = micro_xs(i_nuclide) % index_temp - !if (i_temp > 0) then - !i_energy = micro_xs(i_nuclide) % index_grid - !f = micro_xs(i_nuclide) % interp_factor - - !associate (xs => nuclides(i_nuclide) % reactions(m) % - !xs(i_temp)) - !if (i_energy >= xs % threshold) then - score = micro_xs(i_nuclide) % n4n * atom_density * flux - - - !else - ! This block is reached if multipole is turned on and we're in - ! the resolved range. For (n,gamma), use absorption - - ! fission. For everything else, assume it's zero. - - ! score = ZERO - - else - do l = 1, materials(p % material) % n_nuclides - ! Get atom density - atom_density_ = materials(p % material) % atom_density(l) - - ! Get index in nuclides array - i_nuc = materials(p % material) % nuclide(l) - - score = score + micro_xs(i_nuc) % n4n * atom_density_ * flux - end do - end if - - case (SCORE_NP) - - score = ZERO - - if (i_nuclide > 0) then - !if (nuclides(i_nuclide)%reaction_index%has_key(score_bin)) then - !m = nuclides(i_nuclide)%reaction_index%get_key(score_bin) - - ! Retrieve temperature and energy grid index and interpolation - ! factor - !i_temp = micro_xs(i_nuclide) % index_temp - !if (i_temp > 0) then - !i_energy = micro_xs(i_nuclide) % index_grid - !f = micro_xs(i_nuclide) % interp_factor - - !associate (xs => nuclides(i_nuclide) % reactions(m) % - !xs(i_temp)) - !if (i_energy >= xs % threshold) then - score = micro_xs(i_nuclide) % np * atom_density * flux - - - !else - ! This block is reached if multipole is turned on and we're in - ! the resolved range. For (n,gamma), use absorption - - ! fission. For everything else, assume it's zero. - - ! score = ZERO - - else - do l = 1, materials(p % material) % n_nuclides - ! Get atom density - atom_density_ = materials(p % material) % atom_density(l) - - ! Get index in nuclides array - i_nuc = materials(p % material) % nuclide(l) - - score = score + micro_xs(i_nuc) % np * atom_density_ * flux - end do - end if - - case (SCORE_NALPHA) - - score = ZERO - - if (i_nuclide > 0) then - !if (nuclides(i_nuclide)%reaction_index%has_key(score_bin)) then - !m = nuclides(i_nuclide)%reaction_index%get_key(score_bin) - - ! Retrieve temperature and energy grid index and interpolation - ! factor - !i_temp = micro_xs(i_nuclide) % index_temp - !if (i_temp > 0) then - !i_energy = micro_xs(i_nuclide) % index_grid - !f = micro_xs(i_nuclide) % interp_factor - - !associate (xs => nuclides(i_nuclide) % reactions(m) % - !xs(i_temp)) - !if (i_energy >= xs % threshold) then - score = micro_xs(i_nuclide) % nalpha * atom_density * flux - - - !else - ! This block is reached if multipole is turned on and we're in - ! the resolved range. For (n,gamma), use absorption - - ! fission. For everything else, assume it's zero. - - ! score = ZERO - - else - do l = 1, materials(p % material) % n_nuclides - ! Get atom density - atom_density_ = materials(p % material) % atom_density(l) - - ! Get index in nuclides array - i_nuc = materials(p % material) % nuclide(l) - - score = score + micro_xs(i_nuc) % nalpha * atom_density_ * flux - end do - end if - - - - case (SCORE_NGAMMA) + case (N_2N) if (t % estimator == ESTIMATOR_ANALOG) then - if (survival_biasing) then - ! No absorption events actually occur if survival biasing is on -- - ! just use weight absorbed in survival biasing - score = p % absorb_wgt * flux + ! Check if event MT matches + if (p % event_MT /= N_2N) cycle SCORE_LOOP + score = p % last_wgt * flux + + else + if (i_nuclide > 0) then + score = micro_xs(i_nuclide) % n2n * atom_density * flux else - ! Skip any event where the particle wasn't absorbed - if (p % event == EVENT_SCATTER) cycle SCORE_LOOP - ! All fission and absorption events will contribute here, so we - ! can just use the particle's weight entering the collision - score = p % last_wgt * flux + score = material_xs % n2n * flux end if + end if + + case (N_3N) + if (t % estimator == ESTIMATOR_ANALOG) then + ! Check if event MT matches + if (p % event_MT /= N_3N) cycle SCORE_LOOP + score = p % last_wgt * flux + + else + if (i_nuclide > 0) then + score = micro_xs(i_nuclide) % n3n * atom_density * flux + else + score = material_xs % n3n * flux + end if + end if + + case (N_4N) + if (t % estimator == ESTIMATOR_ANALOG) then + ! Check if event MT matches + if (p % event_MT /= N_4N) cycle SCORE_LOOP + score = p % last_wgt * flux + + else + if (i_nuclide > 0) then + score = micro_xs(i_nuclide) % n4n * atom_density * flux + else + score = material_xs % n4n * flux + end if + end if + + case (N_P) + if (t % estimator == ESTIMATOR_ANALOG) then + ! Check if event MT matches + if (p % event_MT /= N_P) cycle SCORE_LOOP + score = p % last_wgt * flux + + else + if (i_nuclide > 0) then + score = micro_xs(i_nuclide) % np * atom_density * flux + else + score = material_xs % np * flux + end if + end if + + case (N_A) + if (t % estimator == ESTIMATOR_ANALOG) then + ! Check if event MT matches + if (p % event_MT /= N_A) cycle SCORE_LOOP + score = p % last_wgt * flux + + else + if (i_nuclide > 0) then + score = micro_xs(i_nuclide) % nalpha * atom_density * flux + else + score = material_xs % nalpha * flux + end if + end if + + case (N_GAMMA) + if (t % estimator == ESTIMATOR_ANALOG) then + ! Check if event MT matches + if (p % event_MT /= N_GAMMA) cycle SCORE_LOOP + score = p % last_wgt * flux else if (i_nuclide > 0) then @@ -1357,8 +1215,6 @@ contains end if end if - - case default if (t % estimator == ESTIMATOR_ANALOG) then ! Any other score is assumed to be a MT number. Thus, we just need @@ -1395,14 +1251,8 @@ contains end associate else ! This block is reached if multipole is turned on and we're in - ! the resolved range. For (n,gamma), use absorption - - ! fission. For everything else, assume it's zero. - if (score_bin == N_GAMMA) then - score = (micro_xs(i_nuclide) % absorption - & - micro_xs(i_nuclide) % fission) * atom_density * flux - else - score = ZERO - end if + ! the resolved range. Assume xs is zero. + score = ZERO end if end if @@ -1434,16 +1284,8 @@ contains end associate else ! This block is reached if multipole is turned on and - ! we're in the resolved range. For (n,gamma), use - ! absorption - fission. For everything else, assume it's - ! zero. - if (score_bin == N_GAMMA) then - score = (micro_xs(i_nuc) % absorption & - - micro_xs(i_nuc) % fission) & - * atom_density_ * flux - else - score = ZERO - end if + ! we're in the resolved range. Assume xs is zero. + score = ZERO end if end if end do @@ -2949,7 +2791,7 @@ contains type(Particle), intent(in) :: p real(8), intent(in) :: distance - integer :: check_value + integer :: i integer :: i_tally integer :: i_filt @@ -3030,7 +2872,6 @@ contains if (p % material /= MATERIAL_VOID) then ! Get pointer to current material mat => materials(p % material) - check_value= mat % mat_nuclide_list(i_nuclide) ! Determine index of nuclide in Material % atom_density array j = mat % mat_nuclide_index(i_nuclide) diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index c283c3f16c..1778474e04 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -88,7 +88,6 @@ module tally_header ! reset property - allows a tally to be reset after every batch logical :: reset = .false. - logical :: has_threshold_rxn = .false. ! Number of realizations of tally random variables integer :: n_realizations = 0 From 1a5ebc7aaba19542d29912acaededcd99930f0e0 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 22 Dec 2017 15:14:23 +0700 Subject: [PATCH 09/14] Don't pre-calculate depletion reactions during inactive batches --- src/cross_section.F90 | 124 ++++++++++++++++++++++++------------------ 1 file changed, 71 insertions(+), 53 deletions(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 7066d96425..5e3d992c36 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -15,6 +15,7 @@ module cross_section use sab_header, only: SAlphaBeta, sab_tables use settings use simulation_header + use tally_header, only: active_tallies implicit none @@ -38,6 +39,9 @@ contains real(8) :: atom_density ! atom density of a nuclide real(8) :: sab_frac ! fraction of atoms affected by S(a,b) logical :: check_sab ! should we check for S(a,b) table? + logical :: in_active ! are we in active batches? + + in_active = (active_tallies % size() > 0) ! Set all material macroscopic cross sections to zero material_xs % total = ZERO @@ -107,7 +111,7 @@ contains .or. i_sab /= micro_xs(i_nuclide) % index_sab & .or. sab_frac /= micro_xs(i_nuclide) % sab_frac) then call calculate_nuclide_xs(i_nuclide, i_sab, p % E, i_grid, & - p % sqrtkT, sab_frac) + p % sqrtkT, sab_frac, in_active) end if ! ====================================================================== @@ -136,29 +140,31 @@ contains material_xs % nu_fission = material_xs % nu_fission + & atom_density * micro_xs(i_nuclide) % nu_fission - ! Add contributions to material macroscopic n2n cross section - material_xs % n2n = material_xs % n2n + & - atom_density * micro_xs(i_nuclide) % n2n + if (in_active) then + ! Add contributions to material macroscopic n2n cross section + material_xs % n2n = material_xs % n2n + & + atom_density * micro_xs(i_nuclide) % n2n - ! Add contributions to material macroscopic n3n cross section - material_xs % n3n = material_xs % n3n + & - atom_density * micro_xs(i_nuclide) % n3n + ! Add contributions to material macroscopic n3n cross section + material_xs % n3n = material_xs % n3n + & + atom_density * micro_xs(i_nuclide) % n3n - ! Add contributions to material macroscopic n4n cross section - material_xs % n4n = material_xs % n4n + & - atom_density * micro_xs(i_nuclide) % n4n + ! Add contributions to material macroscopic n4n cross section + material_xs % n4n = material_xs % n4n + & + atom_density * micro_xs(i_nuclide) % n4n - ! Add contributions to material macroscopic ngamma cross section - material_xs % ngamma = material_xs % ngamma + & - atom_density * micro_xs(i_nuclide) % ngamma + ! Add contributions to material macroscopic ngamma cross section + material_xs % ngamma = material_xs % ngamma + & + atom_density * micro_xs(i_nuclide) % ngamma - ! Add contributions to material macroscopic np cross section - material_xs % np = material_xs % np + & - atom_density * micro_xs(i_nuclide) % np + ! Add contributions to material macroscopic np cross section + material_xs % np = material_xs % np + & + atom_density * micro_xs(i_nuclide) % np - ! Add contributions to material macroscopic nalpha cross section - material_xs % nalpha = material_xs % nalpha + & - atom_density * micro_xs(i_nuclide) % nalpha + ! Add contributions to material macroscopic nalpha cross section + material_xs % nalpha = material_xs % nalpha + & + atom_density * micro_xs(i_nuclide) % nalpha + end if end do end associate @@ -170,7 +176,7 @@ contains !=============================================================================== subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_log_union, sqrtkT, & - sab_frac) + sab_frac, in_active) 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 @@ -178,6 +184,7 @@ contains ! material union energy grid real(8), intent(in) :: sqrtkT ! square root of kT, material dependent real(8), intent(in) :: sab_frac ! fraction of atoms affected by S(a,b) + logical, intent(in) :: in_active logical :: use_mp ! true if XS can be calculated with windowed multipole integer :: i_temp ! index for temperature @@ -185,19 +192,14 @@ contains integer :: i_low ! lower logarithmic mapping index integer :: i_high ! upper logarithmic mapping index integer :: i_rxn ! reaction index - integer :: j - real(8) :: val + integer :: j ! index in DEPLETION_RX + real(8) :: val ! temporary xs value real(8) :: f ! interp factor on nuclide energy grid real(8) :: kT ! temperature in eV real(8) :: sigT, sigA, sigF ! Intermediate multipole variables integer, parameter :: DEPLETION_RX(6) = [N_2N, N_3N, N_4N, N_GAMMA, N_P, N_A] ! Initialize cached cross sections to zero - micro_xs(i_nuclide) % n2n = ZERO - micro_xs(i_nuclide) % n3n = ZERO - micro_xs(i_nuclide) % n4n = ZERO - micro_xs(i_nuclide) % np = ZERO - micro_xs(i_nuclide) % nalpha = ZERO micro_xs(i_nuclide) % thermal = ZERO micro_xs(i_nuclide) % thermal_elastic = ZERO @@ -227,7 +229,15 @@ contains micro_xs(i_nuclide) % fission = ZERO micro_xs(i_nuclide) % nu_fission = ZERO end if - micro_xs(i_nuclide) % ngamma = sigF - sigA + + if (in_active) then + micro_xs(i_nuclide) % ngamma = sigA - sigF + micro_xs(i_nuclide) % n2n = ZERO + micro_xs(i_nuclide) % n3n = ZERO + micro_xs(i_nuclide) % n4n = ZERO + micro_xs(i_nuclide) % np = ZERO + micro_xs(i_nuclide) % nalpha = ZERO + end if ! Ensure these values are set ! Note, the only time either is used is in one of 4 places: @@ -321,31 +331,39 @@ contains end associate ! Depletion-related reactions - do j = 1, 6 - i_rxn = nuc % rxn_index_MT(DEPLETION_RX(j)) - if (i_rxn > 0) then - associate (xs => nuc % reactions(i_rxn) % xs(i_temp)) - if (i_grid >= xs % threshold) then - val = (ONE - f) * xs % value(i_grid - xs % threshold + 1) + & - f * xs % value(i_grid - xs % threshold + 2) - select case (DEPLETION_RX(j)) - case (N_2N) - micro_xs(i_nuclide) % n2n = val - case (N_3N) - micro_xs(i_nuclide) % n3n = val - case (N_4N) - micro_xs(i_nuclide) % n4n = val - case (N_GAMMA) - micro_xs(i_nuclide) % ngamma = val - case (N_P) - micro_xs(i_nuclide) % np = val - case (N_A) - micro_xs(i_nuclide) % nalpha = val - end select - end if - end associate - end if - end do + if (in_active) then + do j = 1, 6 + i_rxn = nuc % rxn_index_MT(DEPLETION_RX(j)) + if (i_rxn > 0) then + associate (xs => nuc % reactions(i_rxn) % xs(i_temp)) + if (i_grid >= xs % threshold) then + val = (ONE - f) * xs % value(i_grid - xs % threshold + 1) + & + f * xs % value(i_grid - xs % threshold + 2) + else + val = ZERO + end if + end associate + else + val = ZERO + end if + + select case (DEPLETION_RX(j)) + case (N_2N) + micro_xs(i_nuclide) % n2n = val + case (N_3N) + micro_xs(i_nuclide) % n3n = val + case (N_4N) + micro_xs(i_nuclide) % n4n = val + case (N_GAMMA) + micro_xs(i_nuclide) % ngamma = val + case (N_P) + micro_xs(i_nuclide) % np = val + case (N_A) + micro_xs(i_nuclide) % nalpha = val + end select + end do + end if + end if ! Initialize sab treatment to false From 12ce15978d123a588c1f2fad6be5e997a47cb0c8 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 22 Dec 2017 16:04:59 +0700 Subject: [PATCH 10/14] Fix repr bug in plots.py --- openmc/plots.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/plots.py b/openmc/plots.py index ad0d48d114..a60b2a3094 100644 --- a/openmc/plots.py +++ b/openmc/plots.py @@ -439,7 +439,7 @@ class Plot(IDManagerMixin): string += '{: <16}=\t{}\n'.format('\tWidth', self._width) string += '{: <16}=\t{}\n'.format('\tOrigin', self._origin) string += '{: <16}=\t{}\n'.format('\tPixels', self._origin) - string += '{: <16}=\t{}\n'.format('\tColor by', self._color) + string += '{: <16}=\t{}\n'.format('\tColor by', self._color_by) string += '{: <16}=\t{}\n'.format('\tBackground', self._background) string += '{: <16}=\t{}\n'.format('\tMask components', self._mask_components) From 26eedca6620843e5272f75996d1b9b4f7d0f7fab Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 22 Dec 2017 21:17:38 +0700 Subject: [PATCH 11/14] Revert to reaction_index name --- src/cross_section.F90 | 2 +- src/nuclide_header.F90 | 7 +++---- src/tallies/tally.F90 | 10 +++++----- 3 files changed, 9 insertions(+), 10 deletions(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 5e3d992c36..ebd8a342a0 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -333,7 +333,7 @@ contains ! Depletion-related reactions if (in_active) then do j = 1, 6 - i_rxn = nuc % rxn_index_MT(DEPLETION_RX(j)) + i_rxn = nuc % reaction_index(DEPLETION_RX(j)) if (i_rxn > 0) then associate (xs => nuc % reactions(i_rxn) % xs(i_temp)) if (i_grid >= xs % threshold) then diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 14e16b268c..f2a53faa58 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -90,7 +90,7 @@ module nuclide_header ! Array that maps MT values to index in reactions; used at tally-time. Note ! that ENDF-102 does not assign any MT values above 891. - integer :: rxn_index_MT(891) + integer :: reaction_index(891) ! Fission energy release class(Function1D), allocatable :: fission_q_prompt ! prompt neutrons, gammas @@ -575,7 +575,7 @@ contains n_temperature = size(this % kTs) allocate(this % sum_xs(n_temperature)) - this % rxn_index_MT(:) = 0 + this % reaction_index(:) = 0 do i = 1, n_temperature ! Allocate and initialize derived cross sections n_grid = size(this % grid(i) % energy) @@ -594,9 +594,8 @@ contains i_fission = 0 do i = 1, size(this % reactions) - call MTs % push_back(this % reactions(i) % MT) - this % rxn_index_MT(this % reactions(i) % MT) = i + this % reaction_index(this % reactions(i) % MT) = i associate (rx => this % reactions(i)) ! Skip total inelastic level scattering, gas production cross sections diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 373c695109..42aaf8e679 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -247,7 +247,7 @@ contains ! multiplicities of one. score = p % last_wgt * flux else - m = nuclides(p % event_nuclide) % rxn_index_MT(p % event_MT) + m = nuclides(p % event_nuclide) % reaction_index(p % event_MT) ! Get yield and apply to score associate (rxn => nuclides(p % event_nuclide) % reactions(m)) @@ -273,7 +273,7 @@ contains ! multiplicities of one. score = p % last_wgt * flux else - m = nuclides(p % event_nuclide) % rxn_index_MT(p % event_MT) + m = nuclides(p % event_nuclide) % reaction_index(p % event_MT) ! Get yield and apply to score associate (rxn => nuclides(p % event_nuclide) % reactions(m)) @@ -299,7 +299,7 @@ contains ! multiplicities of one. score = p % last_wgt * flux else - m = nuclides(p % event_nuclide) % rxn_index_MT(p % event_MT) + m = nuclides(p % event_nuclide) % reaction_index(p % event_MT) ! Get yield and apply to score associate (rxn => nuclides(p%event_nuclide)%reactions(m)) @@ -1232,7 +1232,7 @@ contains score = ZERO if (i_nuclide > 0) then - m = nuclides(i_nuclide) % rxn_index_MT(score_bin) + m = nuclides(i_nuclide) % reaction_index(score_bin) if (m /= 0) then ! Retrieve temperature and energy grid index and interpolation ! factor @@ -1265,7 +1265,7 @@ contains ! Get index in nuclides array i_nuc = materials(p % material) % nuclide(l) - m = nuclides(i_nuc) % rxn_index_MT(score_bin) + m = nuclides(i_nuc) % reaction_index(score_bin) if (m /= 0) then ! Retrieve temperature and energy grid index and ! interpolation factor From e137187510efee5befba8bb353b09b9b134460d1 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 23 Dec 2017 12:53:54 +0700 Subject: [PATCH 12/14] Don't cache macroscopic cross sections for depletion reactions --- src/cross_section.F90 | 32 -------------------- src/nuclide_header.F90 | 24 ++++++--------- src/tallies/tally.F90 | 66 ++++++++++++++++++++++++++++++++++++++---- 3 files changed, 69 insertions(+), 53 deletions(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index ebd8a342a0..0ac0835383 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -49,12 +49,6 @@ contains material_xs % absorption = ZERO material_xs % fission = ZERO material_xs % nu_fission = ZERO - material_xs % n2n = ZERO - material_xs % n3n = ZERO - material_xs % n4n = ZERO - material_xs % ngamma = ZERO - material_xs % np = ZERO - material_xs % nalpha = ZERO ! Exit subroutine if material is void if (p % material == MATERIAL_VOID) return @@ -139,32 +133,6 @@ contains ! Add contributions to material macroscopic nu-fission cross section material_xs % nu_fission = material_xs % nu_fission + & atom_density * micro_xs(i_nuclide) % nu_fission - - if (in_active) then - ! Add contributions to material macroscopic n2n cross section - material_xs % n2n = material_xs % n2n + & - atom_density * micro_xs(i_nuclide) % n2n - - ! Add contributions to material macroscopic n3n cross section - material_xs % n3n = material_xs % n3n + & - atom_density * micro_xs(i_nuclide) % n3n - - ! Add contributions to material macroscopic n4n cross section - material_xs % n4n = material_xs % n4n + & - atom_density * micro_xs(i_nuclide) % n4n - - ! Add contributions to material macroscopic ngamma cross section - material_xs % ngamma = material_xs % ngamma + & - atom_density * micro_xs(i_nuclide) % ngamma - - ! Add contributions to material macroscopic np cross section - material_xs % np = material_xs % np + & - atom_density * micro_xs(i_nuclide) % np - - ! Add contributions to material macroscopic nalpha cross section - material_xs % nalpha = material_xs % nalpha + & - atom_density * micro_xs(i_nuclide) % nalpha - end if end do end associate diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index f2a53faa58..661ea48504 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -115,15 +115,15 @@ module nuclide_header real(8) :: total real(8) :: elastic ! If sab_frac is not 1 or 0, then this value is ! averaged over bound and non-bound nuclei - real(8) :: absorption - real(8) :: fission - real(8) :: nu_fission - real(8) :: n2n - real(8) :: n3n - real(8) :: n4n - real(8) :: ngamma - real(8) :: np - real(8) :: nalpha + real(8) :: absorption ! absorption (disappearance) + real(8) :: fission ! fission + real(8) :: nu_fission ! neutron production from fission + real(8) :: n2n ! (n,2n) + real(8) :: n3n ! (n,3n) + real(8) :: n4n ! (n,4n) + real(8) :: ngamma ! (n,gamma) + real(8) :: np ! (n,p) + real(8) :: nalpha ! (n,alpha) real(8) :: thermal ! Bound thermal elastic & inelastic scattering real(8) :: thermal_elastic ! Bound thermal elastic scattering @@ -154,12 +154,6 @@ module nuclide_header real(8) :: absorption ! macroscopic absorption xs real(8) :: fission ! macroscopic fission xs real(8) :: nu_fission ! macroscopic production xs - real(8) :: n2n ! macroscopic (n,2n) xs - real(8) :: n3n ! macroscopic (n,3n) xs - real(8) :: n4n ! macroscopic (n,4n) xs - real(8) :: ngamma ! macroscopic (n,gamma) xs - real(8) :: np ! macroscopic (n,p) xs - real(8) :: nalpha ! macroscopic (n,alpha) xs end type MaterialMacroXS !=============================================================================== diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 42aaf8e679..d96ddb7cb7 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -1141,7 +1141,16 @@ contains if (i_nuclide > 0) then score = micro_xs(i_nuclide) % n2n * atom_density * flux else - score = material_xs % n2n * flux + score = ZERO + if (p % material /= MATERIAL_VOID) then + associate (mat => materials(p % material)) + do l = 1, materials(p % material) % n_nuclides + i_nuc = mat % nuclide(l) + atom_density_ = mat % atom_density(l) + score = score + micro_xs(i_nuc) % n2n * atom_density_ * flux + end do + end associate + end if end if end if @@ -1155,7 +1164,16 @@ contains if (i_nuclide > 0) then score = micro_xs(i_nuclide) % n3n * atom_density * flux else - score = material_xs % n3n * flux + score = ZERO + if (p % material /= MATERIAL_VOID) then + associate (mat => materials(p % material)) + do l = 1, materials(p % material) % n_nuclides + i_nuc = mat % nuclide(l) + atom_density_ = mat % atom_density(l) + score = score + micro_xs(i_nuc) % n3n * atom_density_ * flux + end do + end associate + end if end if end if @@ -1169,7 +1187,16 @@ contains if (i_nuclide > 0) then score = micro_xs(i_nuclide) % n4n * atom_density * flux else - score = material_xs % n4n * flux + score = ZERO + if (p % material /= MATERIAL_VOID) then + associate (mat => materials(p % material)) + do l = 1, materials(p % material) % n_nuclides + i_nuc = mat % nuclide(l) + atom_density_ = mat % atom_density(l) + score = score + micro_xs(i_nuc) % n4n * atom_density_ * flux + end do + end associate + end if end if end if @@ -1183,7 +1210,16 @@ contains if (i_nuclide > 0) then score = micro_xs(i_nuclide) % np * atom_density * flux else - score = material_xs % np * flux + score = ZERO + if (p % material /= MATERIAL_VOID) then + associate (mat => materials(p % material)) + do l = 1, materials(p % material) % n_nuclides + i_nuc = mat % nuclide(l) + atom_density_ = mat % atom_density(l) + score = score + micro_xs(i_nuc) % np * atom_density_ * flux + end do + end associate + end if end if end if @@ -1197,7 +1233,16 @@ contains if (i_nuclide > 0) then score = micro_xs(i_nuclide) % nalpha * atom_density * flux else - score = material_xs % nalpha * flux + score = ZERO + if (p % material /= MATERIAL_VOID) then + associate (mat => materials(p % material)) + do l = 1, materials(p % material) % n_nuclides + i_nuc = mat % nuclide(l) + atom_density_ = mat % atom_density(l) + score = score + micro_xs(i_nuc) % nalpha * atom_density_ * flux + end do + end associate + end if end if end if @@ -1211,7 +1256,16 @@ contains if (i_nuclide > 0) then score = micro_xs(i_nuclide) % ngamma * atom_density * flux else - score = material_xs % ngamma * flux + score = ZERO + if (p % material /= MATERIAL_VOID) then + associate (mat => materials(p % material)) + do l = 1, materials(p % material) % n_nuclides + i_nuc = mat % nuclide(l) + atom_density_ = mat % atom_density(l) + score = score + micro_xs(i_nuc) % ngamma * atom_density_ * flux + end do + end associate + end if end if end if From ae2d8721438c5a5641e884c92519d9585d8438f4 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 23 Dec 2017 13:11:17 +0700 Subject: [PATCH 13/14] Simplify depletion xs cache by making it an array in NuclideMicroXS --- src/constants.F90 | 3 + src/cross_section.F90 | 39 ++++-------- src/nuclide_header.F90 | 10 ++- src/tallies/tally.F90 | 139 ++++++----------------------------------- 4 files changed, 39 insertions(+), 152 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index 853a40c8d4..3968781e98 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -219,6 +219,9 @@ module constants N_D0 = 650, N_DC = 699, N_T0 = 700, N_TC = 749, N_3HE0 = 750, & N_3HEC = 799, N_A0 = 800, N_AC = 849, N_2N0 = 875, N_2NC = 891 + ! Depletion reactions + integer, parameter :: DEPLETION_RX(6) = [N_2N, N_3N, N_4N, N_GAMMA, N_P, N_A] + ! ACE table types integer, parameter :: & ACE_NEUTRON = 1, & ! continuous-energy neutron diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 0ac0835383..0bda4a4f11 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -165,7 +165,6 @@ contains real(8) :: f ! interp factor on nuclide energy grid real(8) :: kT ! temperature in eV real(8) :: sigT, sigA, sigF ! Intermediate multipole variables - integer, parameter :: DEPLETION_RX(6) = [N_2N, N_3N, N_4N, N_GAMMA, N_P, N_A] ! Initialize cached cross sections to zero micro_xs(i_nuclide) % thermal = ZERO @@ -199,12 +198,11 @@ contains end if if (in_active) then - micro_xs(i_nuclide) % ngamma = sigA - sigF - micro_xs(i_nuclide) % n2n = ZERO - micro_xs(i_nuclide) % n3n = ZERO - micro_xs(i_nuclide) % n4n = ZERO - micro_xs(i_nuclide) % np = ZERO - micro_xs(i_nuclide) % nalpha = ZERO + ! Initialize all reaction cross sections to zero + micro_xs(i_nuclide) % reaction(:) = ZERO + + ! Only non-zero reaction is (n,gamma) + micro_xs(i_nuclide) % reaction(4) = sigA - sigF end if ! Ensure these values are set @@ -301,34 +299,21 @@ contains ! Depletion-related reactions if (in_active) then do j = 1, 6 + ! Initialize reaction xs to zero + micro_xs(i_nuclide) % reaction(j) = ZERO + + ! If reaction is present and energy is greater than threshold, set + ! the reaction xs appropriately i_rxn = nuc % reaction_index(DEPLETION_RX(j)) if (i_rxn > 0) then associate (xs => nuc % reactions(i_rxn) % xs(i_temp)) if (i_grid >= xs % threshold) then - val = (ONE - f) * xs % value(i_grid - xs % threshold + 1) + & + micro_xs(i_nuclide) % reaction(j) = (ONE - f) * & + xs % value(i_grid - xs % threshold + 1) + & f * xs % value(i_grid - xs % threshold + 2) - else - val = ZERO end if end associate - else - val = ZERO end if - - select case (DEPLETION_RX(j)) - case (N_2N) - micro_xs(i_nuclide) % n2n = val - case (N_3N) - micro_xs(i_nuclide) % n3n = val - case (N_4N) - micro_xs(i_nuclide) % n4n = val - case (N_GAMMA) - micro_xs(i_nuclide) % ngamma = val - case (N_P) - micro_xs(i_nuclide) % np = val - case (N_A) - micro_xs(i_nuclide) % nalpha = val - end select end do end if diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 661ea48504..38e89b6483 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -118,15 +118,13 @@ module nuclide_header real(8) :: absorption ! absorption (disappearance) real(8) :: fission ! fission real(8) :: nu_fission ! neutron production from fission - real(8) :: n2n ! (n,2n) - real(8) :: n3n ! (n,3n) - real(8) :: n4n ! (n,4n) - real(8) :: ngamma ! (n,gamma) - real(8) :: np ! (n,p) - real(8) :: nalpha ! (n,alpha) real(8) :: thermal ! Bound thermal elastic & inelastic scattering real(8) :: thermal_elastic ! Bound thermal elastic scattering + ! Cross sections for depletion reactions (note that these are not stored in + ! macroscopic cache) + real(8) :: reaction(size(DEPLETION_RX)) + ! Indicies and factors needed to compute cross sections from the data tables integer :: index_grid ! Index on nuclide energy grid integer :: index_temp ! Temperature index for nuclide diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index d96ddb7cb7..c5219a390c 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -1131,15 +1131,31 @@ contains end if end if - case (N_2N) + case (N_2N, N_3N, N_4N, N_GAMMA, N_P, N_A) if (t % estimator == ESTIMATOR_ANALOG) then ! Check if event MT matches - if (p % event_MT /= N_2N) cycle SCORE_LOOP + if (p % event_MT /= score_bin) cycle SCORE_LOOP score = p % last_wgt * flux else + ! Determine index in NuclideMicroXS % reaction array + select case (score_bin) + case (N_2N) + m = 1 + case (N_3N) + m = 2 + case (N_4N) + m = 3 + case (N_GAMMA) + m = 4 + case (N_P) + m = 5 + case (N_A) + m = 6 + end select + if (i_nuclide > 0) then - score = micro_xs(i_nuclide) % n2n * atom_density * flux + score = micro_xs(i_nuclide) % reaction(m) * atom_density * flux else score = ZERO if (p % material /= MATERIAL_VOID) then @@ -1147,122 +1163,7 @@ contains do l = 1, materials(p % material) % n_nuclides i_nuc = mat % nuclide(l) atom_density_ = mat % atom_density(l) - score = score + micro_xs(i_nuc) % n2n * atom_density_ * flux - end do - end associate - end if - end if - end if - - case (N_3N) - if (t % estimator == ESTIMATOR_ANALOG) then - ! Check if event MT matches - if (p % event_MT /= N_3N) cycle SCORE_LOOP - score = p % last_wgt * flux - - else - if (i_nuclide > 0) then - score = micro_xs(i_nuclide) % n3n * atom_density * flux - else - score = ZERO - if (p % material /= MATERIAL_VOID) then - associate (mat => materials(p % material)) - do l = 1, materials(p % material) % n_nuclides - i_nuc = mat % nuclide(l) - atom_density_ = mat % atom_density(l) - score = score + micro_xs(i_nuc) % n3n * atom_density_ * flux - end do - end associate - end if - end if - end if - - case (N_4N) - if (t % estimator == ESTIMATOR_ANALOG) then - ! Check if event MT matches - if (p % event_MT /= N_4N) cycle SCORE_LOOP - score = p % last_wgt * flux - - else - if (i_nuclide > 0) then - score = micro_xs(i_nuclide) % n4n * atom_density * flux - else - score = ZERO - if (p % material /= MATERIAL_VOID) then - associate (mat => materials(p % material)) - do l = 1, materials(p % material) % n_nuclides - i_nuc = mat % nuclide(l) - atom_density_ = mat % atom_density(l) - score = score + micro_xs(i_nuc) % n4n * atom_density_ * flux - end do - end associate - end if - end if - end if - - case (N_P) - if (t % estimator == ESTIMATOR_ANALOG) then - ! Check if event MT matches - if (p % event_MT /= N_P) cycle SCORE_LOOP - score = p % last_wgt * flux - - else - if (i_nuclide > 0) then - score = micro_xs(i_nuclide) % np * atom_density * flux - else - score = ZERO - if (p % material /= MATERIAL_VOID) then - associate (mat => materials(p % material)) - do l = 1, materials(p % material) % n_nuclides - i_nuc = mat % nuclide(l) - atom_density_ = mat % atom_density(l) - score = score + micro_xs(i_nuc) % np * atom_density_ * flux - end do - end associate - end if - end if - end if - - case (N_A) - if (t % estimator == ESTIMATOR_ANALOG) then - ! Check if event MT matches - if (p % event_MT /= N_A) cycle SCORE_LOOP - score = p % last_wgt * flux - - else - if (i_nuclide > 0) then - score = micro_xs(i_nuclide) % nalpha * atom_density * flux - else - score = ZERO - if (p % material /= MATERIAL_VOID) then - associate (mat => materials(p % material)) - do l = 1, materials(p % material) % n_nuclides - i_nuc = mat % nuclide(l) - atom_density_ = mat % atom_density(l) - score = score + micro_xs(i_nuc) % nalpha * atom_density_ * flux - end do - end associate - end if - end if - end if - - case (N_GAMMA) - if (t % estimator == ESTIMATOR_ANALOG) then - ! Check if event MT matches - if (p % event_MT /= N_GAMMA) cycle SCORE_LOOP - score = p % last_wgt * flux - - else - if (i_nuclide > 0) then - score = micro_xs(i_nuclide) % ngamma * atom_density * flux - else - score = ZERO - if (p % material /= MATERIAL_VOID) then - associate (mat => materials(p % material)) - do l = 1, materials(p % material) % n_nuclides - i_nuc = mat % nuclide(l) - atom_density_ = mat % atom_density(l) - score = score + micro_xs(i_nuc) % ngamma * atom_density_ * flux + score = score + micro_xs(i_nuc) % reaction(m) * atom_density_ * flux end do end associate end if From a793b9a0878984c3e9b89d557d699c2d52f37c88 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 23 Dec 2017 16:02:48 +0700 Subject: [PATCH 14/14] Don't calculate depletion reaction cross sections unless tallies need them --- src/cross_section.F90 | 12 ++++-------- src/input_xml.F90 | 6 ++++++ src/simulation.F90 | 6 ++++-- src/simulation_header.F90 | 7 ++++--- src/tallies/tally.F90 | 3 +++ src/tallies/tally_header.F90 | 1 + 6 files changed, 22 insertions(+), 13 deletions(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 0bda4a4f11..2b3b32a9c0 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -39,9 +39,6 @@ contains real(8) :: atom_density ! atom density of a nuclide real(8) :: sab_frac ! fraction of atoms affected by S(a,b) logical :: check_sab ! should we check for S(a,b) table? - logical :: in_active ! are we in active batches? - - in_active = (active_tallies % size() > 0) ! Set all material macroscopic cross sections to zero material_xs % total = ZERO @@ -105,7 +102,7 @@ contains .or. i_sab /= micro_xs(i_nuclide) % index_sab & .or. sab_frac /= micro_xs(i_nuclide) % sab_frac) then call calculate_nuclide_xs(i_nuclide, i_sab, p % E, i_grid, & - p % sqrtkT, sab_frac, in_active) + p % sqrtkT, sab_frac) end if ! ====================================================================== @@ -144,7 +141,7 @@ contains !=============================================================================== subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_log_union, sqrtkT, & - sab_frac, in_active) + sab_frac) 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 @@ -152,7 +149,6 @@ contains ! material union energy grid real(8), intent(in) :: sqrtkT ! square root of kT, material dependent real(8), intent(in) :: sab_frac ! fraction of atoms affected by S(a,b) - logical, intent(in) :: in_active logical :: use_mp ! true if XS can be calculated with windowed multipole integer :: i_temp ! index for temperature @@ -197,7 +193,7 @@ contains micro_xs(i_nuclide) % nu_fission = ZERO end if - if (in_active) then + if (need_depletion_rx) then ! Initialize all reaction cross sections to zero micro_xs(i_nuclide) % reaction(:) = ZERO @@ -297,7 +293,7 @@ contains end associate ! Depletion-related reactions - if (in_active) then + if (need_depletion_rx) then do j = 1, 6 ! Initialize reaction xs to zero micro_xs(i_nuclide) % reaction(j) = ZERO diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 91cd57ae29..4dfdb178d9 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3048,12 +3048,15 @@ contains &please remove") case ('n2n', '(n,2n)') t % score_bins(j) = N_2N + t % depletion_rx = .true. case ('n3n', '(n,3n)') t % score_bins(j) = N_3N + t % depletion_rx = .true. case ('n4n', '(n,4n)') t % score_bins(j) = N_4N + t % depletion_rx = .true. case ('absorption') t % score_bins(j) = SCORE_ABSORPTION @@ -3239,8 +3242,10 @@ contains t % score_bins(j) = N_NC case ('(n,gamma)') t % score_bins(j) = N_GAMMA + t % depletion_rx = .true. case ('(n,p)') t % score_bins(j) = N_P + t % depletion_rx = .true. case ('(n,d)') t % score_bins(j) = N_D case ('(n,t)') @@ -3249,6 +3254,7 @@ contains t % score_bins(j) = N_3HE case ('(n,a)') t % score_bins(j) = N_A + t % depletion_rx = .true. case ('(n,2a)') t % score_bins(j) = N_2A case ('(n,3a)') diff --git a/src/simulation.F90 b/src/simulation.F90 index ade209e361..dde2332db5 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -451,8 +451,9 @@ contains end if end if - ! Reset current batch + ! Reset global variables current_batch = 0 + need_depletion_rx = .false. ! Set flag indicating initialization is done simulation_initialized = .true. @@ -535,7 +536,8 @@ contains if (check_overlaps) call print_overlap_check() end if - ! Reset initialization flag + ! Reset flags + need_depletion_rx = .false. simulation_initialized = .false. end subroutine openmc_simulation_finalize diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index 5041f0793c..3aaa213a90 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -20,10 +20,11 @@ module simulation_header ! ============================================================================ ! SIMULATION VARIABLES - integer :: current_batch ! current batch - integer :: current_gen ! current generation within a batch - integer :: total_gen = 0 ! total number of generations simulated + integer :: current_batch ! current batch + integer :: current_gen ! current generation within a batch + integer :: total_gen = 0 ! total number of generations simulated logical(C_BOOL), bind(C) :: simulation_initialized = .false. + logical :: need_depletion_rx ! need to calculate depletion reaction rx? ! ============================================================================ ! TALLY PRECISION TRIGGER VARIABLES diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index c5219a390c..a05cf39848 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -4381,6 +4381,9 @@ contains elseif (t % type == TALLY_SURFACE) then call active_surface_tallies % push_back(i) end if + + ! Check if tally contains depletion reactions and if so, set flag + if (t % depletion_rx) need_depletion_rx = .true. end if end associate end do diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index 1778474e04..5c0256dfbc 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -48,6 +48,7 @@ module tally_header integer :: estimator = ESTIMATOR_TRACKLENGTH ! collision, track-length real(8) :: volume ! volume of region logical :: active = .false. + logical :: depletion_rx = .false. ! has depletion reactions, e.g. (n,2n) integer, allocatable :: filter(:) ! index in filters array ! The stride attribute is used for determining the index in the results