Merge pull request #744 from smharper/temperature_bug

Ensure temperature index is assigned with wmp
This commit is contained in:
Paul Romano 2016-10-30 14:24:06 -05:00 committed by GitHub
commit ed9a360bfd
2 changed files with 42 additions and 39 deletions

View file

@ -157,36 +157,7 @@ contains
if (E >= nuc % multipole % start_E/1.0e6_8 .and. &
E <= nuc % multipole % end_E/1.0e6_8) then
use_mp = .true.
else
! If using multipole data but outside the RRR, pick the nearest
! temperature. Note that there is no tolerance here, so this
! temperature could be very far off!
kT = sqrtkT**2
i_temp = minloc(abs(nuclides(i_nuclide) % kTs - kT), dim=1)
end if
else
kT = sqrtkT**2
select case (temperature_method)
case (TEMPERATURE_NEAREST)
! If using nearest temperature, do linear search on temperature
do i_temp = 1, size(nuc % kTs)
if (abs(nuc % kTs(i_temp) - kT) < K_BOLTZMANN * &
temperature_tolerance) exit
end do
case (TEMPERATURE_INTERPOLATION)
! Find temperatures that bound the actual temperature
do i_temp = 1, size(nuc % kTs) - 1
if (nuc % kTs(i_temp) <= kT .and. kT < nuc % kTs(i_temp + 1)) exit
end do
! Randomly sample between temperature i and i+1
f = (kT - nuc % kTs(i_temp)) / &
(nuc % kTs(i_temp + 1) - nuc % kTs(i_temp))
if (f > prn()) i_temp = i_temp + 1
end select
end if
! Evaluate multipole or interpolate
@ -213,22 +184,44 @@ contains
! 3. tally.F90 - score_general - For tallying on MTxxx reactions.
! 4. cross_section.F90 - calculate_urr_xs - For unresolved purposes.
! It is worth noting that none of these occur in the resolved
! resonance range, so the value here does not matter.
micro_xs(i_nuclide) % index_temp = i_temp
! resonance range, so the value here does not matter. index_temp is
! set to -1 to force a segfault in case a developer messes up and tries
! to use it with multipole.
micro_xs(i_nuclide) % index_temp = -1
micro_xs(i_nuclide) % index_grid = 0
micro_xs(i_nuclide) % interp_factor = ZERO
else
! Find the appropriate temperature index.
kT = sqrtkT**2
select case (temperature_method)
case (TEMPERATURE_NEAREST)
i_temp = minloc(abs(nuclides(i_nuclide) % kTs - kT), dim=1)
case (TEMPERATURE_INTERPOLATION)
! Find temperatures that bound the actual temperature
do i_temp = 1, size(nuc % kTs) - 1
if (nuc % kTs(i_temp) <= kT .and. kT < nuc % kTs(i_temp + 1)) exit
end do
! Randomly sample between temperature i and i+1
f = (kT - nuc % kTs(i_temp)) / &
(nuc % kTs(i_temp + 1) - nuc % kTs(i_temp))
if (f > prn()) i_temp = i_temp + 1
end select
associate (grid => nuc % grid(i_temp), xs => nuc % sum_xs(i_temp))
! Determine the energy grid index using a logarithmic mapping to reduce
! the energy range over which a binary search needs to be performed
! Determine the energy grid index using a logarithmic mapping to
! reduce the energy range over which a binary search needs to be
! performed
if (E < grid % energy(1)) then
i_grid = 1
elseif (E > grid % energy(size(grid % energy))) then
i_grid = size(grid % energy) - 1
else
! Determine bounding indices based on which equal log-spaced interval
! the energy is in
! Determine bounding indices based on which equal log-spaced
! interval the energy is in
i_low = grid % grid_index(i_log_union)
i_high = grid % grid_index(i_log_union + 1) + 1

View file

@ -97,10 +97,10 @@ contains
if (nuc % fissionable) then
if (run_mode == MODE_EIGENVALUE) then
call sample_fission(i_nuclide, i_reaction)
call sample_fission(i_nuclide, p % E, i_reaction)
call create_fission_sites(p, i_nuclide, i_reaction, fission_bank, n_bank)
elseif (run_mode == MODE_FIXEDSOURCE .and. create_fission_neutrons) then
call sample_fission(i_nuclide, i_reaction)
call sample_fission(i_nuclide, p % E, i_reaction)
call create_fission_sites(p, i_nuclide, i_reaction, &
p % secondary_bank, p % n_secondary)
end if
@ -202,8 +202,9 @@ contains
! SAMPLE_FISSION
!===============================================================================
subroutine sample_fission(i_nuclide, i_reaction)
subroutine sample_fission(i_nuclide, E, i_reaction)
integer, intent(in) :: i_nuclide ! index in nuclides array
real(8), intent(in) :: E ! incident neutron energy
integer, intent(out) :: i_reaction ! index in nuc % reactions array
integer :: i
@ -220,13 +221,22 @@ contains
! If we're in the URR, by default use the first fission reaction. We also
! default to the first reaction if we know that there are no partial fission
! reactions
if (micro_xs(i_nuclide) % use_ptable .or. &
.not. nuc % has_partial_fission) then
i_reaction = nuc % index_fission(1)
return
end if
! Check to see if we are in a windowed multipole range. WMP only supports
! the first fission reaction.
if (nuc % mp_present) then
if (E >= nuc % multipole % start_E/1.0e6_8 .and. &
E <= nuc % multipole % end_E/1.0e6_8) then
i_reaction = nuc % index_fission(1)
return
end if
end if
! Get grid index and interpolatoin factor and sample fission cdf
i_temp = micro_xs(i_nuclide) % index_temp
i_grid = micro_xs(i_nuclide) % index_grid