From 7f588e878b60b89960d124307859608a2739420a Mon Sep 17 00:00:00 2001 From: Jon Walsh Date: Mon, 3 Apr 2017 12:30:41 -0700 Subject: [PATCH 1/3] sample target energy directly! --- src/physics.F90 | 45 +++++++++++++++++++++------------------------ 1 file changed, 21 insertions(+), 24 deletions(-) diff --git a/src/physics.F90 b/src/physics.F90 index a946c161b2..df7ffd6793 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -947,33 +947,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 From 008f2ec2204e1cf2339121e26413e423a981190b Mon Sep 17 00:00:00 2001 From: Jon Walsh Date: Tue, 4 Apr 2017 07:46:39 -0700 Subject: [PATCH 2/3] disallow negative resonance scattering energy bounds --- src/physics.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/physics.F90 b/src/physics.F90 index df7ffd6793..913533da86 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -875,9 +875,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 From 7491dc8fa5bcbbc65495ca45daa82d673caca4d7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 26 Apr 2017 14:45:35 -0500 Subject: [PATCH 3/3] Update resonance scattering test result --- tests/test_resonance_scattering/results_true.dat | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_resonance_scattering/results_true.dat b/tests/test_resonance_scattering/results_true.dat index 2dfadb7dd0..153e8b0250 100644 --- a/tests/test_resonance_scattering/results_true.dat +++ b/tests/test_resonance_scattering/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.457296E+00 1.246018E-02 +1.432684E+00 1.233834E-02