crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Handle res scat case when i_E_low == i_E_up (assume cxs)

Browse source
This commit is contained in:
Paul Romano 2017-03-27 22:38:56 -05:00
parent f15f740abf
commit 05440f5f2a
1 changed files with 84 additions and 113 deletions
Download patch file Download diff file Expand all files Collapse all files

197
src/physics.F90
View file

@ -830,8 +830,6 @@ contains
integer :: i_E_rel ! index to trial relative energy
integer :: n_grid ! number of energies on 0K grid
logical :: reject ! resample if true
integer :: sampling_method ! method of target velocity sampling
awr = nuc % awr
@ -880,7 +878,7 @@ contains
wcf = xs_0K / xs_eff
wgt = wcf * wgt
case (RES_SCAT_DBRC)
case (RES_SCAT_DBRC, RES_SCAT_ARES)
E_red = sqrt((awr * E) / kT)
E_low = (((E_red - FOUR)**2) * kT) / awr
E_up = (((E_red + FOUR)**2) * kT) / awr
@ -905,125 +903,98 @@ contains
i_E_up = binary_search(nuc % energy_0K, n_grid, E_up)
end if
! interpolate xs since we're not exactly at the energy indices
xs_low = nuc % elastic_0K(i_E_low)
m = (nuc % elastic_0K(i_E_low + 1) - xs_low) &
& / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
xs_low = xs_low + m * (E_low - nuc % energy_0K(i_E_low))
xs_up = nuc % elastic_0K(i_E_up)
m = (nuc % elastic_0K(i_E_up + 1) - xs_up) &
& / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up))
xs_up = xs_up + m * (E_up - nuc % energy_0K(i_E_up))
! get max 0K xs value over range of practical relative energies
xs_max = max(xs_low, &
& maxval(nuc % elastic_0K(i_E_low + 1 : i_E_up - 1)), xs_up)
reject = .true.
! sample target velocities until one is accepted by the DBRC
do
! sample target velocity with the constant cross section (cxs) approx.
if (i_E_up == i_E_low) then
! Handle degenerate case -- if the upper/lower bounds occur for the same
! index, then using cxs is probably a good approximation
call sample_cxs_target_velocity(nuc, v_target, E, uvw, kT)
! perform Doppler broadening rejection correction (dbrc)
E_rel = dot_product((v_neut - v_target), (v_neut - v_target))
xs_0K = elastic_xs_0K(E_rel, nuc)
R_dbrc = xs_0K / xs_max
if (prn() < R_dbrc) reject = .false.
if (.not. reject) exit
end do
case (RES_SCAT_ARES)
E_red = sqrt((awr * E) / kT)
E_low = (((E_red - FOUR)**2) * kT) / awr
E_up = (((E_red + FOUR)**2) * kT) / awr
! find lower and upper energy bound indices
! lower index
n_grid = size(nuc % energy_0K)
if (E_low < nuc % energy_0K(1)) then
i_E_low = 1
elseif (E_low > nuc % energy_0K(n_grid)) then
i_E_low = n_grid - 1
else
i_E_low = binary_search(nuc % energy_0K, n_grid, E_low)
end if
if (sampling_method == RES_SCAT_DBRC) then
! interpolate xs since we're not exactly at the energy indices
xs_low = nuc % elastic_0K(i_E_low)
m = (nuc % elastic_0K(i_E_low + 1) - xs_low) &
/ (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
xs_low = xs_low + m * (E_low - nuc % energy_0K(i_E_low))
xs_up = nuc % elastic_0K(i_E_up)
m = (nuc % elastic_0K(i_E_up + 1) - xs_up) &
/ (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up))
xs_up = xs_up + m * (E_up - nuc % energy_0K(i_E_up))
! upper index
if (E_up < nuc % energy_0K(1)) then
i_E_up = 1
elseif (E_up > nuc % energy_0K(n_grid)) then
i_E_up = n_grid - 1
else
i_E_up = binary_search(nuc % energy_0K, n_grid, E_up)
end if
! get max 0K xs value over range of practical relative energies
xs_max = max(xs_low, &
maxval(nuc % elastic_0K(i_E_low + 1 : i_E_up)), xs_up)
! interpolate xs CDF since we're not exactly at the energy indices
! cdf value at lower bound attainable energy
if (i_E_low > 1) then
m = (nuc % xs_cdf(i_E_low) - nuc % xs_cdf(i_E_low - 1)) &
& / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
cdf_low = nuc % xs_cdf(i_E_low - 1) &
& + m * (E_low - nuc % energy_0K(i_E_low))
else
m = nuc % xs_cdf(i_E_low) &
& / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
cdf_low = m * (E_low - nuc % energy_0K(i_E_low))
if (E_low <= nuc % energy_0K(1)) cdf_low = ZERO
end if
DBRC_REJECT_LOOP: do
! sample target velocity with the constant cross section (cxs) approx.
call sample_cxs_target_velocity(nuc, v_target, E, uvw, kT)
! cdf value at upper bound attainable energy
m = (nuc % xs_cdf(i_E_up) - nuc % xs_cdf(i_E_up - 1)) &
& / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up))
cdf_up = nuc % xs_cdf(i_E_up - 1) &
& + m * (E_up - nuc % energy_0K(i_E_up))
! perform Doppler broadening rejection correction (dbrc)
E_rel = dot_product((v_neut - v_target), (v_neut - v_target))
xs_0K = elastic_xs_0K(E_rel, nuc)
R_dbrc = xs_0K / xs_max
if (prn() < R_dbrc) exit DBRC_REJECT_LOOP
end do DBRC_REJECT_LOOP
! values used to sample the Maxwellian
E_mode = kT
p_mode = TWO * sqrt(E_mode / pi) * sqrt((ONE / kT)**3) &
& * exp(-E_mode / kT)
E_t_max = 16.0_8 * E_mode
reject = .true.
do
! perform Maxwellian rejection sampling
E_t = E_t_max * prn()**2
p_t = TWO * sqrt(E_t / pi) * sqrt((ONE / kT)**3) &
& * exp(-E_t / kT)
R_speed = p_t / p_mode
if (prn() < R_speed) then
! sample a relative energy using the xs cdf
cdf_rel = cdf_low + prn() * (cdf_up - cdf_low)
i_E_rel = binary_search(nuc % xs_cdf(i_E_low-1:i_E_up), &
& i_E_up - i_E_low + 2, cdf_rel)
E_rel = nuc % energy_0K(i_E_low + i_E_rel - 1)
m = (nuc % xs_cdf(i_E_low + i_E_rel - 1) &
& - nuc % xs_cdf(i_E_low + i_E_rel - 2)) &
& / (nuc % energy_0K(i_E_low + i_E_rel) &
& - nuc % energy_0K(i_E_low + i_E_rel - 1))
E_rel = E_rel + (cdf_rel - nuc % xs_cdf(i_E_low + i_E_rel - 2)) / m
! perform rejection sampling on cosine between
! neutron and target velocities
mu = (E_t + awr * (E - E_rel)) / (TWO * sqrt(awr * E * E_t))
if (abs(mu) < ONE) then
! set and accept target velocity
E_t = E_t / awr
v_target = sqrt(E_t) * rotate_angle(uvw, mu)
reject = .false.
elseif (sampling_method == RES_SCAT_ARES) then
! interpolate xs CDF since we're not exactly at the energy indices
! cdf value at lower bound attainable energy
if (i_E_low > 1) then
m = (nuc % xs_cdf(i_E_low) - nuc % xs_cdf(i_E_low - 1)) &
/ (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
cdf_low = nuc % xs_cdf(i_E_low - 1) &
+ m * (E_low - nuc % energy_0K(i_E_low))
else
m = nuc % xs_cdf(i_E_low) &
/ (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
cdf_low = m * (E_low - nuc % energy_0K(i_E_low))
if (E_low <= nuc % energy_0K(1)) cdf_low = ZERO
end if
end if
if (.not. reject) exit
end do
! cdf value at upper bound attainable energy
m = (nuc % xs_cdf(i_E_up) - nuc % xs_cdf(i_E_up - 1)) &
/ (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up))
cdf_up = nuc % xs_cdf(i_E_up - 1) &
+ m * (E_up - nuc % energy_0K(i_E_up))
! values used to sample the Maxwellian
E_mode = kT
p_mode = TWO * sqrt(E_mode / pi) * sqrt((ONE / kT)**3) &
* exp(-E_mode / kT)
E_t_max = 16.0_8 * E_mode
ARES_REJECT_LOOP: do
! perform Maxwellian rejection sampling
E_t = E_t_max * prn()**2
p_t = TWO * sqrt(E_t / pi) * sqrt((ONE / kT)**3) &
* exp(-E_t / kT)
R_speed = p_t / p_mode
if (prn() < R_speed) then
! sample a relative energy using the xs cdf
cdf_rel = cdf_low + prn() * (cdf_up - cdf_low)
i_E_rel = binary_search(nuc % xs_cdf(i_E_low-1:i_E_up), &
i_E_up - i_E_low + 2, cdf_rel)
E_rel = nuc % energy_0K(i_E_low + i_E_rel - 1)
m = (nuc % xs_cdf(i_E_low + i_E_rel - 1) &
- nuc % xs_cdf(i_E_low + i_E_rel - 2)) &
/ (nuc % energy_0K(i_E_low + i_E_rel) &
- nuc % energy_0K(i_E_low + i_E_rel - 1))
E_rel = E_rel + (cdf_rel - nuc % xs_cdf(i_E_low + i_E_rel - 2)) / m
! perform rejection sampling on cosine between
! neutron and target velocities
mu = (E_t + awr * (E - E_rel)) / (TWO * sqrt(awr * E * E_t))
if (abs(mu) < ONE) then
! set and accept target velocity
E_t = E_t / awr
v_target = sqrt(E_t) * rotate_angle(uvw, mu)
exit ARES_REJECT_LOOP
end if
end if
end do ARES_REJECT_LOOP
end if
end if
end select
end subroutine sample_target_velocity