From 989569ac5e4adb3eb8ea989768c87d5e3e1943fd Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Wed, 3 Feb 2016 13:25:51 -0500 Subject: [PATCH] Removing need for spectra (putting back in math) --- src/distribution_multivariate.F90 | 2 +- src/distribution_univariate.F90 | 2 +- src/energy_distribution.F90 | 10 +-- src/math.F90 | 106 +++++++++++++++++++++++++++- src/physics.F90 | 2 +- src/physics_mg.F90 | 2 +- src/source.F90 | 2 +- src/spectra.F90 | 112 ------------------------------ 8 files changed, 115 insertions(+), 123 deletions(-) delete mode 100644 src/spectra.F90 diff --git a/src/distribution_multivariate.F90 b/src/distribution_multivariate.F90 index e7be36db84..694fac3013 100644 --- a/src/distribution_multivariate.F90 +++ b/src/distribution_multivariate.F90 @@ -3,7 +3,7 @@ module distribution_multivariate use constants, only: ONE, TWO, PI use distribution_univariate, only: Distribution use random_lcg, only: prn - use spectra, only: rotate_angle + use math, only: rotate_angle implicit none diff --git a/src/distribution_univariate.F90 b/src/distribution_univariate.F90 index 0ec2959c3b..f3e4fdee2e 100644 --- a/src/distribution_univariate.F90 +++ b/src/distribution_univariate.F90 @@ -4,7 +4,7 @@ module distribution_univariate MAX_LINE_LEN, MAX_WORD_LEN use error, only: fatal_error use random_lcg, only: prn - use spectra, only: maxwell_spectrum, watt_spectrum + use math, only: maxwell_spectrum, watt_spectrum use string, only: to_lower use xml_interface diff --git a/src/energy_distribution.F90 b/src/energy_distribution.F90 index 7d136ab10e..2cfc1b1841 100644 --- a/src/energy_distribution.F90 +++ b/src/energy_distribution.F90 @@ -1,11 +1,11 @@ module energy_distribution - use constants, only: ZERO, ONE, TWO, PI, HISTOGRAM, LINEAR_LINEAR - use endf_header, only: Tab1 + use constants, only: ZERO, ONE, TWO, PI, HISTOGRAM, LINEAR_LINEAR + use endf_header, only: Tab1 use interpolation, only: interpolate_tab1 - use random_lcg, only: prn - use search, only: binary_search - use spectra, only: maxwell_spectrum, watt_spectrum + use math, only: maxwell_spectrum, watt_spectrum + use random_lcg, only: prn + use search, only: binary_search !=============================================================================== ! ENERGYDISTRIBUTION (abstract) defines an energy distribution that is a diff --git a/src/math.F90 b/src/math.F90 index 9ed82dabb2..aedf182358 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -1,6 +1,7 @@ module math use constants + use random_lcg, only: prn implicit none @@ -580,7 +581,8 @@ contains end function expand_harmonic !=============================================================================== -! EVALUATE_LEGENDRE +! EVALUATE_LEGENDRE Find the value of f(x) given a set of Legendre coefficients +! and the value of x !=============================================================================== pure function evaluate_legendre(data, x) result(val) real(8), intent(in) :: data(:) @@ -596,4 +598,106 @@ contains end function evaluate_legendre +!=============================================================================== +! ROTATE_ANGLE rotates direction cosines through a polar angle whose cosine is +! mu and through an azimuthal angle sampled uniformly. Note that this is done +! with direct sampling rather than rejection as is done in MCNP and SERPENT. +!=============================================================================== + + function rotate_angle(uvw0, mu, phi) result(uvw) + real(8), intent(in) :: uvw0(3) ! directional cosine + real(8), intent(in) :: mu ! cosine of angle in lab or CM + real(8), optional :: phi ! azimuthal angle + real(8) :: uvw(3) ! rotated directional cosine + + real(8) :: phi_ ! azimuthal angle + real(8) :: sinphi ! sine of azimuthal angle + real(8) :: cosphi ! cosine of azimuthal angle + real(8) :: a ! sqrt(1 - mu^2) + real(8) :: b ! sqrt(1 - w^2) + real(8) :: u0 ! original cosine in x direction + real(8) :: v0 ! original cosine in y direction + real(8) :: w0 ! original cosine in z direction + + ! Copy original directional cosines + u0 = uvw0(1) + v0 = uvw0(2) + w0 = uvw0(3) + + ! Sample azimuthal angle in [0,2pi) if none provided + if (present(phi)) then + phi_ = phi + else + phi_ = TWO * PI * prn() + end if + + ! Precompute factors to save flops + sinphi = sin(phi_) + cosphi = cos(phi_) + a = sqrt(max(ZERO, ONE - mu*mu)) + b = sqrt(max(ZERO, ONE - w0*w0)) + + ! Need to treat special case where sqrt(1 - w**2) is close to zero by + ! expanding about the v component rather than the w component + if (b > 1e-10) then + uvw(1) = mu*u0 + a*(u0*w0*cosphi - v0*sinphi)/b + uvw(2) = mu*v0 + a*(v0*w0*cosphi + u0*sinphi)/b + uvw(3) = mu*w0 - a*b*cosphi + else + b = sqrt(ONE - v0*v0) + uvw(1) = mu*u0 + a*(u0*v0*cosphi + w0*sinphi)/b + uvw(2) = mu*v0 - a*b*cosphi + uvw(3) = mu*w0 + a*(v0*w0*cosphi - u0*sinphi)/b + end if + + end function rotate_angle + +!=============================================================================== +! MAXWELL_SPECTRUM samples an energy from the Maxwell fission distribution based +! on a direct sampling scheme. The probability distribution function for a +! Maxwellian is given as p(x) = 2/(T*sqrt(pi))*sqrt(x/T)*exp(-x/T). This PDF can +! be sampled using rule C64 in the Monte Carlo Sampler LA-9721-MS. +!=============================================================================== + + function maxwell_spectrum(T) result(E_out) + + real(8), intent(in) :: T ! tabulated function of incoming E + real(8) :: E_out ! sampled energy + + real(8) :: r1, r2, r3 ! random numbers + real(8) :: c ! cosine of pi/2*r3 + + r1 = prn() + r2 = prn() + r3 = prn() + + ! determine cosine of pi/2*r + c = cos(PI/TWO*r3) + + ! determine outgoing energy + E_out = -T*(log(r1) + log(r2)*c*c) + + end function maxwell_spectrum + +!=============================================================================== +! WATT_SPECTRUM samples the outgoing energy from a Watt energy-dependent fission +! spectrum. Although fitted parameters exist for many nuclides, generally the +! continuous tabular distributions (LAW 4) should be used in lieu of the Watt +! spectrum. This direct sampling scheme is an unpublished scheme based on the +! original Watt spectrum derivation (See F. Brown's MC lectures). +!=============================================================================== + + function watt_spectrum(a, b) result(E_out) + + real(8), intent(in) :: a ! Watt parameter a + real(8), intent(in) :: b ! Watt parameter b + real(8) :: E_out ! energy of emitted neutron + + real(8) :: w ! sampled from Maxwellian + + w = maxwell_spectrum(a) + E_out = w + a*a*b/4. + (TWO*prn() - ONE)*sqrt(a*a*b*w) + + end function watt_spectrum + end module math diff --git a/src/physics.F90 b/src/physics.F90 index cb216467cd..13a311d5e3 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -9,6 +9,7 @@ module physics use global use interpolation, only: interpolate_tab1 use material_header, only: Material + use math use mesh, only: get_mesh_indices use nuclide_header use output, only: write_message @@ -19,7 +20,6 @@ module physics use search, only: binary_search use secondary_uncorrelated, only: UncorrelatedAngleEnergy use string, only: to_str - use spectra implicit none diff --git a/src/physics_mg.F90 b/src/physics_mg.F90 index a8a9084171..048f3ce54f 100644 --- a/src/physics_mg.F90 +++ b/src/physics_mg.F90 @@ -8,6 +8,7 @@ module physics_mg use macroxs_header, only: MacroXS_Base, MacroXSContainer use macroxs, only: sample_fission_energy, sample_scatter use material_header, only: Material + use math, only: rotate_angle use mesh, only: get_mesh_indices use output, only: write_message use particle_header, only: Particle @@ -16,7 +17,6 @@ module physics_mg use random_lcg, only: prn use scattdata_header use string, only: to_str - use spectra, only: rotate_angle implicit none diff --git a/src/source.F90 b/src/source.F90 index 601232192a..1fb9de5c58 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -14,7 +14,7 @@ module source use random_lcg, only: prn, set_particle_seed, prn_set_stream use search, only: binary_search use string, only: to_str - use spectra + use math use state_point, only: read_source_bank, write_source_bank #ifdef MPI diff --git a/src/spectra.F90 b/src/spectra.F90 deleted file mode 100644 index 7f96bb69e5..0000000000 --- a/src/spectra.F90 +++ /dev/null @@ -1,112 +0,0 @@ -module spectra - -use constants, only: ZERO, ONE, TWO, PI -use random_lcg, only: prn - -implicit none - -contains - -!=============================================================================== -! MAXWELL_SPECTRUM samples an energy from the Maxwell fission distribution based -! on a direct sampling scheme. The probability distribution function for a -! Maxwellian is given as p(x) = 2/(T*sqrt(pi))*sqrt(x/T)*exp(-x/T). This PDF can -! be sampled using rule C64 in the Monte Carlo Sampler LA-9721-MS. -!=============================================================================== - - function maxwell_spectrum(T) result(E_out) - - real(8), intent(in) :: T ! tabulated function of incoming E - real(8) :: E_out ! sampled energy - - real(8) :: r1, r2, r3 ! random numbers - real(8) :: c ! cosine of pi/2*r3 - - r1 = prn() - r2 = prn() - r3 = prn() - - ! determine cosine of pi/2*r - c = cos(PI/TWO*r3) - - ! determine outgoing energy - E_out = -T*(log(r1) + log(r2)*c*c) - - end function maxwell_spectrum - -!=============================================================================== -! WATT_SPECTRUM samples the outgoing energy from a Watt energy-dependent fission -! spectrum. Although fitted parameters exist for many nuclides, generally the -! continuous tabular distributions (LAW 4) should be used in lieu of the Watt -! spectrum. This direct sampling scheme is an unpublished scheme based on the -! original Watt spectrum derivation (See F. Brown's MC lectures). -!=============================================================================== - - function watt_spectrum(a, b) result(E_out) - - real(8), intent(in) :: a ! Watt parameter a - real(8), intent(in) :: b ! Watt parameter b - real(8) :: E_out ! energy of emitted neutron - - real(8) :: w ! sampled from Maxwellian - - w = maxwell_spectrum(a) - E_out = w + a*a*b/4.0_8 + (TWO*prn() - ONE)*sqrt(a*a*b*w) - - end function watt_spectrum - -!=============================================================================== -! ROTATE_ANGLE rotates direction cosines through a polar angle whose cosine is -! mu and through an azimuthal angle sampled uniformly. Note that this is done -! with direct sampling rather than rejection as is done in MCNP and SERPENT. -!=============================================================================== - - function rotate_angle(uvw0, mu, phi) result(uvw) - real(8), intent(in) :: uvw0(3) ! directional cosine - real(8), intent(in) :: mu ! cosine of angle in lab or CM - real(8), optional :: phi ! azimuthal angle - real(8) :: uvw(3) ! rotated directional cosine - - real(8) :: phi_ ! azimuthal angle - real(8) :: sinphi ! sine of azimuthal angle - real(8) :: cosphi ! cosine of azimuthal angle - real(8) :: a ! sqrt(1 - mu^2) - real(8) :: b ! sqrt(1 - w^2) - real(8) :: u0 ! original cosine in x direction - real(8) :: v0 ! original cosine in y direction - real(8) :: w0 ! original cosine in z direction - - ! Copy original directional cosines - u0 = uvw0(1) - v0 = uvw0(2) - w0 = uvw0(3) - - ! Sample azimuthal angle in [0,2pi) if none provided - if (present(phi)) then - phi_ = phi - else - phi_ = TWO * PI * prn() - end if - - ! Precompute factors to save flops - sinphi = sin(phi_) - cosphi = cos(phi_) - a = sqrt(max(ZERO, ONE - mu*mu)) - b = sqrt(max(ZERO, ONE - w0*w0)) - - ! Need to treat special case where sqrt(1 - w**2) is close to zero by - ! expanding about the v component rather than the w component - if (b > 1e-10) then - uvw(1) = mu*u0 + a*(u0*w0*cosphi - v0*sinphi)/b - uvw(2) = mu*v0 + a*(v0*w0*cosphi + u0*sinphi)/b - uvw(3) = mu*w0 - a*b*cosphi - else - b = sqrt(ONE - v0*v0) - uvw(1) = mu*u0 + a*(u0*v0*cosphi + w0*sinphi)/b - uvw(2) = mu*v0 - a*b*cosphi - uvw(3) = mu*w0 + a*(v0*w0*cosphi - u0*sinphi)/b - end if - - end function rotate_angle - -end module spectra \ No newline at end of file