From f31c7e5f0a795c8d3956dbabccde5684fa49fed4 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 18 Dec 2011 15:47:23 -0500 Subject: [PATCH] Added most logic for unresolved resonance probability tables. Does not work yet. --- src/constants.F90 | 9 +++ src/cross_section.F90 | 19 ++++++ src/cross_section_header.F90 | 1 + src/global.F90 | 3 + src/physics.F90 | 119 ++++++++++++++++++++++++++++++++++- 5 files changed, 150 insertions(+), 1 deletion(-) diff --git a/src/constants.F90 b/src/constants.F90 index 4fdaf1725a..5188035f83 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -208,6 +208,15 @@ module constants ASCII = 1, & ! ASCII cross section file BINARY = 2 ! Binary cross section file + ! Probability table parameters + integer, parameter :: & + URR_CUM_PROB = 1, & + URR_TOTAL = 2, & + URR_ELASTIC = 3, & + URR_FISSION = 4, & + URR_N_GAMMA = 5, & + URR_HEATING = 6 + ! Maximum number of partial fission reactions integer, parameter :: PARTIAL_FISSION_MAX = 4 diff --git a/src/cross_section.F90 b/src/cross_section.F90 index fde72e51b9..4332615abf 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -1107,6 +1107,25 @@ contains nuc % urr_data % multiply_smooth = .true. end if + ! if the inelastic competition flag indicates that the inelastic cross + ! section should be determined from a normal reaction cross section, we need + ! to set up a pointer to that reaction + nuc % urr_inelastic = NONE + if (nuc % urr_data % inelastic_flag > 0) then + do i = 1, nuc % n_reaction + if (nuc % reactions(i) % MT == nuc % urr_data % inelastic_flag) then + nuc % urr_inelastic = i + end if + end do + + ! Abort if no corresponding inelastic reaction was found + if (nuc % urr_inelastic == NONE) then + message = "Could not find inelastic reaction specified on " & + // "unresolved resonance probability table." + call fatal_error() + end if + end if + ! allocate incident energies and probability tables N = nuc % urr_data % n_energy M = nuc % urr_data % n_prob diff --git a/src/cross_section_header.F90 b/src/cross_section_header.F90 index 3eaf6d4462..2cc7b481c7 100644 --- a/src/cross_section_header.F90 +++ b/src/cross_section_header.F90 @@ -113,6 +113,7 @@ module cross_section_header ! Unresolved resonance data logical :: urr_present + integer :: urr_inelastic type(UrrData), pointer :: urr_data => null() ! Reactions diff --git a/src/global.F90 b/src/global.F90 index 4cfb6a6758..9c65efaf8a 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -73,6 +73,9 @@ module global integer :: n_grid ! number of points on unionized grid real(8), allocatable :: e_grid(:) ! energies on unionized grid + ! Unreoslved resonance probablity tables + logical :: urr_ptables_on = .false. + ! ============================================================================ ! TALLY-RELATED VARIABLES diff --git a/src/physics.F90 b/src/physics.F90 index 0db20f28c9..7f1487ddb6 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -1,7 +1,7 @@ module physics use constants - use cross_section_header, only: Nuclide, Reaction, DistEnergy + use cross_section_header, only: Nuclide, Reaction, DistEnergy, UrrData use endf, only: reaction_name, is_fission, is_scatter use error, only: fatal_error, warning use fission, only: nu_total, nu_prompt, nu_delayed @@ -279,6 +279,16 @@ contains if (index_sab > 0) call calculate_sab_xs(p, index_nuclide, index_sab) + ! if the particle is in the unresolved resonance range and there are + ! probability tables, we need to determine cross sections from the table + + if (urr_ptables_on .and. nuc % urr_present) then + if (p % E > nuc % urr_data % energy(1) .and. & + p % E < nuc % urr_data % energy(nuc % urr_data % n_energy)) then + call calculate_urr_xs(p, index_nuclide) + end if + end if + ! Set last evaluated energy micro_xs(index_nuclide) % last_E = p % E @@ -368,6 +378,113 @@ contains end subroutine calculate_sab_xs +!=============================================================================== +! CALCULATE_URR_XS determines cross sections in the unresolved resonance range +! from probability tables +!=============================================================================== + + subroutine calculate_urr_xs(p, index_nuclide) + + type(Particle), pointer :: p + integer, intent(in) :: index_nuclide ! index into nuclides array + + integer :: i_energy ! index for energy + integer :: i_table ! index for table + real(8) :: f ! interpolation factor + real(8) :: r ! pseudo-random number + real(8) :: elastic ! smooth elastic cross section + real(8) :: absorption ! smooth absorption cross section + real(8) :: fission ! smooth fission cross section + real(8) :: inelastic ! smooth inelastic cross section + type(UrrData), pointer :: urr => null() + type(Nuclide), pointer :: nuc => null() + type(Reaction), pointer :: rxn => null() + + ! copy cross-sections already calculated + elastic = micro_xs(index_nuclide) % elastic + absorption = micro_xs(index_nuclide) % absorption + fission = micro_xs(index_nuclide) % fission + + ! get pointer to probability table + nuc => nuclides(index_nuclide) + urr => nuc % urr_data + + ! determine energy table + i_energy = 1 + do + if (p % E > urr % energy(i_energy)) exit + i_energy = i_energy + 1 + end do + + ! determine interpolation factor on table + f = (p % E - urr % energy(i_energy)) / & + (urr % energy(i_energy + 1) - urr % energy(i_energy)) + + ! sample probability table using the cumulative distribution + r = prn() + i_table = 1 + do + if (urr % prob(i_energy, URR_CUM_PROB, i_table) > r) exit + i_table = i_table + 1 + end do + + ! determine elastic, fission, and capture cross sections from probability + ! table + if (urr % interp == LINEAR_LINEAR) then + micro_xs(index_nuclide) % elastic = & + (ONE - f) * urr % prob(i_energy, URR_ELASTIC, i_table) + & + f * urr % prob(i_energy + 1, URR_ELASTIC, i_table) + micro_xs(index_nuclide) % fission = & + (ONE - f) * urr % prob(i_energy, URR_FISSION, i_table) + & + f * urr % prob(i_energy + 1, URR_FISSION, i_table) + micro_xs(index_nuclide) % absorption = & + micro_xs(index_nuclide) % fission + & + (ONE - f) * urr % prob(i_energy, URR_N_GAMMA, i_table) + & + f * urr % prob(i_energy + 1, URR_N_GAMMA, i_table) + elseif (urr % interp == LOG_LOG) then + message = "Log-log interpolation on probability table not yet supported." + call fatal_error() + end if + + ! Determine treatment of inelastic scattering + inelastic = ZERO + if (urr % inelastic_flag > 0) then + ! Get pointer to inelastic scattering reaction + rxn => nuc % reactions(nuc % urr_inelastic) + + ! Get index on energy grid and interpolation factor + i_energy = micro_xs(index_nuclide) % index_grid + f = micro_xs(index_nuclide) % interp_factor + + ! Determine inelastic scattering cross section + if (i_energy >= rxn % IE) then + inelastic = (ONE - f) * rxn % sigma(i_energy - rxn%IE + 1) + & + f * rxn % sigma(i_energy - rxn%IE + 2) + end if + end if + + ! Multiply by smooth cross-section if needed + if (urr % multiply_smooth) then + micro_xs(index_nuclide) % elastic = elastic * & + micro_xs(index_nuclide) % elastic + micro_xs(index_nuclide) % absorption = absorption * & + micro_xs(index_nuclide) % absorption + micro_xs(index_nuclide) % fission = fission * & + micro_xs(index_nuclide) % fission + end if + + ! Determine nu-fission cross section + if (nuc % fissionable) then + micro_xs(index_nuclide) % nu_fission = nu_total(nuc, p % E) * & + micro_xs(index_nuclide) % fission + end if + + ! Calculate total cross section as sum of partials + micro_xs(index_nuclide) % total = micro_xs(index_nuclide) % elastic + & + inelastic + micro_xs(index_nuclide) % absorption + + end subroutine calculate_urr_xs + !=============================================================================== ! FIND_ENERGY_INDEX determines the index on the union energy grid and the ! interpolation factor for a particle at a certain energy