crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-26 05:05:30 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Spread discrete angles in continuous S(a,b) treatment

Browse source
This commit is contained in:
Paul Romano 2017-01-05 10:06:25 -06:00
parent f5c05e266c
commit 254eb6148e
1 changed files with 27 additions and 13 deletions
Download patch file Download diff file Expand all files Collapse all files

40
src/physics.F90
View file

@ -546,6 +546,7 @@ contains
real(8) :: c_j, c_j1 ! cumulative probability
real(8) :: frac ! interpolation factor on outgoing energy
real(8) :: r1 ! RNG for outgoing energy
real(8) :: mu_left, mu_right ! adjacent mu values
i_temp = micro_xs(i_nuclide) % index_temp_sab
@ -737,23 +738,36 @@ contains
E = E_1 + (E - E_i1_1) * (E_J - E_1) / (E_i1_J - E_i1_1)
end if
! Find angular distribution for closest outgoing energy bin
if (r1 - c_j < c_j1 - r1) then
j = j
else
j = j + 1
end if
! Sample outgoing cosine bin
k = 1 + int(prn() * sab % n_inelastic_mu)
! Will use mu from the randomly chosen incoming and closest outgoing
! energy bins
mu = sab % inelastic_data(l) % mu(k, j)
! Rather than use the sampled discrete mu directly, it is smeared over
! a bin of width min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the
! discrete mu value itself.
associate (mu_l => sab % inelastic_data(l) % mu)
f = (r1 - c_j)/(c_j1 - c_j)
! Determine (k-1)th mu value
if (k == 1) then
mu_left = -ONE
else
mu_left = mu_l(k-1, j) + f*(mu_l(k-1, j+1) - mu_l(k-1,j))
end if
! Determine kth mu value
mu = mu_l(k, j) + f*(mu_l(k, j+1) - mu_l(k, j))
! Determine (k+1)th mu value
if (k == sab % n_inelastic_mu) then
mu_right = ONE - mu
else
mu_right = mu_l(k+1, j) + f*(mu_l(k+1, j+1) - mu_l(k+1,j)) - mu
end if
end associate
! Smear angle
mu = mu + min(mu - mu_left, mu_right - mu)*(prn() - HALF)
else
call fatal_error("Invalid secondary energy mode on S(a,b) table " &
// trim(sab_tables(i_sab) % name))
end if ! (inelastic secondary energy treatment)
end if ! (elastic or inelastic)
end associate