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

Merge branch 'improved-ares' into develop

Browse source
This commit is contained in:
Paul Romano 2017-04-26 14:53:01 -05:00
commit f0214f4c12
2 changed files with 25 additions and 28 deletions
Download patch file Download diff file Expand all files Collapse all files

51
src/physics.F90
View file

@ -879,9 +879,9 @@ contains
wgt = wcf * wgt
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
E_red = sqrt(awr * E / kT)
E_low = max(ZERO, E_red - FOUR)**2 * kT / awr
E_up = (E_red + FOUR)**2 * kT / awr
! find lower and upper energy bound indices
! lower index
@ -951,33 +951,30 @@ contains
+ m * (E_up - nuc % energy_0K(i_E_up))
ARES_REJECT_LOOP: do
! perform Maxwellian rejection sampling
xi = prn()
E_t = 16.0_8 * kT * xi**2
R = FOUR * xi * exp(ONE - E_t/kT)
if (prn() < R) 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
! directly sample Maxwellian
E_t = -kT * log(prn())
! perform rejection sampling on cosine between
! neutron and target velocities
mu = (E_t + awr * (E - E_rel)) / (TWO * sqrt(awr * E * E_t))
! 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
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
! 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 do ARES_REJECT_LOOP
end if

2
tests/test_resonance_scattering/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
1.457296E+00 1.246018E-02
1.432684E+00 1.233834E-02