crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 13:45:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge pull request #949 from paulromano/tally-optimizations

Browse source 
Optimizations to improve performance for reaction rate tallies
This commit is contained in:
Adam Nelson 2017-12-27 06:34:21 -05:00 • committed by GitHub
commit 2d7cccaa27
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
10 changed files with 203 additions and 137 deletions
Download patch file Download diff file Expand all files Collapse all files

2
openmc/plots.py
View file

@ -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)

3
src/constants.F90
View file

@ -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

69
src/cross_section.F90
View file

@ -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
@ -154,10 +155,17 @@ 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 ! 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
! Initialize cached cross sections to 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.
@ -185,6 +193,14 @@ contains
micro_xs(i_nuclide) % nu_fission = ZERO
end if
if (need_depletion_rx) then
! 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
! Note, the only time either is used is in one of 4 places:
! 1. physics.F90 - scatter - For inelastic scatter.
@ -250,12 +266,6 @@ contains
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) % 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)
@ -276,8 +286,33 @@ contains
! Calculate microscopic nuclide nu-fission cross section
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
! Depletion-related reactions
if (need_depletion_rx) 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
micro_xs(i_nuclide) % reaction(j) = (ONE - f) * &
xs % value(i_grid - xs % threshold + 1) + &
f * xs % value(i_grid - xs % threshold + 2)
end if
end associate
end if
end do
end if
end if
! Initialize sab treatment to false
@ -573,28 +608,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

6
src/input_xml.F90
View file

@ -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)')

58
src/material_header.F90
View file

@ -36,12 +36,11 @@ module material_header
real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm
real(8) :: density_gpcc ! total density in g/cm^3
! 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
! 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(:)
! S(a,b) data
integer :: n_sab = 0 ! number of S(a,b) tables
@ -63,6 +62,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
@ -79,8 +79,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
@ -88,23 +88,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
@ -113,6 +113,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

20
src/nuclide_header.F90
View file

@ -87,8 +87,10 @@ module nuclide_header
! Reactions
type(Reaction), allocatable :: reactions(:)
type(DictIntInt) :: reaction_index ! map MT values to index in reactions
! array; used at tally-time
! 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 :: reaction_index(891)
! Fission energy release
class(Function1D), allocatable :: fission_q_prompt ! prompt neutrons, gammas
@ -113,12 +115,16 @@ 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) :: absorption ! absorption (disappearance)
real(8) :: fission ! fission
real(8) :: nu_fission ! neutron production from fission
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
@ -561,7 +567,7 @@ contains
n_temperature = size(this % kTs)
allocate(this % sum_xs(n_temperature))
this % reaction_index(:) = 0
do i = 1, n_temperature
! Allocate and initialize derived cross sections
n_grid = size(this % grid(i) % energy)
@ -581,7 +587,7 @@ contains
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 % reaction_index(this % reactions(i) % MT) = i
associate (rx => this % reactions(i))
! Skip total inelastic level scattering, gas production cross sections

25
src/simulation.F90
View file

@ -18,6 +18,7 @@ module simulation
#endif
use error, only: fatal_error, write_message
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))
@ -444,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.
@ -459,8 +467,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 +478,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 +493,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
@ -526,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

7
src/simulation_header.F90
View file

@ -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

149
src/tallies/tally.F90
View file

@ -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) % 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) % reaction_index % get(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) % reaction_index % get(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))
@ -1131,6 +1131,45 @@ contains
end if
end if
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 /= 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) % reaction(m) * 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) % reaction(m) * atom_density_ * flux
end do
end associate
end if
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
@ -1148,8 +1187,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) % reaction_index(score_bin)
if (m /= 0) then
! Retrieve temperature and energy grid index and interpolation
! factor
i_temp = micro_xs(i_nuclide) % index_temp
@ -1167,14 +1206,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
@ -1187,8 +1220,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) % reaction_index(score_bin)
if (m /= 0) then
! Retrieve temperature and energy grid index and
! interpolation factor
i_temp = micro_xs(i_nuc) % index_temp
@ -1206,16 +1239,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
@ -2377,10 +2402,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 +2444,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
! ======================================================================
@ -2812,19 +2828,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
@ -2977,19 +2985,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
@ -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

1
src/tallies/tally_header.F90
View file

@ -49,6 +49,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