From 1602fe672a8de40376bf09b22461139e00b083a0 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 20 Dec 2011 13:53:57 -0500 Subject: [PATCH] Moved cross section calculation routines into separate module. --- src/DEPENDENCIES | 11 ++ src/OBJECTS | 1 + src/cross_section.F90 | 408 ++++++++++++++++++++++++++++++++++++++++++ src/physics.F90 | 394 +--------------------------------------- 4 files changed, 422 insertions(+), 392 deletions(-) create mode 100644 src/cross_section.F90 diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index f385b0b940..977875c1cb 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -14,6 +14,16 @@ ace.o: string.o ace_header.o: constants.o ace_header.o: endf_header.o +cross_section.o: ace_header.o +cross_section.o: constants.o +cross_section.o: error.o +cross_section.o: fission.o +cross_section.o: global.o +cross_section.o: material_header.o +cross_section.o: particle_header.o +cross_section.o: random_lcg.o +cross_section.o: search.o + datatypes.o: datatypes_header.o doppler.o: constants.o @@ -147,6 +157,7 @@ particle_header.o: geometry_header.o physics.o: ace_header.o physics.o: constants.o +physics.o: cross_section.o physics.o: endf.o physics.o: error.o physics.o: fission.o diff --git a/src/OBJECTS b/src/OBJECTS index 6703d72439..bbf3b2a257 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -2,6 +2,7 @@ objects = \ ace.o \ ace_header.o \ bank_header.o \ +cross_section.o \ datatypes.o \ datatypes_header.o \ doppler.o \ diff --git a/src/cross_section.F90 b/src/cross_section.F90 new file mode 100644 index 0000000000..79ac18d00c --- /dev/null +++ b/src/cross_section.F90 @@ -0,0 +1,408 @@ +module cross_section + + use ace_header, only: Nuclide, SAB_Table, Reaction, UrrData + use constants + use error, only: fatal_error + use fission, only: nu_total + use global + use material_header, only: Material + use particle_header, only: Particle + use random_lcg, only: prn + use search, only: binary_search + + implicit none + +contains + +!=============================================================================== +! CALCULATE_XS determines the macroscopic cross sections for the material the +! particle is currently traveling through. +!=============================================================================== + + subroutine calculate_xs(p) + + type(Particle), pointer :: p + + integer :: i ! loop index over nuclides + integer :: index_nuclide ! index into nuclides array + integer :: index_sab ! index into sab_tables array + real(8) :: atom_density ! atom density of a nuclide + real(8) :: sab_threshold ! threshold for S(a,b) table + type(Material), pointer :: mat => null() ! current material + + ! Set all material macroscopic cross sections to zero + material_xs % total = ZERO + material_xs % elastic = ZERO + material_xs % absorption = ZERO + material_xs % fission = ZERO + material_xs % nu_fission = ZERO + + mat => materials(p % material) + + ! Find energy index on unionized grid + call find_energy_index(p) + + ! Check if there's an S(a,b) table for this material + if (mat % has_sab_table) then + sab_threshold = sab_tables(mat % sab_table) % threshold_inelastic + else + sab_threshold = ZERO + end if + + ! Add contribution from each nuclide in material + do i = 1, mat % n_nuclides + ! Determine microscopic cross sections for this nuclide + index_nuclide = mat % nuclide(i) + + ! Determine whether to use S(a,b) based on energy of particle and whether + ! the nuclide matches the S(a,b) table + if (p % E < sab_threshold .and. i == mat % sab_nuclide) then + index_sab = mat % sab_table + else + index_sab = 0 + end if + + ! Calculate microscopic cross section for this nuclide + if (p % E /= micro_xs(index_nuclide) % last_E) then + call calculate_nuclide_xs(p, index_nuclide, index_sab) + end if + + ! Copy atom density of nuclide in material + atom_density = mat % atom_density(i) + + ! Add contributions to material macroscopic total cross section + material_xs % total = material_xs % total + & + atom_density * micro_xs(index_nuclide) % total + + ! Add contributions to material macroscopic scattering cross section + material_xs % elastic = material_xs % elastic + & + atom_density * micro_xs(index_nuclide) % elastic + + ! Add contributions to material macroscopic absorption cross section + material_xs % absorption = material_xs % absorption + & + atom_density * micro_xs(index_nuclide) % absorption + + ! Add contributions to material macroscopic fission cross section + material_xs % fission = material_xs % fission + & + atom_density * micro_xs(index_nuclide) % fission + + ! Add contributions to material macroscopic nu-fission cross section + material_xs % nu_fission = material_xs % nu_fission + & + atom_density * micro_xs(index_nuclide) % nu_fission + end do + + end subroutine calculate_xs + +!=============================================================================== +! CALCULATE_NUCLIDE_XS determines microscopic cross sections for a nuclide of a +! given index in the nuclides array at the energy of the given particle +!=============================================================================== + + subroutine calculate_nuclide_xs(p, index_nuclide, index_sab) + + type(Particle), pointer :: p + integer, intent(in) :: index_nuclide ! index into nuclides array + integer, intent(in) :: index_sab ! index into sab_tables array + + integer :: IE ! index on nuclide energy grid + real(8) :: f ! interp factor on nuclide energy grid + type(Nuclide), pointer :: nuc => null() + + ! Set pointer to nuclide + nuc => nuclides(index_nuclide) + + ! TODO: If not using unionized energy grid, we need to find the index on the + ! nuclide energy grid using lethargy mapping or whatever other technique + + ! get index on nuclide energy grid + IE = nuc % grid_index(p % IE) + + ! check for rare case where two energy points are the same + if (nuc % energy(IE) == nuc % energy(IE+1)) IE = IE + 1 + + ! calculate interpolation factor + f = (p%E - nuc%energy(IE))/(nuc%energy(IE+1) - nuc%energy(IE)) + + micro_xs(index_nuclide) % index_grid = IE + micro_xs(index_nuclide) % interp_factor = f + + ! Initialize sab treatment to false + micro_xs(index_nuclide) % use_sab = .false. + micro_xs(index_nuclide) % elastic_sab = ZERO + + ! Initialize nuclide cross-sections to zero + micro_xs(index_nuclide) % fission = ZERO + micro_xs(index_nuclide) % nu_fission = ZERO + + ! Calculate microscopic nuclide total cross section + micro_xs(index_nuclide) % total = & + (ONE-f) * nuc % total(IE) + f * nuc % total(IE+1) + + ! Calculate microscopic nuclide total cross section + micro_xs(index_nuclide) % elastic = & + (ONE-f) * nuc % elastic(IE) + f * nuc % elastic(IE+1) + + ! Calculate microscopic nuclide absorption cross section + micro_xs(index_nuclide) % absorption = & + (ONE-f) * nuc % absorption(IE) + f * nuc % absorption(IE+1) + + if (nuc % fissionable) then + ! Calculate microscopic nuclide total cross section + micro_xs(index_nuclide) % fission = & + (ONE-f) * nuc % fission(IE) + f * nuc % fission(IE+1) + + ! Calculate microscopic nuclide nu-fission cross section + micro_xs(index_nuclide) % nu_fission = & + (ONE-f) * nuc % nu_fission(IE) + f * nuc % nu_fission(IE+1) + end if + + ! If there is S(a,b) data for this nuclide, we need to do a few + ! things. Since the total cross section was based on non-S(a,b) data, we + ! need to correct it by subtracting the non-S(a,b) elastic cross section and + ! then add back in the calculated S(a,b) elastic+inelastic cross section. + + 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 + + end subroutine calculate_nuclide_xs + +!=============================================================================== +! CALCULATE_SAB_XS determines the elastic and inelastic scattering +! cross-sections in the thermal energy range. These cross sections replace +! whatever data were taken from the normal Nuclide table. +!=============================================================================== + + subroutine calculate_sab_xs(p, index_nuclide, index_sab) + + type(Particle), pointer :: p + integer, intent(in) :: index_nuclide ! index into nuclides array + integer, intent(in) :: index_sab ! index into sab_tables array + + integer :: IE_sab ! index on S(a,b) energy grid + real(8) :: f_sab ! interp factor on S(a,b) energy grid + real(8) :: inelastic ! S(a,b) inelastic cross section + real(8) :: elastic ! S(a,b) elastic cross section + type(SAB_Table), pointer :: sab => null() + + ! Set flag that S(a,b) treatment should be used for scattering + micro_xs(index_nuclide) % use_sab = .true. + + ! Get pointer to S(a,b) table + sab => sab_tables(index_sab) + + ! Get index and interpolation factor for inelastic grid + if (p%E < sab % inelastic_e_in(1)) then + IE_sab = 1 + f_sab = ZERO + else + IE_sab = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, p%E) + f_sab = (p%E - sab%inelastic_e_in(IE_sab)) / & + (sab%inelastic_e_in(IE_sab+1) - sab%inelastic_e_in(IE_sab)) + end if + + ! Calculate S(a,b) inelastic scattering cross section + inelastic = (ONE-f_sab) * sab % inelastic_sigma(IE_sab) + f_sab * & + sab % inelastic_sigma(IE_sab + 1) + + ! Check for elastic data + if (p % E < sab % threshold_elastic) then + ! Determine whether elastic scattering is given in the coherent or + ! incoherent approximation. For coherent, the cross section is + ! represented as P/E whereas for incoherent, it is simply P + + if (sab % elastic_mode == SAB_ELASTIC_EXACT) then + if (p % E < sab % elastic_e_in(1)) then + ! If energy is below that of the lowest Bragg peak, the elastic + ! cross section will be zero + elastic = ZERO + else + IE_sab = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E) + elastic = sab % elastic_P(IE_sab) / p % E + end if + else + ! Determine index on elastic energy grid + if (p % E < sab % elastic_e_in(1)) then + IE_sab = 1 + else + IE_sab = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E) + end if + + ! Get interpolation factor for elastic grid + f_sab = (p%E - sab%elastic_e_in(IE_sab))/(sab%elastic_e_in(IE_sab+1) - & + sab%elastic_e_in(IE_sab)) + + ! Calculate S(a,b) elastic scattering cross section + elastic = (ONE-f_sab) * sab % elastic_P(IE_sab) + f_sab * & + sab % elastic_P(IE_sab + 1) + end if + else + ! No elastic data + elastic = ZERO + end if + + ! Correct total and elastic cross sections + micro_xs(index_nuclide) % total = micro_xs(index_nuclide) % total - & + micro_xs(index_nuclide) % elastic + inelastic + elastic + micro_xs(index_nuclide) % elastic = inelastic + elastic + + ! Store S(a,b) elastic cross section for sampling later + micro_xs(index_nuclide) % elastic_sab = elastic + + 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 +!=============================================================================== + + subroutine find_energy_index(p) + + type(Particle), pointer :: p + + integer :: IE ! index on union energy grid + real(8) :: E ! energy of particle + real(8) :: interp ! interpolation factor + + ! copy particle's energy + E = p % E + + ! if particle's 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 < e_grid(1)) then + IE = 1 + elseif (E > e_grid(n_grid)) then + IE = n_grid - 1 + else + IE = binary_search(e_grid, n_grid, E) + end if + + ! calculate the interpolation factor -- note this will be outside of [0,1) + ! for a particle outside the energy range of the union grid + interp = (E - e_grid(IE))/(e_grid(IE+1) - e_grid(IE)) + + ! set particle attributes + p % IE = IE + p % interp = interp + + end subroutine find_energy_index + +end module cross_section diff --git a/src/physics.F90 b/src/physics.F90 index e9e86a862b..86650d3141 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -1,7 +1,8 @@ module physics - use ace_header, only: Nuclide, Reaction, DistEnergy, UrrData + use ace_header, only: Nuclide, Reaction, DistEnergy use constants + use cross_section, only: calculate_xs, find_energy_index use endf, only: reaction_name, is_fission, is_scatter use error, only: fatal_error, warning use fission, only: nu_total, nu_prompt, nu_delayed @@ -135,397 +136,6 @@ contains end subroutine transport -!=============================================================================== -! CALCULATE_XS determines the macroscopic cross sections for the material the -! particle is currently traveling through. -!=============================================================================== - - subroutine calculate_xs(p) - - type(Particle), pointer :: p - - integer :: i ! loop index over nuclides - integer :: index_nuclide ! index into nuclides array - integer :: index_sab ! index into sab_tables array - real(8) :: atom_density ! atom density of a nuclide - real(8) :: sab_threshold ! threshold for S(a,b) table - type(Material), pointer :: mat => null() ! current material - - ! Set all material macroscopic cross sections to zero - material_xs % total = ZERO - material_xs % elastic = ZERO - material_xs % absorption = ZERO - material_xs % fission = ZERO - material_xs % nu_fission = ZERO - - mat => materials(p % material) - - ! Find energy index on unionized grid - call find_energy_index(p) - - ! Check if there's an S(a,b) table for this material - if (mat % has_sab_table) then - sab_threshold = sab_tables(mat % sab_table) % threshold_inelastic - else - sab_threshold = ZERO - end if - - ! Add contribution from each nuclide in material - do i = 1, mat % n_nuclides - ! Determine microscopic cross sections for this nuclide - index_nuclide = mat % nuclide(i) - - ! Determine whether to use S(a,b) based on energy of particle and whether - ! the nuclide matches the S(a,b) table - if (p % E < sab_threshold .and. i == mat % sab_nuclide) then - index_sab = mat % sab_table - else - index_sab = 0 - end if - - ! Calculate microscopic cross section for this nuclide - if (p % E /= micro_xs(index_nuclide) % last_E) then - call calculate_nuclide_xs(p, index_nuclide, index_sab) - end if - - ! Copy atom density of nuclide in material - atom_density = mat % atom_density(i) - - ! Add contributions to material macroscopic total cross section - material_xs % total = material_xs % total + & - atom_density * micro_xs(index_nuclide) % total - - ! Add contributions to material macroscopic scattering cross section - material_xs % elastic = material_xs % elastic + & - atom_density * micro_xs(index_nuclide) % elastic - - ! Add contributions to material macroscopic absorption cross section - material_xs % absorption = material_xs % absorption + & - atom_density * micro_xs(index_nuclide) % absorption - - ! Add contributions to material macroscopic fission cross section - material_xs % fission = material_xs % fission + & - atom_density * micro_xs(index_nuclide) % fission - - ! Add contributions to material macroscopic nu-fission cross section - material_xs % nu_fission = material_xs % nu_fission + & - atom_density * micro_xs(index_nuclide) % nu_fission - end do - - end subroutine calculate_xs - -!=============================================================================== -! CALCULATE_NUCLIDE_XS determines microscopic cross sections for a nuclide of a -! given index in the nuclides array at the energy of the given particle -!=============================================================================== - - subroutine calculate_nuclide_xs(p, index_nuclide, index_sab) - - type(Particle), pointer :: p - integer, intent(in) :: index_nuclide ! index into nuclides array - integer, intent(in) :: index_sab ! index into sab_tables array - - integer :: IE ! index on nuclide energy grid - real(8) :: f ! interp factor on nuclide energy grid - type(Nuclide), pointer :: nuc => null() - - ! Set pointer to nuclide - nuc => nuclides(index_nuclide) - - ! TODO: If not using unionized energy grid, we need to find the index on the - ! nuclide energy grid using lethargy mapping or whatever other technique - - ! get index on nuclide energy grid - IE = nuc % grid_index(p % IE) - - ! check for rare case where two energy points are the same - if (nuc % energy(IE) == nuc % energy(IE+1)) IE = IE + 1 - - ! calculate interpolation factor - f = (p%E - nuc%energy(IE))/(nuc%energy(IE+1) - nuc%energy(IE)) - - micro_xs(index_nuclide) % index_grid = IE - micro_xs(index_nuclide) % interp_factor = f - - ! Initialize sab treatment to false - micro_xs(index_nuclide) % use_sab = .false. - micro_xs(index_nuclide) % elastic_sab = ZERO - - ! Initialize nuclide cross-sections to zero - micro_xs(index_nuclide) % fission = ZERO - micro_xs(index_nuclide) % nu_fission = ZERO - - ! Calculate microscopic nuclide total cross section - micro_xs(index_nuclide) % total = & - (ONE-f) * nuc % total(IE) + f * nuc % total(IE+1) - - ! Calculate microscopic nuclide total cross section - micro_xs(index_nuclide) % elastic = & - (ONE-f) * nuc % elastic(IE) + f * nuc % elastic(IE+1) - - ! Calculate microscopic nuclide absorption cross section - micro_xs(index_nuclide) % absorption = & - (ONE-f) * nuc % absorption(IE) + f * nuc % absorption(IE+1) - - if (nuc % fissionable) then - ! Calculate microscopic nuclide total cross section - micro_xs(index_nuclide) % fission = & - (ONE-f) * nuc % fission(IE) + f * nuc % fission(IE+1) - - ! Calculate microscopic nuclide nu-fission cross section - micro_xs(index_nuclide) % nu_fission = & - (ONE-f) * nuc % nu_fission(IE) + f * nuc % nu_fission(IE+1) - end if - - ! If there is S(a,b) data for this nuclide, we need to do a few - ! things. Since the total cross section was based on non-S(a,b) data, we - ! need to correct it by subtracting the non-S(a,b) elastic cross section and - ! then add back in the calculated S(a,b) elastic+inelastic cross section. - - 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 - - end subroutine calculate_nuclide_xs - -!=============================================================================== -! CALCULATE_SAB_XS determines the elastic and inelastic scattering -! cross-sections in the thermal energy range. These cross sections replace -! whatever data were taken from the normal Nuclide table. -!=============================================================================== - - subroutine calculate_sab_xs(p, index_nuclide, index_sab) - - type(Particle), pointer :: p - integer, intent(in) :: index_nuclide ! index into nuclides array - integer, intent(in) :: index_sab ! index into sab_tables array - - integer :: IE_sab ! index on S(a,b) energy grid - real(8) :: f_sab ! interp factor on S(a,b) energy grid - real(8) :: inelastic ! S(a,b) inelastic cross section - real(8) :: elastic ! S(a,b) elastic cross section - type(SAB_Table), pointer :: sab => null() - - ! Set flag that S(a,b) treatment should be used for scattering - micro_xs(index_nuclide) % use_sab = .true. - - ! Get pointer to S(a,b) table - sab => sab_tables(index_sab) - - ! Get index and interpolation factor for inelastic grid - if (p%E < sab % inelastic_e_in(1)) then - IE_sab = 1 - f_sab = ZERO - else - IE_sab = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, p%E) - f_sab = (p%E - sab%inelastic_e_in(IE_sab)) / & - (sab%inelastic_e_in(IE_sab+1) - sab%inelastic_e_in(IE_sab)) - end if - - ! Calculate S(a,b) inelastic scattering cross section - inelastic = (ONE-f_sab) * sab % inelastic_sigma(IE_sab) + f_sab * & - sab % inelastic_sigma(IE_sab + 1) - - ! Check for elastic data - if (p % E < sab % threshold_elastic) then - ! Determine whether elastic scattering is given in the coherent or - ! incoherent approximation. For coherent, the cross section is - ! represented as P/E whereas for incoherent, it is simply P - - if (sab % elastic_mode == SAB_ELASTIC_EXACT) then - if (p % E < sab % elastic_e_in(1)) then - ! If energy is below that of the lowest Bragg peak, the elastic - ! cross section will be zero - elastic = ZERO - else - IE_sab = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E) - elastic = sab % elastic_P(IE_sab) / p % E - end if - else - ! Determine index on elastic energy grid - if (p % E < sab % elastic_e_in(1)) then - IE_sab = 1 - else - IE_sab = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E) - end if - - ! Get interpolation factor for elastic grid - f_sab = (p%E - sab%elastic_e_in(IE_sab))/(sab%elastic_e_in(IE_sab+1) - & - sab%elastic_e_in(IE_sab)) - - ! Calculate S(a,b) elastic scattering cross section - elastic = (ONE-f_sab) * sab % elastic_P(IE_sab) + f_sab * & - sab % elastic_P(IE_sab + 1) - end if - else - ! No elastic data - elastic = ZERO - end if - - ! Correct total and elastic cross sections - micro_xs(index_nuclide) % total = micro_xs(index_nuclide) % total - & - micro_xs(index_nuclide) % elastic + inelastic + elastic - micro_xs(index_nuclide) % elastic = inelastic + elastic - - ! Store S(a,b) elastic cross section for sampling later - micro_xs(index_nuclide) % elastic_sab = elastic - - 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 -!=============================================================================== - - subroutine find_energy_index(p) - - type(Particle), pointer :: p - - integer :: IE ! index on union energy grid - real(8) :: E ! energy of particle - real(8) :: interp ! interpolation factor - - ! copy particle's energy - E = p % E - - ! if particle's 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 < e_grid(1)) then - IE = 1 - elseif (E > e_grid(n_grid)) then - IE = n_grid - 1 - else - IE = binary_search(e_grid, n_grid, E) - end if - - ! calculate the interpolation factor -- note this will be outside of [0,1) - ! for a particle outside the energy range of the union grid - interp = (E - e_grid(IE))/(e_grid(IE+1) - e_grid(IE)) - - ! set particle attributes - p % IE = IE - p % interp = interp - - end subroutine find_energy_index - !=============================================================================== ! COLLISION samples a nuclide and reaction and then calls the appropriate ! routine for that reaction