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) 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 0bdadd87ee..2b3b32a9c0 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 @@ -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 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index e95762a0b3..d80f419526 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/material_header.F90 b/src/material_header.F90 index 9d93131acf..a3ce32db3d 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -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 diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 676e1b500d..38e89b6483 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -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 diff --git a/src/simulation.F90 b/src/simulation.F90 index b04e967f27..e8997fc9b8 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -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 diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index 825b88129f..11be9bf9ca 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 8b1e7ca9a2..a05cf39848 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) % 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 diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index 21b4c8a586..fa397632f7 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -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