From 35a689c3550c1a38abd3388606f834c0ceb3026a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 12 Nov 2015 10:52:56 -0600 Subject: [PATCH 01/11] Started refactor angle and energy distributions. --- src/angle_distribution.F90 | 55 +++++ src/energy_distribution.F90 | 407 +++++++++++++++++++++++++++++++++ src/secondary_correlated.F90 | 145 ++++++++++++ src/secondary_header.F90 | 55 +++++ src/secondary_kalbach.F90 | 163 +++++++++++++ src/secondary_uncorrelated.F90 | 29 +++ 6 files changed, 854 insertions(+) create mode 100644 src/angle_distribution.F90 create mode 100644 src/energy_distribution.F90 create mode 100644 src/secondary_correlated.F90 create mode 100644 src/secondary_header.F90 create mode 100644 src/secondary_kalbach.F90 create mode 100644 src/secondary_uncorrelated.F90 diff --git a/src/angle_distribution.F90 b/src/angle_distribution.F90 new file mode 100644 index 0000000000..448bafd8f7 --- /dev/null +++ b/src/angle_distribution.F90 @@ -0,0 +1,55 @@ +module angle_distribution + + use constants, only: ZERO, ONE + use distribution_univariate, only: DistributionContainer + use random_lcg, only: prn + use search, only: binary_search + + implicit none + private + + type, public :: AngleDistribution + real(8), allocatable :: energy(:) + type(DistributionContainer), allocatable :: distribution(:) + contains + procedure :: sample => angle_sample + end type AngleDistribution + +contains + + function angle_sample(this, E) result(mu) + class(AngleDistribution), intent(in) :: this + real(8), intent(in) :: E + real(8) :: mu + + integer :: i + integer :: n + real(8) :: r + + ! determine number of incoming energies + n = size(this%energy) + + ! find energy bin and calculate interpolation factor -- if the energy is + ! outside the range of the tabulated energies, choose the first or last bins + if (E < this%energy(1)) then + i = 1 + r = ZERO + elseif (E > this%energy(n)) then + i = n - 1 + r = ONE + else + i = binary_search(this%energy, n, E) + r = (E - this%energy(i))/(this%energy(i+1) - this%energy(i)) + end if + + ! Sample between the ith and (i+1)th bin + if (r > prn()) i = i + 1 + + ! Sample i-th distribution + mu = this%distribution(i)%obj%sample() + + ! Make sure mu is in range [-1,1] + if (abs(mu) > ONE) mu = sign(ONE, mu) + end function angle_sample + +end module angle_distribution diff --git a/src/energy_distribution.F90 b/src/energy_distribution.F90 new file mode 100644 index 0000000000..2a616e0f1f --- /dev/null +++ b/src/energy_distribution.F90 @@ -0,0 +1,407 @@ +module energy_distribution + + use constants, only: ZERO, ONE, TWO, PI, HISTOGRAM, LINEAR_LINEAR + use endf_header, only: Tab1 + use error, only: fatal_error + use interpolation, only: interpolate_tab1 + 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 +! function of the incident energy of a projectile +!=============================================================================== + + type, abstract :: EnergyDistribution + contains + procedure(iSampleEnergy), deferred :: sample + end type EnergyDistribution + + abstract interface + function iSampleEnergy(this, E_in) result(E_out) + import EnergyDistribution + class(EnergyDistribution), intent(in) :: this + real(8), intent(in) :: E_in + real(8) :: E_out + end function iSampleEnergy + end interface + +!=============================================================================== +! Derived classes +!=============================================================================== + + type Array1D + real(8), allocatable :: data(:) + end type Array1D + + type, extends(EnergyDistribution) :: TabularEquiprobable + integer :: n_region + integer, allocatable :: breakpoints(:) + integer, allocatable :: interpolation(:) + real(8), allocatable :: energy_in(:) + real(8), allocatable :: energy_out(:,:) + contains + procedure :: sample => equiprobable_sample + end type TabularEquiprobable + + type, extends(EnergyDistribution) :: LevelInelastic + real(8) :: threshold + real(8) :: mass_ratio + contains + procedure :: sample => level_inelastic_sample + end type LevelInelastic + + type CTTable + integer :: interpolation + integer :: n_discrete + real(8), allocatable :: e_out(:) + real(8), allocatable :: p(:) + real(8), allocatable :: c(:) + end type CTTable + + type, extends(EnergyDistribution) :: ContinuousTabular + integer :: n_region + integer, allocatable :: breakpoints(:) + integer, allocatable :: interpolation(:) + real(8), allocatable :: energy_in(:) + type(CTTable), allocatable :: energy_out(:) + contains + procedure :: sample => continuous_sample + end type ContinuousTabular + + type, extends(EnergyDistribution) :: MaxwellEnergy + type(Tab1) :: theta + real(8) :: u + contains + procedure :: sample => maxwellenergy_sample + end type MaxwellEnergy + + type, extends(EnergyDistribution) :: Evaporation + type(Tab1) :: theta + real(8) :: u + contains + procedure :: sample => evaporation_sample + end type Evaporation + + type, extends(EnergyDistribution) :: WattEnergy + type(Tab1) :: a + type(Tab1) :: b + real(8) :: u + contains + procedure :: sample => watt_sample + end type WattEnergy + + type, extends(EnergyDistribution) :: NBodyPhaseSpace + integer :: n_bodies + real(8) :: mass_ratio + real(8) :: A + real(8) :: Q + contains + procedure :: sample => nbody_sample + end type NBodyPhaseSpace + +contains + + function equiprobable_sample(this, E_in) result(E_out) + class(TabularEquiprobable), intent(in) :: this + real(8), intent(in) :: E_in + real(8) :: E_out + + integer :: i, k, l + integer :: n_energy_in + integer :: n_energy_out + real(8) :: r ! interpolation factor on incoming energy + real(8) :: E_i_1, E_i_K ! endpoints on outgoing grid i + real(8) :: E_i1_1, E_i1_K ! endpoints on outgoing grid i+1 + real(8) :: E_1, E_K ! endpoints interpolated between i and i+1 + real(8) :: E_l_k, E_l_k1 ! adjacent E on outgoing grid l + + ! Determine number of incoming/outgoing energies + n_energy_in = size(this%energy_in) + n_energy_out = size(this%energy_out, 1) + + ! read number of interpolation regions, incoming energies, and outgoing + ! energies + ! TODO: Move this error to input + if (this%n_region > 0) then + call fatal_error("Multiple interpolation regions not supported while & + &attempting to sample equiprobable energy bins.") + end if + + ! determine index on incoming energy grid and interpolation factor + i = binary_search(this%energy_in, size(this%energy_in), E_in) + r = (E_in - this%energy_in(i)) / & + (this%energy_in(i+1) - this%energy_in(i)) + + ! Sample outgoing energy bin + k = 1 + int(n_energy_out * prn()) + + ! Determine E_1 and E_K + E_i_1 = this%energy_out(1, i) + E_i_K = this%energy_out(n_energy_out, i) + + E_i1_1 = this%energy_out(1, i+1) + E_i1_K = this%energy_out(n_energy_out, i+1) + + E_1 = E_i_1 + r*(E_i1_1 - E_i_1) + E_K = E_i_K + r*(E_i1_K - E_i_K) + + ! Randomly select between the outgoing table for incoming energy E_i and + ! E_(i+1) + if (prn() < r) then + l = i + 1 + else + l = i + end if + + ! Determine E_l_k and E_l_k+1 + E_l_k = this%energy_out(k, l) + E_l_k1 = this%energy_out(k+1, l) + + ! Determine E' (denoted here as E_out) + E_out = E_l_k + prn()*(E_l_k1 - E_l_k) + + ! Now interpolate between incident energy bins i and i + 1 + if (l == i) then + E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1) + else + E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1) + end if + end function equiprobable_sample + + function level_inelastic_sample(this, E_in) result(E_out) + class(LevelInelastic), intent(in) :: this + real(8), intent(in) :: E_in + real(8) :: E_out + + E_out = this%mass_ratio*(E_in - this%threshold) + end function level_inelastic_sample + + function continuous_sample(this, E_in) result(E_out) + class(ContinuousTabular), intent(in) :: this + real(8), intent(in) :: E_in + real(8) :: E_out + + integer :: i, k, l + integer :: n_energy_in + integer :: n_energy_out + real(8) :: r ! interpolation factor on incoming energy + real(8) :: r1 ! random number on [0,1) + real(8) :: frac ! interpolation factor on outgoing energy + real(8) :: E_i_1, E_i_K ! endpoints on outgoing grid i + real(8) :: E_i1_1, E_i1_K ! endpoints on outgoing grid i+1 + real(8) :: E_1, E_K ! endpoints interpolated between i and i+1 + real(8) :: E_l_k, E_l_k1 ! adjacent E on outgoing grid l + real(8) :: p_l_k, p_l_k1 ! adjacent p on outgoing grid l + real(8) :: c_k, c_k1 ! cumulative probability + logical :: histogram_interp + + ! read number of interpolation regions and incoming energies + if (this%n_region == 1) then + histogram_interp = (this%interpolation(1) == 1) + else if (this%n_region > 1) then + call fatal_error("Multiple interpolation regions not supported while & + &attempting to sample continuous tabular distribution.") + else + histogram_interp = .false. + end if + + ! find energy bin and calculate interpolation factor -- if the energy is + ! outside the range of the tabulated energies, choose the first or last bins + n_energy_in = size(this%energy_in) + if (E_in < this%energy_in(1)) then + i = 1 + r = ZERO + elseif (E_in > this%energy_in(n_energy_in)) then + i = n_energy_in - 1 + r = ONE + else + i = binary_search(this%energy_in, n_energy_in, E_in) + r = (E_in - this%energy_in(i)) / & + (this%energy_in(i+1) - this%energy_in(i)) + end if + + ! Sample between the ith and (i+1)th bin + if (histogram_interp) then + l = i + else + if (r > prn()) then + l = i + 1 + else + l = i + end if + end if + + ! interpolation for energy E1 and EK + n_energy_out = size(this%energy_out(i)%e_out) + E_i_1 = this%energy_out(i)%e_out(1) + E_i_K = this%energy_out(i)%e_out(n_energy_out) + + n_energy_out = size(this%energy_out(i+1)%e_out) + E_i1_1 = this%energy_out(i+1)%e_out(1) + E_i1_K = this%energy_out(i+1)%e_out(n_energy_out) + + E_1 = E_i_1 + r*(E_i1_1 - E_i_1) + E_K = E_i_K + r*(E_i1_K - E_i_K) + + ! TODO: Write error at initizliation + if (this%energy_out(l)%n_discrete > 0) then + ! discrete lines present + call fatal_error("Discrete lines in continuous tabular distributed not & + &yet supported") + end if + + ! determine outgoing energy bin + n_energy_out = size(this%energy_out(l)%e_out) + r1 = prn() + c_k = this%energy_out(l)%c(1) + do k = 1, n_energy_out - 1 + c_k1 = this%energy_out(l)%c(k+1) + if (r1 < c_k1) exit + c_k = c_k1 + end do + + ! check to make sure k is <= NP - 1 + k = min(k, n_energy_out - 1) + + E_l_k = this%energy_out(l)%e_out(k) + p_l_k = this%energy_out(l)%p(k) + if (this%energy_out(l)%interpolation == HISTOGRAM) then + ! Histogram interpolation + if (p_l_k > ZERO) then + E_out = E_l_k + (r1 - c_k)/p_l_k + else + E_out = E_l_k + end if + + elseif (this%energy_out(l)%interpolation == LINEAR_LINEAR) then + ! Linear-linear interpolation + E_l_k1 = this%energy_out(l)%e_out(k+1) + p_l_k1 = this%energy_out(l)%p(k+1) + + frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k) + if (frac == ZERO) then + E_out = E_l_k + (r1 - c_k)/p_l_k + else + E_out = E_l_k + (sqrt(max(ZERO, p_l_k*p_l_k + & + TWO*frac*(r1 - c_k))) - p_l_k)/frac + end if + end if + + ! Now interpolate between incident energy bins i and i + 1 + if (.not. histogram_interp) then + if (l == i) then + E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1) + else + E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1) + end if + end if + end function continuous_sample + + function maxwellenergy_sample(this, E_in) result(E_out) + class(MaxwellEnergy), intent(in) :: this + real(8), intent(in) :: E_in + real(8) :: E_out + + real(8) :: theta + + ! Get temperature corresponding to incoming energy + theta = interpolate_tab1(this%theta, E_in) + + do + ! sample maxwell fission spectrum + E_out = maxwell_spectrum(theta) + + ! accept energy based on restriction energy + if (E_out <= E_in - this%u) exit + end do + end function maxwellenergy_sample + + function evaporation_sample(this, E_in) result(E_out) + class(Evaporation), intent(in) :: this + real(8), intent(in) :: E_in + real(8) :: E_out + + real(8) :: theta + real(8) :: x, y, v + + ! Get temperature corresponding to incoming energy + theta = interpolate_tab1(this%theta, E_in) + + y = (E_in - this%U)/theta + v = 1 - exp(-y) + + ! sample outgoing energy based on evaporation spectrum probability + ! density function + do + x = -log((ONE - v*prn())*(ONE - v*prn())) + if (x <= y) exit + end do + + E_out = x*theta + end function evaporation_sample + + function watt_sample(this, E_in) result(E_out) + class(WattEnergy), intent(in) :: this + real(8), intent(in) :: E_in + real(8) :: E_out + + real(8) :: a, b + + ! determine Watt parameter 'a' from tabulated function + a = interpolate_tab1(this%a, E_in) + + ! determine Watt parameter 'b' from tabulated function + b = interpolate_tab1(this%b, E_in) + + do + ! Sample energy-dependent Watt fission spectrum + E_out = watt_spectrum(a, b) + + ! accept energy based on restriction energy + if (E_out <= E_in - this%u) exit + end do + end function watt_sample + + function nbody_sample(this, E_in) result(E_out) + class(NBodyPhaseSpace), intent(in) :: this + real(8), intent(in) :: E_in + real(8) :: E_out + + real(8) :: Ap + real(8) :: E_max + real(8) :: x, y, v + real(8) :: r1, r2, r3, r4, r5, r6 + + ! determine E_max parameter + Ap = this%mass_ratio + E_max = (Ap - ONE)/Ap * (this%A/(this%A + ONE)*E_in + this%Q) + + ! x is essentially a Maxwellian distribution + x = maxwell_spectrum(ONE) + + select case (this%n_bodies) + case (3) + y = maxwell_spectrum(ONE) + case (4) + r1 = prn() + r2 = prn() + r3 = prn() + y = -log(r1*r2*r3) + case (5) + r1 = prn() + r2 = prn() + r3 = prn() + r4 = prn() + r5 = prn() + r6 = prn() + y = -log(r1*r2*r3*r4) - log(r5) * cos(PI/TWO*r6)**2 + end select + + ! now determine v and E_out + v = x/(x+y) + E_out = E_max * v + end function nbody_sample + +end module energy_distribution diff --git a/src/secondary_correlated.F90 b/src/secondary_correlated.F90 new file mode 100644 index 0000000000..c3c82b588d --- /dev/null +++ b/src/secondary_correlated.F90 @@ -0,0 +1,145 @@ +module secondary_correlated + + use constants, only: ZERO, ONE, TWO, HISTOGRAM, LINEAR_LINEAR + use distribution_univariate, only: Tabular + use error, only: fatal_error + use secondary_header, only: SecondaryDistribution + use random_lcg, only: prn + use search, only: binary_search + + type AngleEnergyTable + type(Tabular) :: energy + type(Tabular), allocatable :: angle(:) + end type AngleEnergyTable + + type, extends(SecondaryDistribution) :: CorrelatedAngleEnergy + integer :: n_region + integer, allocatable :: breakpoints(:) + integer, allocatable :: interpolation(:) + real(8), allocatable :: energy_in(:) + type(AngleEnergyTable), allocatable :: table(:) + contains + procedure :: sample => correlated_sample + end type CorrelatedAngleEnergy + +contains + + subroutine correlated_sample(this, E_in, E_out, mu) + class(CorrelatedAngleEnergy), intent(in) :: this + real(8), intent(in) :: E_in + real(8), intent(out) :: E_out + real(8), intent(out) :: mu + + integer :: i, k, l + integer :: n_energy_in + integer :: n_energy_out + real(8) :: r ! interpolation factor on incoming energy + real(8) :: r1 ! random number on [0,1) + real(8) :: frac ! interpolation factor on outgoing energy + real(8) :: E_i_1, E_i_K ! endpoints on outgoing grid i + real(8) :: E_i1_1, E_i1_K ! endpoints on outgoing grid i+1 + real(8) :: E_1, E_K ! endpoints interpolated between i and i+1 + real(8) :: E_l_k, E_l_k1 ! adjacent E on outgoing grid l + real(8) :: p_l_k, p_l_k1 ! adjacent p on outgoing grid l + real(8) :: c_k, c_k1 ! cumulative probability + + ! TODO: Write error during initialization + if (this%n_region > 1) then + call fatal_error("Multiple interpolation regions not supported while & + &attempting to sample Kalbach-Mann distribution.") + end if + + ! find energy bin and calculate interpolation factor -- if the energy is + ! outside the range of the tabulated energies, choose the first or last bins + n_energy_in = size(this%energy_in) + if (E_in < this%energy_in(1)) then + i = 1 + r = ZERO + elseif (E_in > this%energy_in(n_energy_in)) then + i = n_energy_in - 1 + r = ONE + else + i = binary_search(this%energy_in, n_energy_in, E_in) + r = (E_in - this%energy_in(i)) / & + (this%energy_in(i+1) - this%energy_in(i)) + end if + + ! Sample between the ith and (i+1)th bin + if (r > prn()) then + l = i + 1 + else + l = i + end if + + ! interpolation for energy E1 and EK + n_energy_out = size(this%table(i)%energy%x) + E_i_1 = this%table(i)%energy%x(1) + E_i_K = this%table(i)%energy%x(n_energy_out) + + n_energy_out = size(this%table(i+1)%energy%x) + E_i1_1 = this%table(i+1)%energy%x(1) + E_i1_K = this%table(i+1)%energy%x(n_energy_out) + + E_1 = E_i_1 + r*(E_i1_1 - E_i_1) + E_K = E_i_K + r*(E_i1_K - E_i_K) + +!!$ ! TODO: Write error at initizliation +!!$ if (this%table(l)%n_discrete > 0) then +!!$ ! discrete lines present +!!$ call fatal_error("Discrete lines in continuous tabular distributed not & +!!$ &yet supported") +!!$ end if + + ! determine outgoing energy bin + n_energy_out = size(this%table(l)%energy%x) + r1 = prn() + c_k = this%table(l)%energy%c(1) + do k = 1, n_energy_out - 1 + c_k1 = this%table(l)%energy%c(k+1) + if (r1 < c_k1) exit + c_k = c_k1 + end do + + ! check to make sure k is <= NP - 1 + k = min(k, n_energy_out - 1) + + E_l_k = this%table(l)%energy%x(k) + p_l_k = this%table(l)%energy%p(k) + if (this%table(l)%energy%interpolation == HISTOGRAM) then + ! Histogram interpolation + if (p_l_k > ZERO) then + E_out = E_l_k + (r1 - c_k)/p_l_k + else + E_out = E_l_k + end if + + elseif (this%table(l)%energy%interpolation == LINEAR_LINEAR) then + ! Linear-linear interpolation + E_l_k1 = this%table(l)%energy%x(k+1) + p_l_k1 = this%table(l)%energy%p(k+1) + + frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k) + if (frac == ZERO) then + E_out = E_l_k + (r1 - c_k)/p_l_k + else + E_out = E_l_k + (sqrt(max(ZERO, p_l_k*p_l_k + & + TWO*frac*(r1 - c_k))) - p_l_k)/frac + end if + end if + + ! Now interpolate between incident energy bins i and i + 1 + if (l == i) then + E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1) + else + E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1) + end if + + ! Find correlated angular distribution for closest outgoing energy bin + if (r1 - c_k < c_k1 - r1) then + mu = this%table(l)%angle(k)%sample() + else + mu = this%table(l)%angle(k + 1)%sample() + end if + end subroutine correlated_sample + +end module secondary_correlated diff --git a/src/secondary_header.F90 b/src/secondary_header.F90 new file mode 100644 index 0000000000..1f8ad3aed1 --- /dev/null +++ b/src/secondary_header.F90 @@ -0,0 +1,55 @@ +module secondary_header + + use endf_header, only: Tab1 + use interpolation, only: interpolate_tab1 + use random_lcg, only: prn + + type, abstract :: AngleEnergy + contains + procedure(iSampleAngleEnergy), deferred :: sample + end type AngleEnergy + + type :: AngleEnergyContainer + class(AngleEnergy), allocatable :: obj + end type AngleEnergyContainer + + type :: SecondaryDistribution + type(Tab1), allocatable :: applicability(:) + type(AngleEnergyContainer), allocatable :: distribution(:) + contains + procedure :: sample => secondary_sample + end type SecondaryDistribution + + abstract interface + subroutine iSampleAngleEnergy(this, E_in, E_out, mu) + import AngleEnergy + class(AngleEnergy), intent(in) :: this + real(8), intent(in) :: E_in + real(8), intent(out) :: E_out + real(8), intent(out) :: mu + end subroutine iSampleAngleEnergy + end interface + +contains + + subroutine secondary_sample(this, E_in, E_out, mu) + class(SecondaryDistribution), intent(in) :: this + real(8), intent(in) :: E_in + real(8), intent(out) :: E_out + real(8), intent(out) :: mu + + real(8) :: p_valid + + do i = 1, size(this%applicability) + ! Determine probability that i-th energy distribution is sampled + p_valid = interpolate_tab1(this%applicability(i), E_in) + + ! If i-th distribution is sampled, sample energy from the distribution + if (prn() <= p_valid) then + call this%distribution(i)%obj%sample(E_in, E_out, mu) + exit + end if + end do + end subroutine secondary_sample + +end module secondary_header diff --git a/src/secondary_kalbach.F90 b/src/secondary_kalbach.F90 new file mode 100644 index 0000000000..d0d2f32bf2 --- /dev/null +++ b/src/secondary_kalbach.F90 @@ -0,0 +1,163 @@ +module secondary_kalbach + + use constants, only: ZERO, ONE, TWO, HISTOGRAM, LINEAR_LINEAR + use error, only: fatal_error + use secondary_header, only: SecondaryDistribution + use random_lcg, only: prn + use search, only: binary_search + + type KalbachMannTable + integer :: n_discrete + integer :: interpolation + real(8), allocatable :: e_out(:) + real(8), allocatable :: p(:) + real(8), allocatable :: c(:) + real(8), allocatable :: r(:) + real(8), allocatable :: a(:) + end type KalbachMannTable + + type, extends(SecondaryDistribution) :: KalbachMann + integer :: n_region + integer, allocatable :: breakpoints(:) + integer, allocatable :: interpolation(:) + real(8), allocatable :: energy_in(:) + type(KalbachMannTable), allocatable :: table(:) + contains + procedure :: sample => kalbachmann_sample + end type KalbachMann + +contains + + subroutine kalbachmann_sample(this, E_in, E_out, mu) + class(KalbachMann), intent(in) :: this + real(8), intent(in) :: E_in + real(8), intent(out) :: E_out + real(8), intent(out) :: mu + + integer :: i, k, l + integer :: n_energy_in + integer :: n_energy_out + real(8) :: r ! interpolation factor on incoming energy + real(8) :: r1 ! random number on [0,1) + real(8) :: frac ! interpolation factor on outgoing energy + real(8) :: E_i_1, E_i_K ! endpoints on outgoing grid i + real(8) :: E_i1_1, E_i1_K ! endpoints on outgoing grid i+1 + real(8) :: E_1, E_K ! endpoints interpolated between i and i+1 + real(8) :: E_l_k, E_l_k1 ! adjacent E on outgoing grid l + real(8) :: p_l_k, p_l_k1 ! adjacent p on outgoing grid l + real(8) :: c_k, c_k1 ! cumulative probability + real(8) :: km_r, km_a ! Kalbach-Mann parameters + real(8) :: T + + ! TODO: Write error during initialization + if (this%n_region > 1) then + call fatal_error("Multiple interpolation regions not supported while & + &attempting to sample Kalbach-Mann distribution.") + end if + + ! find energy bin and calculate interpolation factor -- if the energy is + ! outside the range of the tabulated energies, choose the first or last bins + n_energy_in = size(this%energy_in) + if (E_in < this%energy_in(1)) then + i = 1 + r = ZERO + elseif (E_in > this%energy_in(n_energy_in)) then + i = n_energy_in - 1 + r = ONE + else + i = binary_search(this%energy_in, n_energy_in, E_in) + r = (E_in - this%energy_in(i)) / & + (this%energy_in(i+1) - this%energy_in(i)) + end if + + ! Sample between the ith and (i+1)th bin + if (r > prn()) then + l = i + 1 + else + l = i + end if + + ! interpolation for energy E1 and EK + n_energy_out = size(this%table(i)%e_out) + E_i_1 = this%table(i)%e_out(1) + E_i_K = this%table(i)%e_out(n_energy_out) + + n_energy_out = size(this%table(i+1)%e_out) + E_i1_1 = this%table(i+1)%e_out(1) + E_i1_K = this%table(i+1)%e_out(n_energy_out) + + E_1 = E_i_1 + r*(E_i1_1 - E_i_1) + E_K = E_i_K + r*(E_i1_K - E_i_K) + + ! TODO: Write error at initizliation + if (this%table(l)%n_discrete > 0) then + ! discrete lines present + call fatal_error("Discrete lines in continuous tabular distributed not & + &yet supported") + end if + + ! determine outgoing energy bin + n_energy_out = size(this%table(l)%e_out) + r1 = prn() + c_k = this%table(l)%c(1) + do k = 1, n_energy_out - 1 + c_k1 = this%table(l)%c(k+1) + if (r1 < c_k1) exit + c_k = c_k1 + end do + + ! check to make sure k is <= NP - 1 + k = min(k, n_energy_out - 1) + + E_l_k = this%table(l)%e_out(k) + p_l_k = this%table(l)%p(k) + if (this%table(l)%interpolation == HISTOGRAM) then + ! Histogram interpolation + if (p_l_k > ZERO) then + E_out = E_l_k + (r1 - c_k)/p_l_k + else + E_out = E_l_k + end if + + ! Determine Kalbach-Mann parameters + km_r = this%table(l)%r(k) + km_a = this%table(l)%a(k) + + elseif (this%table(l)%interpolation == LINEAR_LINEAR) then + ! Linear-linear interpolation + E_l_k1 = this%table(l)%e_out(k+1) + p_l_k1 = this%table(l)%p(k+1) + + frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k) + if (frac == ZERO) then + E_out = E_l_k + (r1 - c_k)/p_l_k + else + E_out = E_l_k + (sqrt(max(ZERO, p_l_k*p_l_k + & + TWO*frac*(r1 - c_k))) - p_l_k)/frac + end if + + ! Determine Kalbach-Mann parameters + km_r = this%table(l)%r(k) + (E_out - E_l_k)/(E_l_k1 - E_l_k) * & + (this%table(l)%r(k+1) - this%table(l)%r(k)) + km_a = this%table(l)%a(k) + (E_out - E_l_k)/(E_l_k1 - E_l_k) * & + (this%table(l)%a(k+1) - this%table(l)%a(k)) + end if + + ! Now interpolate between incident energy bins i and i + 1 + if (l == i) then + E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1) + else + E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1) + end if + + ! Sampled correlated angle from Kalbach-Mann parameters + if (prn() > km_r) then + T = (TWO*prn() - ONE) * sinh(km_a) + mu = log(T + sqrt(T*T + ONE))/km_a + else + r1 = prn() + mu = log(r1*exp(km_a) + (ONE - r1)*exp(-km_a))/km_a + end if + end subroutine kalbachmann_sample + +end module secondary_kalbach diff --git a/src/secondary_uncorrelated.F90 b/src/secondary_uncorrelated.F90 new file mode 100644 index 0000000000..7d0046aa4f --- /dev/null +++ b/src/secondary_uncorrelated.F90 @@ -0,0 +1,29 @@ +module secondary_uncorrelated + + use angle_distribution, only: AngleDistribution + use energy_distribution, only: EnergyDistribution + use secondary_header, only: SecondaryDistribution + + type, extends(SecondaryDistribution) :: UncorrelatedAngleEnergy + type(AngleDistribution) :: angle + class(EnergyDistribution), allocatable :: energy + contains + procedure :: sample => uncorrelated_sample + end type UncorrelatedAngleEnergy + +contains + + subroutine uncorrelated_sample(this, E_in, E_out, mu) + class(UncorrelatedAngleEnergy), intent(in) :: this + real(8), intent(in) :: E_in + real(8), intent(out) :: E_out + real(8), intent(out) :: mu + + ! Sample cosine of scattering angle + mu = this%angle%sample(E_in) + + ! Sample outgoing energy + E_out = this%energy%sample(E_in) + end subroutine uncorrelated_sample + +end module secondary_uncorrelated From dcb675fbacd20c562410088af72db2c7e3798471 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 17 Nov 2015 15:30:18 -0600 Subject: [PATCH 02/11] Full use of new object-oriented secondary angle-energy distributions --- src/ace.F90 | 764 ++++++++++++++----------- src/ace_header.F90 | 77 +-- src/distribution_univariate.F90 | 29 +- src/endf_header.F90 | 34 ++ src/energy_distribution.F90 | 27 +- src/interpolation.F90 | 1 - src/output.F90 | 31 +- src/physics.F90 | 962 +------------------------------- src/search.F90 | 1 - src/secondary_correlated.F90 | 62 +- src/secondary_header.F90 | 19 +- src/secondary_kalbach.F90 | 18 +- src/secondary_uncorrelated.F90 | 13 +- 13 files changed, 581 insertions(+), 1457 deletions(-) diff --git a/src/ace.F90 b/src/ace.F90 index e97338b558..1407a4b132 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -1,17 +1,23 @@ module ace - use ace_header, only: Nuclide, Reaction, SAlphaBeta, XsListing, & - DistEnergy + use ace_header, only: Nuclide, Reaction, SAlphaBeta, XsListing use constants - use endf, only: reaction_name, is_fission, is_disappearance - use error, only: fatal_error, warning - use fission, only: nu_total + use distribution_univariate, only: Uniform, Equiprobable, Tabular + use endf, only: reaction_name, is_fission, is_disappearance + use energy_distribution, only: TabularEquiprobable, LevelInelastic, & + ContinuousTabular, MaxwellEnergy, Evaporation, WattEnergy, NBodyPhaseSpace + use error, only: fatal_error, warning + use fission, only: nu_total use global - use list_header, only: ListInt - use material_header, only: Material - use output, only: write_message - use set_header, only: SetChar - use string, only: to_str, to_lower + use list_header, only: ListInt + use material_header, only: Material + use output, only: write_message + use set_header, only: SetChar + use secondary_header, only: AngleEnergy + use secondary_correlated, only: CorrelatedAngleEnergy + use secondary_kalbach, only: KalbachMann + use secondary_uncorrelated, only: UncorrelatedAngleEnergy + use string, only: to_str, to_lower implicit none @@ -372,8 +378,8 @@ contains else call read_nu_data(nuc) call read_reactions(nuc) - call read_angular_dist(nuc) call read_energy_dist(nuc) + call read_angular_dist(nuc) call read_unr_res(nuc) end if @@ -516,9 +522,10 @@ contains integer :: LED ! location of energy distribution locators integer :: LDIS ! location of all energy distributions integer :: LOCC ! location of energy distributions for given MT + integer :: LAW + integer :: IDAT integer :: lc ! locator integer :: length ! length of data to allocate - type(DistEnergy), pointer :: edist JXS2 = JXS(2) JXS24 = JXS(24) @@ -661,11 +668,15 @@ contains ! Loop over all delayed neutron precursor groups do i = 1, NPCR ! find location of energy distribution data - LOCC = int(XSS(LED + i - 1)) + LOCC = nint(XSS(LED + i - 1)) + + ! Determine law and location of data + LAW = nint(XSS(LDIS + LOCC)) + IDAT = nint(XSS(LDIS + LOCC + 1)) ! read energy distribution data - edist => nuc % nu_d_edist(i) - call get_energy_dist(edist, LOCC, .true.) + call get_energy_dist(nuc%nu_d_edist(i)%obj, LAW, LDIS, IDAT, & + ZERO, ZERO) end do ! ======================================================================= @@ -738,8 +749,8 @@ contains rxn%multiplicity = 1 rxn%threshold = 1 rxn%scatter_in_cm = .true. - rxn%has_angle_dist = .false. - rxn%has_energy_dist = .false. + allocate(rxn%secondary%distribution(1)) + allocate(UncorrelatedAngleEnergy :: rxn%secondary%distribution(1)%obj) end associate ! Add contribution of elastic scattering to total cross section @@ -754,10 +765,6 @@ contains do i = 1, NMT associate (rxn => nuc % reactions(i+1)) - ! set defaults - rxn % has_angle_dist = .false. - rxn % has_energy_dist = .false. - ! read MT number, Q-value, and neutrons produced rxn % MT = int(XSS(LMT + i - 1)) rxn % Q_value = XSS(JXS4 + i - 1) @@ -887,86 +894,101 @@ contains subroutine read_angular_dist(nuc) type(Nuclide), intent(inout) :: nuc - integer :: JXS8 ! location of angular distribution locators - integer :: JXS9 ! location of angular distributions integer :: LOCB ! location of angular distribution for given MT integer :: NE ! number of incoming energies integer :: NP ! number of points for cosine distribution - integer :: LC ! locator integer :: i ! index in reactions array integer :: j ! index over incoming energies - integer :: length ! length of data array to allocate - - JXS8 = JXS(8) - JXS9 = JXS(9) + integer :: k ! index over energy distributions + integer :: interp + integer, allocatable :: LC(:) ! locator + real(8), allocatable :: x(:) + real(8), allocatable :: p(:) ! loop over all reactions with secondary neutrons -- NXS(5) does not include ! elastic scattering do i = 1, NXS(5) + 1 associate (rxn => nuc%reactions(i)) - ! find location of angular distribution - LOCB = int(XSS(JXS8 + i - 1)) - if (LOCB == -1) then - ! Angular distribution data are specified through LAWi = 44 in the DLW - ! block - cycle - elseif (LOCB == 0) then - ! No angular distribution data are given for this reaction, isotropic - ! scattering is assumed (in CM if TY < 0 and in LAB if TY > 0) - cycle - end if - rxn % has_angle_dist = .true. + LOCB = int(XSS(JXS(8) + i - 1)) - ! allocate space for incoming energies and locations - NE = int(XSS(JXS9 + LOCB - 1)) - rxn % adist % n_energy = NE - allocate(rxn % adist % energy(NE)) - allocate(rxn % adist % type(NE)) - allocate(rxn % adist % location(NE)) + ! Angular distribution given as part of a correlated angle-energy distribution + if (LOCB == -1) cycle - ! read incoming energy grid and location of nucs - XSS_index = JXS9 + LOCB - rxn % adist % energy = get_real(NE) - rxn % adist % location = get_int(NE) + ! No angular distribution data are given for this reaction, isotropic + ! scattering is assumed (in CM if TY < 0 and in LAB if TY > 0) + if (LOCB == 0) cycle - ! determine dize of data block - length = 0 - do j = 1, NE - LC = rxn % adist % location(j) - if (LC == 0) then - ! isotropic - rxn % adist % type(j) = ANGLE_ISOTROPIC - elseif (LC > 0) then - ! 32 equiprobable bins - rxn % adist % type(j) = ANGLE_32_EQUI - length = length + 33 - elseif (LC < 0) then - ! tabular distribution - rxn % adist % type(j) = ANGLE_TABULAR - NP = int(XSS(JXS9 + abs(LC))) - length = length + 2 + 3*NP - end if - end do + ! Loop over each separate energy distribution. Even though there is only + ! "one" angular distribution, it is repeated as many times as there are + ! energy distributions for this reaction since the + ! UncorrelatedAngleEnergy type holds one angle and energy distribution. + do k = 1, size(rxn%secondary%distribution) + select type (aedist => rxn%secondary%distribution(k)%obj) + type is (UncorrelatedAngleEnergy) + ! allocate space for incoming energies and locations + NE = int(XSS(JXS(9) + LOCB - 1)) + allocate(aedist%angle%energy(NE)) + allocate(aedist%angle%distribution(NE)) + allocate(LC(NE)) - ! allocate angular distribution data and read - allocate(rxn % adist % data(length)) + ! read incoming energy grid and location of nucs + XSS_index = JXS(9) + LOCB + aedist%angle%energy(:) = get_real(NE) + LC(:) = get_int(NE) - ! read angular distribution -- currently this does not actually parse the - ! angular distribution tables for each incoming energy, that must be done - ! on-the-fly - XSS_index = JXS9 + LOCB + 2 * NE - rxn % adist % data = get_real(length) + ! determine dize of data block + do j = 1, NE + if (LC(j) == 0) then + ! isotropic + allocate(Uniform :: aedist%angle%distribution(j)%obj) + select type (adist => aedist%angle%distribution(j)%obj) + type is (Uniform) + adist%a = -ONE + adist%b = ONE + end select - ! change location pointers since they are currently relative to JXS(9) - LC = LOCB + 2 * NE + 1 - do j = 1, NE - ! For consistency, leave location as 0 if type is isotropic. - ! This is not necessary for current correctness, but can avoid - ! future issues - if (rxn % adist % location(j) /= 0) then - rxn % adist % location(j) = abs(rxn % adist % location(j)) - LC - end if + elseif (LC(j) > 0) then + ! 32 equiprobable bins + allocate(Equiprobable :: aedist%angle%distribution(j)%obj) + select type (adist => aedist%angle%distribution(j)%obj) + type is (Equiprobable) + allocate(adist%x(33)) + end select + + elseif (LC(j) < 0) then + ! tabular distribution + allocate(Tabular :: aedist%angle%distribution(j)%obj) + end if + end do + + ! read angular distribution -- currently this does not actually parse the + ! angular distribution tables for each incoming energy, that must be done + ! on-the-fly + do j = 1, NE + XSS_index = JXS(9) + abs(LC(j)) - 1 + select type(adist => aedist%angle%distribution(j)%obj) + type is (Equiprobable) + adist%x(:) = get_real(33) + type is (Tabular) + ! determine interpolation and number of points + interp = nint(XSS(XSS_index)) + NP = nint(XSS(XSS_index + 1)) + + ! Get probability density data + XSS_index = XSS_index + 2 + allocate(x(NP), p(NP)) + x(:) = get_real(NP) + p(:) = get_real(NP) + + ! initialize distribution + call adist%initialize(x, p, interp) + deallocate(x, p) + end select + end do + deallocate(LC) + + end select end do end associate end do @@ -981,25 +1003,46 @@ contains subroutine read_energy_dist(nuc) type(Nuclide), intent(inout) :: nuc - integer :: LED ! location of energy distribution locators - integer :: LOCC ! location of energy distributions for given MT integer :: i ! loop index - - LED = JXS(10) + integer :: n + integer :: IDAT ! locator for distribution data + integer :: LNW ! location of next energy law + integer :: LAW ! Type of energy law ! Loop over all reactions do i = 1, NXS(5) - associate (rxn => nuc % reactions(i+1)) ! skip over elastic scattering - rxn % has_energy_dist = .true. + ! Determine how many energy distributions are present for this reaction + LNW = nint(XSS(JXS(10) + i - 1)) + n = 0 + do while (LNW > 0) + n = n + 1 + LNW = nint(XSS(JXS(11) + LNW - 1)) + end do - ! find location of energy distribution data - LOCC = int(XSS(LED + i - 1)) + ! Allocate space for distributions and probability of validity + associate (secondary => nuc%reactions(i + 1)%secondary) + allocate(secondary%applicability(n)) + allocate(secondary%distribution(n)) - ! allocate energy distribution - allocate(rxn % edist) + LNW = nint(XSS(JXS(10) + i - 1)) + n = 0 + do while (LNW > 0) + n = n + 1 - ! read data for energy distribution - call get_energy_dist(rxn % edist, LOCC) + ! Determine energy law and location of data + LAW = nint(XSS(JXS(11) + LNW)) + IDAT = nint(XSS(JXS(11) + LNW + 1)) + + ! Read probability of law validity + call secondary%applicability(n)%from_ace(XSS, JXS(11) + LNW + 2) + + ! Read energy law data + call get_energy_dist(secondary%distribution(n)%obj, LAW, & + JXS(11), IDAT, nuc%awr, nuc%reactions(i)%Q_value) + + ! Get locator for next distribution + LNW = nint(XSS(JXS(11) + LNW - 1)) + end do end associate end do @@ -1011,281 +1054,324 @@ contains ! single reaction !=============================================================================== - recursive subroutine get_energy_dist(edist, loc_law, delayed_n) - type(DistEnergy), intent(inout) :: edist ! energy distribution - integer, intent(in) :: loc_law ! locator for data - logical, intent(in), optional :: delayed_n ! is this for delayed neutrons? + recursive subroutine get_energy_dist(aedist, law, LDIS, IDAT, awr, Q_value) + class(AngleEnergy), allocatable, intent(inout) :: aedist + integer, intent(in) :: law + integer, intent(in) :: LDIS + integer, intent(in) :: IDAT + real(8), intent(in) :: awr + real(8), intent(in) :: Q_value - integer :: LDIS ! location of all energy distributions - integer :: LNW ! location of next energy distribution if multiple - integer :: LAW ! secondary energy distribution law + integer :: i, j integer :: NR ! number of interpolation regions integer :: NE ! number of incoming energies - integer :: IDAT ! location of first energy distribution for given MT - integer :: lc ! locator - integer :: length ! length of data to allocate - integer :: length_interp_data ! length of interpolation data + integer :: NP ! number of outgoing energies/angles + integer :: interp + integer, allocatable :: L(:) ! locations of distributions for each Ein + integer, allocatable :: LC(:) ! locations of distributions for each Ein + real(8), allocatable :: x(:) + real(8), allocatable :: p(:) - ! determine location of energy distribution - if (present(delayed_n)) then - LDIS = JXS(27) + XSS_index = LDIS + IDAT - 1 + + if (law == 44) then + allocate(KalbachMann :: aedist) + elseif (law == 61) then + allocate(CorrelatedAngleEnergy :: aedist) else - LDIS = JXS(11) + allocate(UncorrelatedAngleEnergy :: aedist) end if - ! locator for next law and information on this law - LNW = int(XSS(LDIS + loc_law - 1)) - LAW = int(XSS(LDIS + loc_law)) - IDAT = int(XSS(LDIS + loc_law + 1)) - NR = int(XSS(LDIS + loc_law + 2)) - edist % law = LAW - edist % p_valid % n_regions = NR + select type (aedist) + type is (UncorrelatedAngleEnergy) + ! ======================================================================== + ! UNCORRELATED ENERGY DISTRIBUTIONS - ! allocate space for ENDF interpolation parameters - if (NR > 0) then - allocate(edist % p_valid % nbt(NR)) - allocate(edist % p_valid % int(NR)) - end if - - ! read ENDF interpolation parameters - XSS_index = LDIS + loc_law + 3 - if (NR > 0) then - edist % p_valid % nbt = int(get_real(NR)) - edist % p_valid % int = int(get_real(NR)) - end if - - ! allocate space for law validity data - NE = int(XSS(LDIS + loc_law + 3 + 2*NR)) - edist % p_valid % n_pairs = NE - allocate(edist % p_valid % x(NE)) - allocate(edist % p_valid % y(NE)) - - length_interp_data = 5 + 2*(NR + NE) - - ! read law validity data - XSS_index = LDIS + loc_law + 4 + 2*NR - edist % p_valid % x = get_real(NE) - edist % p_valid % y = get_real(NE) - - ! Set index to beginning of IDAT array - lc = LDIS + IDAT - 2 - - ! determine length of energy distribution - length = length_energy_dist(lc, LAW, loc_law, length_interp_data) - - ! allocate secondary energy distribution array - allocate(edist % data(length)) - - ! read secondary energy distribution - XSS_index = lc + 1 - edist % data = get_real(length) - - ! read next energy distribution if present - if (LNW > 0) then - allocate(edist % next) - call get_energy_dist(edist % next, LNW) - end if - - end subroutine get_energy_dist - -!=============================================================================== -! LENGTH_ENERGY_DIST determines how many values are contained in an LDAT energy -! distribution array based on the secondary energy law and location in XSS -!=============================================================================== - - function length_energy_dist(lc, law, LOCC, lid) result(length) - integer, intent(in) :: lc ! location in XSS array - integer, intent(in) :: law ! energy distribution law - integer, intent(in) :: LOCC ! location of energy distribution - integer, intent(in) :: lid ! length of interpolation data - integer :: length ! length of energy distribution (LDAT) - - integer :: i ! loop index for incoming energies - integer :: j ! loop index for outgoing energies - integer :: k ! dummy index in XSS - integer :: NR ! number of interpolation regions - integer :: NE ! number of incoming energies - integer :: NP ! number of points in outgoing energy distribution - integer :: NMU ! number of points in outgoing cosine distribution - integer :: NRa ! number of interpolation regions for Watt 'a' - integer :: NEa ! number of energies for Watt 'a' - integer :: NRb ! number of interpolation regions for Watt 'b' - integer :: NEb ! number of energies for Watt 'b' - real(8), allocatable :: L(:) ! locations of distributions for each Ein - - ! initialize length - length = 0 - - select case (law) - case (1) - ! Tabular equiprobable energy bins - NR = int(XSS(lc + 1)) - NE = int(XSS(lc + 2 + 2*NR)) - NP = int(XSS(lc + 3 + 2*NR + NE)) - length = 3 + 2*NR + NE + 3*NP*NE - - case (2) - ! Discrete photon energy - length = 2 - - case (3) - ! Level scattering - length = 2 - - case (4) - ! Continuous tabular distribution - NR = int(XSS(lc + 1)) - NE = int(XSS(lc + 2 + 2*NR)) - allocate(L(NE)) - L(:) = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1)) - - ! Continue with finding data length - length = length + 2 + 2*NR + 2*NE - do i = 1,NE - ! Some older data sets use the same LDAT for multiple Ein tables. - ! If this is the case, we should skip incrementing length when it is - ! not needed. - if (i < NE) then - if (any(L(i) == L(i + 1: NE))) then - ! adjust location for this block - j = lc + 2 + 2*NR + NE + i - XSS(j) = XSS(j) - LOCC - lid - cycle + select case (law) + case (1) + allocate(TabularEquiprobable :: aedist%energy) + select type (edist => aedist%energy) + type is (TabularEquiprobable) + NR = nint(XSS(XSS_index)) + NE = nint(XSS(XSS_index + 1 + 2*NR)) + if (NR > 0) then + call fatal_error("Multiple interpolation regions not yet supported & + &for tabular equiprobable energy distributions.") end if - end if - ! determine length - NP = int(XSS(lc + length + 2)) - length = length + 2 + 3*NP + edist%n_region = NR - ! adjust location for this block - j = lc + 2 + 2*NR + NE + i - XSS(j) = XSS(j) - LOCC - lid + ! Read incoming energies for which outgoing energies are tabulated + allocate(edist%energy_in(NE)) + XSS_index = XSS_index + 2 + 2*NR + edist%energy_in(:) = get_real(NE) + + ! Read outgoing energy tables + NP = nint(XSS(XSS_index)) + allocate(edist%energy_out(NP, NE)) + XSS_index = XSS_index + 1 + do i = 1, NE + edist%energy_out(:, i) = get_real(NP) + end do + end select + + case (3) + allocate(LevelInelastic :: aedist%energy) + select type (edist => aedist%energy) + type is (LevelInelastic) + edist%threshold = XSS(XSS_index) + edist%mass_ratio = XSS(XSS_index + 1) + end select + + case (4) + allocate(ContinuousTabular :: aedist%energy) + select type (edist => aedist%energy) + type is (ContinuousTabular) + NR = nint(XSS(XSS_index)) + XSS_index = XSS_index + 1 + if (NR > 1) then + call fatal_error("Multiple interpolation regions not yet supported & + &for continuous tabular energy distributions.") + end if + edist%n_region = NR + + ! Read breakpoints and interpolation parameters + if (NR > 0) then + allocate(edist%breakpoints(NR)) + allocate(edist%interpolation(NR)) + edist%breakpoints(:) = get_int(NR) + edist%interpolation(:) = get_int(NR) + end if + + ! Read incoming energies for which outgoing energies are tabulated and + ! locators + NE = nint(XSS(XSS_index)) + XSS_index = XSS_index + 1 + allocate(edist%energy_in(NE)) + allocate(L(NE)) + edist%energy_in(:) = get_real(NE) + L(:) = get_int(NE) + + ! Read outgoing energy tables + allocate(edist%energy_out(NE)) + do i = 1, NE + ! Determine interpolation and number of discrete points + XSS_index = LDIS + L(i) - 1 + interp = nint(XSS(XSS_index)) + edist%energy_out(i)%interpolation = mod(interp, 10) + edist%energy_out(i)%n_discrete = (interp - & + edist%energy_out(i)%interpolation)/10 + + ! check for discrete lines present + if (edist%energy_out(i)%n_discrete > 0) then + call fatal_error("Discrete lines in continuous tabular & + &distribution not yet supported") + end if + + ! Determine number of points and allocate space + NP = nint(XSS(XSS_index + 1)) + allocate(edist%energy_out(i)%e_out(NP)) + allocate(edist%energy_out(i)%p(NP)) + allocate(edist%energy_out(i)%c(NP)) + + ! Read tabular PDF for outgoing energy + XSS_index = XSS_index + 2 + edist%energy_out(i)%e_out(:) = get_real(NP) + edist%energy_out(i)%p(:) = get_real(NP) + edist%energy_out(i)%c(:) = get_real(NP) + end do + + deallocate(L) + end select + + case (7) + allocate(MaxwellEnergy :: aedist%energy) + select type (edist => aedist%energy) + type is (MaxwellEnergy) + call edist%theta%from_ace(XSS, XSS_index) + edist%u = XSS(XSS_index + 2 + 2*edist%theta%n_regions + & + 2*edist%theta%n_pairs) + end select + + case (9) + allocate(Evaporation :: aedist%energy) + select type(edist => aedist%energy) + type is (Evaporation) + call edist%theta%from_ace(XSS, XSS_index) + edist%u = XSS(XSS_index + 2 + 2*edist%theta%n_regions + & + 2*edist%theta%n_pairs) + end select + + case (11) + allocate(WattEnergy :: aedist%energy) + select type(edist => aedist%energy) + type is (WattEnergy) + call edist%a%from_ace(XSS, XSS_index) + XSS_index = XSS_index + 2 + 2*edist%a%n_regions + 2*edist%a%n_pairs + call edist%b%from_ace(XSS, XSS_index) + XSS_index = XSS_index + 2 + 2*edist%b%n_regions + 2*edist%b%n_pairs + edist%u = XSS(XSS_index) + end select + + case (66) + allocate(NBodyPhaseSpace :: aedist%energy) + select type(edist => aedist%energy) + type is (NBodyPhaseSpace) + edist%n_bodies = int(XSS(XSS_index)) + edist%mass_ratio = XSS(XSS_index + 1) + edist%A = awr + edist%Q = Q_value + end select + + end select + + type is (KalbachMann) + ! ======================================================================== + ! CORRELATED KALBACH-MANN DISTRIBUTION + + NR = int(XSS(XSS_index)) + NE = int(XSS(XSS_index + 1 + 2*NR)) + if (NR > 0) then + call fatal_error("Multiple interpolation regions not yet supported & + &for Kalbach-Mann energy distributions.") + end if + aedist%n_region = NR + + ! Read incoming energies for which outgoing energies are tabulated and locators + allocate(aedist%energy_in(NE)) + allocate(L(NE)) + XSS_index = XSS_index + 2 + 2*NR + aedist%energy_in(:) = get_real(NE) + L(:) = get_int(NE) + + ! Read outgoing energy tables + allocate(aedist%table(NE)) + do i = 1, NE + ! Determine interpolation and number of discrete points + XSS_index = LDIS + L(i) - 1 + interp = nint(XSS(XSS_index)) + aedist%table(i)%interpolation = mod(interp, 10) + aedist%table(i)%n_discrete = (interp - aedist%table(i)%interpolation)/10 + + ! check for discrete lines present + if (aedist%table(i)%n_discrete > 0) then + call fatal_error("Discrete lines in Kalbach-Mann distribution not & + &yet supported") + end if + + ! Determine number of points and allocate space + NP = nint(XSS(XSS_index + 1)) + allocate(aedist%table(i)%e_out(NP)) + allocate(aedist%table(i)%p(NP)) + allocate(aedist%table(i)%c(NP)) + allocate(aedist%table(i)%r(NP)) + allocate(aedist%table(i)%a(NP)) + + ! Read tabular PDF for outgoing energy + XSS_index = XSS_index + 2 + aedist%table(i)%e_out(:) = get_real(NP) + aedist%table(i)%p(:) = get_real(NP) + aedist%table(i)%c(:) = get_real(NP) + aedist%table(i)%r(:) = get_real(NP) + aedist%table(i)%a(:) = get_real(NP) end do + deallocate(L) - case (5) - ! General evaporation spectrum - NR = int(XSS(lc + 1)) - NE = int(XSS(lc + 2 + 2*NR)) - NP = int(XSS(lc + 3 + 2*NR + 2*NE)) - length = 3 + 2*NR + 2*NE + NP + type is (CorrelatedAngleEnergy) + ! ======================================================================== + ! CORRELATED ANGLE-ENERGY DISTRIBUTION - case (7) - ! Maxwell fission spectrum - NR = int(XSS(lc + 1)) - NE = int(XSS(lc + 2 + 2*NR)) - length = 3 + 2*NR + 2*NE + NR = int(XSS(XSS_index)) + NE = int(XSS(XSS_index + 1 + 2*NR)) + if (NR > 0) then + call fatal_error("Multiple interpolation regions not yet supported & + &for correlated angle-energy distributions.") + end if + aedist%n_region = NR - case (9) - ! Evaporation spectrum - NR = int(XSS(lc + 1)) - NE = int(XSS(lc + 2 + 2*NR)) - length = 3 + 2*NR + 2*NE - - case (11) - ! Watt spectrum - NRa = int(XSS(lc + 1)) - NEa = int(XSS(lc + 2 + 2*NRa)) - NRb = int(XSS(lc + 3 + 2*(NRa+NEa))) - NEb = int(XSS(lc + 4 + 2*(NRa+NEa+NRb))) - length = 5 + 2*(NRa + NEa + NRb + NEb) - - case (44) - ! Kalbach-Mann correlated scattering - NR = int(XSS(lc + 1)) - NE = int(XSS(lc + 2 + 2*NR)) + ! Read incoming energies for which outgoing energies are tabulated and + ! locators + allocate(aedist%energy_in(NE)) allocate(L(NE)) - L(:) = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1)) + XSS_index = XSS_index + 2 + 2*NR + aedist%energy_in(:) = get_real(NE) + L(:) = get_int(NE) - ! Continue with finding data length - length = length + 2 + 2*NR + 2*NE - do i = 1,NE - ! Some older data sets use the same LDAT for multiple Ein tables. - ! If this is the case, we should skip incrementing length when it is - ! not needed. - if (i < NE) then - if (any(L(i) == L(i + 1: NE))) then - ! adjust location for this block - j = lc + 2 + 2*NR + NE + i - XSS(j) = XSS(j) - LOCC - lid - cycle - end if - end if - NP = int(XSS(lc + length + 2)) - length = length + 2 + 5*NP + ! Read outgoing energy tables + allocate(aedist%table(NE)) + do i = 1, NE + ! Determine interpolation and number of discrete points + XSS_index = LDIS + L(i) - 1 + interp = nint(XSS(XSS_index)) + aedist%table(i)%interpolation = mod(interp, 10) + aedist%table(i)%n_discrete = (interp - aedist%table(i)%interpolation)/10 - ! adjust location for this block - j = lc + 2 + 2*NR + NE + i - XSS(j) = XSS(j) - LOCC - lid - end do - deallocate(L) - - case (61) - ! Correlated energy and angle distribution - NR = int(XSS(lc + 1)) - NE = int(XSS(lc + 2 + 2*NR)) - allocate(L(NE)) - L(:) = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1)) - - ! Continue with finding data length - length = length + 2 + 2*NR + 2*NE - do i = 1,NE - ! Some older data sets use the same LDAT for multiple Ein tables. - ! If this is the case, we should skip incrementing length when it is - ! not needed. - if (i < NE) then - if (any(L(i) == L(i + 1: NE))) then - ! adjust locators for energy distribution - j = lc + 2 + 2*NR + NE + i - XSS(j) = XSS(j) - LOCC - lid - cycle - end if + ! check for discrete lines present + if (aedist%table(i)%n_discrete > 0) then + call fatal_error("Discrete lines in correlated angle-energy & + &distribution not yet supported") end if - ! outgoing energy distribution - NP = int(XSS(lc + length + 2)) + ! Determine number of points and allocate space + NP = nint(XSS(XSS_index + 1)) + allocate(aedist%table(i)%e_out(NP)) + allocate(aedist%table(i)%p(NP)) + allocate(aedist%table(i)%c(NP)) + allocate(LC(NP)) - ! adjust locators for angular distribution + ! Read tabular PDF for outgoing energy + XSS_index = XSS_index + 2 + aedist%table(i)%e_out(:) = get_real(NP) + aedist%table(i)%p(:) = get_real(NP) + aedist%table(i)%c(:) = get_real(NP) + LC(:) = get_int(NP) + + ! allocate angular distributions for each incoming/outgoing energy + allocate(aedist%table(i)%angle(NP)) do j = 1, NP - k = lc + length + 2 + 3*NP + j - if (XSS(k) /= 0) XSS(k) = XSS(k) - LOCC - lid + if (LC(j) == 0) then + ! isotropic + allocate(Uniform :: aedist%table(i)%angle(j)%obj) + select type (adist => aedist%table(i)%angle(j)%obj) + type is (Uniform) + adist%a = -ONE + adist%b = ONE + end select + + elseif (LC(j) > 0) then + ! tabular distribution + allocate(Tabular :: aedist%table(i)%angle(j)%obj) + end if end do - length = length + 2 + 4*NP + ! read angular distributions do j = 1, NP - ! outgoing angle distribution -- NMU here is actually - ! referred to as NP in the MCNP documentation - NMU = int(XSS(lc + length + 2)) - length = length + 2 + 3*NMU + XSS_index = LDIS + abs(LC(j)) - 1 + select type(adist => aedist%table(i)%angle(j)%obj) + type is (Tabular) + ! determine interpolation and number of points + interp = nint(XSS(XSS_index)) + NP = nint(XSS(XSS_index + 1)) + + ! Get probability density data + XSS_index = XSS_index + 2 + allocate(x(NP), p(NP)) + x(:) = get_real(NP) + p(:) = get_real(NP) + + ! initialize distribution + call adist%initialize(x, p, interp) + deallocate(x, p) + end select end do + deallocate(LC) - ! adjust locators for energy distribution - j = lc + 2 + 2*NR + NE + i - XSS(j) = XSS(j) - LOCC - lid end do - deallocate(L) - case (66) - ! N-body phase space distribution - length = 2 - case (67) - ! Laboratory energy-angle law - NR = int(XSS(lc + 1)) - NE = int(XSS(lc + 2 + 2*NR)) - ! Before progressing, check to see if data set uses L(I) values - ! in a way inconsistent with the current form of the ACE Format Guide - ! (MCNP5 Manual, Vol 3) - allocate(L(NE)) - L(:) = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1)) - ! Don't currently do anything with L deallocate(L) - ! Continue with finding data length - NMU = int(XSS(lc + 4 + 2*NR + 2*NE)) - length = 4 + 2*(NR + NE + NMU) - end select - end function length_energy_dist + end subroutine get_energy_dist !=============================================================================== ! READ_UNR_RES reads in unresolved resonance probability tables if present. diff --git a/src/ace_header.F90 b/src/ace_header.F90 index d7c298979c..2331f9117c 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -3,42 +3,11 @@ module ace_header use constants, only: MAX_FILE_LEN, ZERO use dict_header, only: DictIntInt use endf_header, only: Tab1 + use secondary_header, only: SecondaryDistribution, AngleEnergyContainer use stl_vector, only: VectorInt implicit none -!=============================================================================== -! DISTANGLE contains data for a tabular secondary angle distribution whether it -! be tabular or 32 equiprobable cosine bins -!=============================================================================== - - type DistAngle - integer :: n_energy ! # of incoming energies - real(8), allocatable :: energy(:) ! incoming energy grid - integer, allocatable :: type(:) ! type of distribution - integer, allocatable :: location(:) ! location of each table - real(8), allocatable :: data(:) ! angular distribution data - end type DistAngle - -!=============================================================================== -! DISTENERGY contains data for a secondary energy distribution for all -! scattering laws -!=============================================================================== - - type DistEnergy - integer :: law ! secondary distribution law - type(Tab1) :: p_valid ! probability of law validity - real(8), allocatable :: data(:) ! energy distribution data - - ! For reactions that may have multiple energy distributions such as (n,2n), - ! this pointer allows multiple laws to be stored - type(DistEnergy), pointer :: next => null() - - ! Type-Bound procedures - contains - procedure :: clear => distenergy_clear ! Deallocates DistEnergy - end type DistEnergy - !=============================================================================== ! REACTION contains the cross-section and secondary energy and angle ! distributions for a single reaction in a continuous-energy ACE-format table @@ -53,10 +22,7 @@ module ace_header logical :: scatter_in_cm ! scattering system in center-of-mass? logical :: multiplicity_with_E = .false. ! Flag to indicate E-dependent multiplicity real(8), allocatable :: sigma(:) ! Cross section values - logical :: has_angle_dist ! Angle distribution present? - logical :: has_energy_dist ! Energy distribution present? - type(DistAngle) :: adist ! Secondary angular distribution - type(DistEnergy), pointer :: edist => null() ! Secondary energy distribution + type(SecondaryDistribution) :: secondary ! Type-Bound procedures contains @@ -137,7 +103,7 @@ module ace_header integer :: n_precursor ! # of delayed neutron precursors real(8), allocatable :: nu_d_data(:) real(8), allocatable :: nu_d_precursor_data(:) - type(DistEnergy), pointer :: nu_d_edist(:) => null() + type(AngleEnergyContainer), allocatable :: nu_d_edist(:) ! Unresolved resonance data logical :: urr_present @@ -284,37 +250,14 @@ module ace_header contains -!=============================================================================== -! DISTENERGY_CLEAR resets and deallocates data in DistEnergy. -!=============================================================================== - - recursive subroutine distenergy_clear(this) - - class(DistEnergy), intent(inout) :: this ! The DistEnergy object to clear - - if (associated(this % next)) then - ! recursively clear this item - call this % next % clear() - deallocate(this % next) - end if - - end subroutine distenergy_clear - !=============================================================================== ! REACTION_CLEAR resets and deallocates data in Reaction. !=============================================================================== subroutine reaction_clear(this) - class(Reaction), intent(inout) :: this ! The Reaction object to clear if (associated(this % multiplicity_E)) deallocate(this % multiplicity_E) - - if (associated(this % edist)) then - call this % edist % clear() - deallocate(this % edist) - end if - end subroutine reaction_clear !=============================================================================== @@ -322,21 +265,11 @@ module ace_header !=============================================================================== subroutine nuclide_clear(this) - - class(Nuclide), intent(inout) :: this ! The Nuclide object to clear + class(Nuclide), intent(inout) :: this integer :: i ! Loop counter - if (associated(this % nu_d_edist)) then - do i = 1, size(this % nu_d_edist) - call this % nu_d_edist(i) % clear() - end do - deallocate(this % nu_d_edist) - end if - - if (associated(this % urr_data)) then - deallocate(this % urr_data) - end if + if (associated(this % urr_data)) deallocate(this % urr_data) if (allocated(this % reactions)) then do i = 1, size(this % reactions) diff --git a/src/distribution_univariate.F90 b/src/distribution_univariate.F90 index b3d77e2f82..f596536d0b 100644 --- a/src/distribution_univariate.F90 +++ b/src/distribution_univariate.F90 @@ -1,7 +1,7 @@ module distribution_univariate - use constants, only: ZERO, HALF, HISTOGRAM, LINEAR_LINEAR, MAX_LINE_LEN, & - MAX_WORD_LEN + use constants, only: ZERO, ONE, HALF, HISTOGRAM, LINEAR_LINEAR, & + MAX_LINE_LEN, MAX_WORD_LEN use error, only: fatal_error use math, only: maxwell_spectrum, watt_spectrum use random_lcg, only: prn @@ -78,6 +78,12 @@ module distribution_univariate procedure :: initialize => tabular_initialize end type Tabular + type, extends(Distribution) :: Equiprobable + real(8), allocatable :: x(:) + contains + procedure :: sample => equiprobable_sample + end type Equiprobable + contains function discrete_sample(this) result(x) @@ -238,6 +244,25 @@ contains this%c(:) = this%c(:)/this%c(n) end subroutine tabular_initialize + function equiprobable_sample(this) result(x) + class(Equiprobable), intent(in) :: this + real(8) :: x + + integer :: i + integer :: n + real(8) :: r + real(8) :: xl, xr + + n = size(this%x) + + r = prn() + i = 1 + int((n - 1)*r) + + xl = this%x(i) + xr = this%x(i+1) + x = xl + ((n - 1)*r - i + ONE) * (xr - xl) + end function equiprobable_sample + subroutine distribution_from_xml(dist, node_dist) class(Distribution), allocatable, intent(inout) :: dist type(Node), pointer :: node_dist diff --git a/src/endf_header.F90 b/src/endf_header.F90 index af62231a5d..7388ea2f54 100644 --- a/src/endf_header.F90 +++ b/src/endf_header.F90 @@ -13,6 +13,40 @@ module endf_header integer :: n_pairs ! # of pairs of (x,y) values real(8), allocatable :: x(:) ! values of abscissa real(8), allocatable :: y(:) ! values of ordinate + contains + procedure :: from_ace end type Tab1 +contains + + subroutine from_ace(this, xss, idx) + class(Tab1), intent(inout) :: this + real(8), intent(in) :: xss(:) + integer, intent(in) :: idx + + integer :: nr, ne + + ! Determine number of regions + nr = nint(xss(idx)) + this%n_regions = nr + + ! Read interpolation region data + if (nr > 0) then + allocate(this%nbt(nr)) + allocate(this%int(nr)) + this%nbt(:) = nint(xss(idx + 1 : idx + nr)) + this%int(:) = nint(xss(idx + nr + 1 : idx + 2*nr)) + end if + + ! Determine number of pairs + ne = int(XSS(idx + 2*nr + 1)) + this%n_pairs = ne + + ! Read (x,y) pairs + allocate(this%x(ne)) + allocate(this%y(ne)) + this%x(:) = xss(idx + 2*nr + 2 : idx + 2*nr + 1 + ne) + this%y(:) = xss(idx + 2*nr + 2 + ne : idx + 2*nr + 1 + 2*ne) + end subroutine from_ace + end module endf_header diff --git a/src/energy_distribution.F90 b/src/energy_distribution.F90 index 2a616e0f1f..a8b975879f 100644 --- a/src/energy_distribution.F90 +++ b/src/energy_distribution.F90 @@ -2,7 +2,6 @@ module energy_distribution use constants, only: ZERO, ONE, TWO, PI, HISTOGRAM, LINEAR_LINEAR use endf_header, only: Tab1 - use error, only: fatal_error use interpolation, only: interpolate_tab1 use math, only: maxwell_spectrum, watt_spectrum use random_lcg, only: prn @@ -27,14 +26,14 @@ module energy_distribution end function iSampleEnergy end interface + type :: EnergyDistributionContainer + class(EnergyDistribution), allocatable :: obj + end type EnergyDistributionContainer + !=============================================================================== ! Derived classes !=============================================================================== - type Array1D - real(8), allocatable :: data(:) - end type Array1D - type, extends(EnergyDistribution) :: TabularEquiprobable integer :: n_region integer, allocatable :: breakpoints(:) @@ -121,14 +120,6 @@ contains n_energy_in = size(this%energy_in) n_energy_out = size(this%energy_out, 1) - ! read number of interpolation regions, incoming energies, and outgoing - ! energies - ! TODO: Move this error to input - if (this%n_region > 0) then - call fatal_error("Multiple interpolation regions not supported while & - &attempting to sample equiprobable energy bins.") - end if - ! determine index on incoming energy grid and interpolation factor i = binary_search(this%energy_in, size(this%energy_in), E_in) r = (E_in - this%energy_in(i)) / & @@ -200,9 +191,6 @@ contains ! read number of interpolation regions and incoming energies if (this%n_region == 1) then histogram_interp = (this%interpolation(1) == 1) - else if (this%n_region > 1) then - call fatal_error("Multiple interpolation regions not supported while & - &attempting to sample continuous tabular distribution.") else histogram_interp = .false. end if @@ -245,13 +233,6 @@ contains E_1 = E_i_1 + r*(E_i1_1 - E_i_1) E_K = E_i_K + r*(E_i1_K - E_i_K) - ! TODO: Write error at initizliation - if (this%energy_out(l)%n_discrete > 0) then - ! discrete lines present - call fatal_error("Discrete lines in continuous tabular distributed not & - &yet supported") - end if - ! determine outgoing energy bin n_energy_out = size(this%energy_out(l)%e_out) r1 = prn() diff --git a/src/interpolation.F90 b/src/interpolation.F90 index 9e28bc086e..5f87870677 100644 --- a/src/interpolation.F90 +++ b/src/interpolation.F90 @@ -2,7 +2,6 @@ module interpolation use constants use endf_header, only: Tab1 - use error, only: fatal_error use search, only: binary_search use string, only: to_str diff --git a/src/output.F90 b/src/output.F90 index 66e6831371..3a48de02da 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -359,21 +359,24 @@ contains write(unit_,*) ' Reaction Q-value COM Law IE size(angle) size(energy)' do i = 1, nuc % n_reaction associate (rxn => nuc % reactions(i)) - ! Determine size of angle distribution - if (rxn % has_angle_dist) then - size_angle = rxn % adist % n_energy * 16 + size(rxn % adist % data) * 8 - else - size_angle = 0 - end if +!!$ ! Determine size of angle distribution +!!$ if (rxn % has_angle_dist) then +!!$ size_angle = rxn % adist % n_energy * 16 + size(rxn % adist % data) * 8 +!!$ else +!!$ size_angle = 0 +!!$ end if + size_angle = 0 - ! Determine size of energy distribution and law - if (rxn % has_energy_dist) then - size_energy = size(rxn % edist % data) * 8 - law = to_str(rxn % edist % law) - else - size_energy = 0 - law = 'None' - end if +!!$ ! Determine size of energy distribution and law +!!$ if (rxn % has_energy_dist) then +!!$ size_energy = size(rxn % edist % data) * 8 +!!$ law = to_str(rxn % edist % law) +!!$ else +!!$ size_energy = 0 +!!$ law = 'None' +!!$ end if + size_energy = 0 + law = 'None' write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,A4,1X,I6,1X,I11,1X,I11)') & reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, & diff --git a/src/physics.F90 b/src/physics.F90 index 217d232182..8519f4d2b0 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -1,6 +1,6 @@ module physics - use ace_header, only: Nuclide, Reaction, DistEnergy + use ace_header, only: Nuclide, Reaction use constants use cross_section, only: elastic_xs_0K use endf, only: reaction_name @@ -17,12 +17,10 @@ module physics use random_lcg, only: prn use search, only: binary_search use string, only: to_str + use secondary_uncorrelated, only: UncorrelatedAngleEnergy implicit none -! TODO: Figure out how to write particle restart files in sample_angle, -! sample_energy, etc. - contains !=============================================================================== @@ -458,7 +456,10 @@ contains vel = sqrt(dot_product(v_n, v_n)) ! Sample scattering angle - mu_cm = sample_angle(rxn, E) + select type (dist => rxn%secondary%distribution(1)%obj) + type is (UncorrelatedAngleEnergy) + mu_cm = dist%angle%sample(E) + end select ! Determine direction cosines in CM uvw_cm = v_n/vel @@ -1195,16 +1196,13 @@ contains integer :: lc ! index before start of energies/nu values integer :: NR ! number of interpolation regions integer :: NE ! number of energies tabulated - integer :: law ! energy distribution law integer :: n_sample ! number of times resampling real(8) :: nu_t ! total nu real(8) :: nu_d ! delayed nu - real(8) :: mu ! fission neutron angular cosine real(8) :: beta ! delayed neutron fraction real(8) :: xi ! random number real(8) :: yield ! delayed neutron precursor yield real(8) :: prob ! cumulative probability - type(DistEnergy), pointer :: edist ! Determine total nu nu_t = nu_total(nuc, p % E) @@ -1248,18 +1246,13 @@ contains ! set the delayed group for the particle born from fission p % delayed_group = j - ! select energy distribution for group j - law = nuc % nu_d_edist(j) % law - edist => nuc % nu_d_edist(j) - ! sample from energy distribution n_sample = 0 do - if (law == 44 .or. law == 61) then - call sample_energy(edist, p % E, E_out, mu) - else - call sample_energy(edist, p % E, E_out) - end if + select type (aedist => nuc%nu_d_edist(j)%obj) + type is (UncorrelatedAngleEnergy) + E_out = aedist%energy%sample(p%E) + end select ! resample if energy is greater than maximum neutron energy if (E_out < energy_max_neutron) exit @@ -1281,14 +1274,9 @@ contains p % delayed_group = 0 ! sample from prompt neutron energy distribution - law = rxn % edist % law n_sample = 0 do - if (law == 44 .or. law == 61) then - call sample_energy(rxn%edist, p % E, E_out, prob) - else - call sample_energy(rxn%edist, p % E, E_out) - end if + call rxn%secondary%sample(p%E, E_out, prob) ! resample if energy is greater than maximum neutron energy if (E_out < energy_max_neutron) exit @@ -1317,41 +1305,18 @@ contains type(Particle), intent(inout) :: p integer :: i ! loop index - integer :: law ! secondary energy distribution law real(8) :: E ! energy in lab (incoming/outgoing) real(8) :: mu ! cosine of scattering angle in lab real(8) :: A ! atomic weight ratio of nuclide real(8) :: E_in ! incoming energy real(8) :: E_cm ! outgoing energy in center-of-mass - real(8) :: Q ! Q-value of reaction real(8) :: yield ! neutron yield ! copy energy of neutron E_in = p % E - ! determine A and Q - A = nuc % awr - Q = rxn % Q_value - - ! determine secondary energy distribution law - law = rxn % edist % law - - ! sample scattering angle - mu = sample_angle(rxn, E_in) - ! sample outgoing energy - if (law == 44 .or. law == 61) then - call sample_energy(rxn%edist, E_in, E, mu) - ! Because of floating-point roundoff, it may be possible for mu to be - ! outside of the range [-1,1). In these cases, we just set mu to exactly - ! -1 or 1 - - if (abs(mu) > ONE) mu = sign(ONE,mu) - elseif (law == 66) then - call sample_energy(rxn%edist, E_in, E, A=A, Q=Q) - else - call sample_energy(rxn%edist, E_in, E) - end if + call rxn%secondary%sample(E_in, E, mu) ! if scattering system is in center-of-mass, transfer cosine of scattering ! angle and outgoing energy from CM to LAB @@ -1359,20 +1324,20 @@ contains E_cm = E ! determine outgoing energy in lab + A = nuc%awr E = E_cm + (E_in + TWO * mu * (A+ONE) * sqrt(E_in * E_cm)) & / ((A+ONE)*(A+ONE)) ! determine outgoing angle in lab mu = mu * sqrt(E_cm/E) + ONE/(A+ONE) * sqrt(E_in/E) - - ! Because of floating-point roundoff, it may be possible for mu to be - ! outside of the range [-1,1). In these cases, we just set mu to exactly - ! -1 or 1 - - if (abs(mu) > ONE) mu = sign(ONE,mu) end if - ! Set outgoing energy and scattering angle + ! Because of floating-point roundoff, it may be possible for mu to be + ! outside of the range [-1,1). In these cases, we just set mu to exactly -1 + ! or 1 + if (abs(mu) > ONE) mu = sign(ONE,mu) + + ! Set outgoing energy and scattering angle p % E = E p % mu = mu @@ -1391,893 +1356,4 @@ contains end subroutine inelastic_scatter -!=============================================================================== -! SAMPLE_ANGLE samples the cosine of the angle between incident and exiting -! particle directions either from 32 equiprobable bins or from a tabular -! distribution. -!=============================================================================== - - function sample_angle(rxn, E) result(mu) - type(Reaction), intent(in) :: rxn ! reaction - real(8), intent(in) :: E ! incoming energy - - real(8) :: xi ! random number on [0,1) - integer :: interp ! type of interpolation - integer :: type ! angular distribution type - integer :: i ! incoming energy bin - integer :: n ! number of incoming energy bins - integer :: lc ! location in data array - integer :: NP ! number of points in cos distribution - integer :: k ! index on cosine grid - real(8) :: r ! interpolation factor on incoming energy - real(8) :: frac ! interpolation fraction on cosine - real(8) :: mu0 ! cosine in bin k - real(8) :: mu1 ! cosine in bin k+1 - real(8) :: mu ! final cosine sampled - real(8) :: c_k ! cumulative frequency at k - real(8) :: c_k1 ! cumulative frequency at k+1 - real(8) :: p0,p1 ! probability distribution - - ! check if reaction has angular distribution -- if not, sample outgoing - ! angle isotropically - if (.not. rxn % has_angle_dist) then - mu = TWO * prn() - ONE - return - end if - - ! determine number of incoming energies - n = rxn % adist % n_energy - - ! find energy bin and calculate interpolation factor -- if the energy is - ! outside the range of the tabulated energies, choose the first or last bins - if (E < rxn % adist % energy(1)) then - i = 1 - r = ZERO - elseif (E > rxn % adist % energy(n)) then - i = n - 1 - r = ONE - else - i = binary_search(rxn % adist % energy, n, E) - r = (E - rxn % adist % energy(i)) / & - (rxn % adist % energy(i+1) - rxn % adist % energy(i)) - end if - - ! Sample between the ith and (i+1)th bin - if (r > prn()) i = i + 1 - - ! check whether this is a 32-equiprobable bin or a tabular distribution - lc = rxn % adist % location(i) - type = rxn % adist % type(i) - if (type == ANGLE_ISOTROPIC) then - mu = TWO * prn() - ONE - elseif (type == ANGLE_32_EQUI) then - ! sample cosine bin - xi = prn() - k = 1 + int(32.0_8*xi) - - ! calculate cosine - mu0 = rxn % adist % data(lc + k) - mu1 = rxn % adist % data(lc + k+1) - mu = mu0 + (32.0_8 * xi - k + ONE) * (mu1 - mu0) - - elseif (type == ANGLE_TABULAR) then - interp = int(rxn % adist % data(lc + 1)) - NP = int(rxn % adist % data(lc + 2)) - - ! determine outgoing cosine bin - xi = prn() - lc = lc + 2 - c_k = rxn % adist % data(lc + 2*NP + 1) - do k = 1, NP - 1 - c_k1 = rxn % adist % data(lc + 2*NP + k+1) - if (xi < c_k1) exit - c_k = c_k1 - end do - - ! check to make sure k is <= NP - 1 - k = min(k, NP - 1) - - p0 = rxn % adist % data(lc + NP + k) - mu0 = rxn % adist % data(lc + k) - if (interp == HISTOGRAM) then - ! Histogram interpolation - if (p0 > ZERO) then - mu = mu0 + (xi - c_k)/p0 - else - mu = mu0 - end if - - elseif (interp == LINEAR_LINEAR) then - ! Linear-linear interpolation - p1 = rxn % adist % data(lc + NP + k+1) - mu1 = rxn % adist % data(lc + k+1) - - frac = (p1 - p0)/(mu1 - mu0) - if (frac == ZERO) then - mu = mu0 + (xi - c_k)/p0 - else - mu = mu0 + (sqrt(max(ZERO, p0*p0 + 2*frac*(xi - c_k))) - p0)/frac - end if - else - ! call write_particle_restart(p) - call fatal_error("Unknown interpolation type: " // trim(to_str(interp))) - end if - - ! Because of floating-point roundoff, it may be possible for mu to be - ! outside of the range [-1,1). In these cases, we just set mu to exactly - ! -1 or 1 - - if (abs(mu) > ONE) mu = sign(ONE,mu) - - else - ! call write_particle_restart(p) - call fatal_error("Unknown angular distribution type: " & - &// trim(to_str(type))) - end if - - end function sample_angle - -!=============================================================================== -! SAMPLE_ENERGY samples an outgoing energy distribution, either for a secondary -! neutron from a collision or for a prompt/delayed fission neutron -!=============================================================================== - - recursive subroutine sample_energy(edist, E_in, E_out, mu_out, A, Q) - type(DistEnergy), intent(in) :: edist - real(8), intent(in) :: E_in ! incoming energy of neutron - real(8), intent(out) :: E_out ! outgoing energy - real(8), intent(inout), optional :: mu_out ! outgoing cosine of angle - real(8), intent(in), optional :: A ! mass number of nuclide - real(8), intent(in), optional :: Q ! Q-value of reaction - - integer :: i ! index on incoming energy grid - integer :: k ! sampled index on outgoing grid - integer :: l ! sampled index on incoming grid - integer :: n_sample ! number of rejections - integer :: lc ! dummy index - integer :: NR ! number of interpolation regions - integer :: NE ! number of energies - integer :: NET ! number of outgoing energies - integer :: INTTp ! combination of INTT and ND - integer :: INTT ! 1 = histogram, 2 = linear-linear - integer :: JJ ! 1 = histogram, 2 = linear-linear - integer :: ND ! number of discrete lines - integer :: NP ! number of points in distribution - - real(8) :: p_valid ! probability of law validity - - real(8) :: E_i_1, E_i_K ! endpoints on outgoing grid i - real(8) :: E_i1_1, E_i1_K ! endpoints on outgoing grid i+1 - real(8) :: E_1, E_K ! endpoints interpolated between i and i+1 - - real(8) :: E_l_k, E_l_k1 ! adjacent E on outgoing grid l - real(8) :: p_l_k, p_l_k1 ! adjacent p on outgoing grid l - real(8) :: c_k, c_k1 ! cumulative probability - - real(8) :: KM_A ! Kalbach-Mann parameter R - real(8) :: KM_R ! Kalbach-Mann parameter R - real(8) :: A_k, A_k1 ! Kalbach-Mann A on outgoing grid l - real(8) :: R_k, R_k1 ! Kalbach-Mann R on outgoing grid l - - real(8) :: Watt_a, Watt_b ! Watt spectrum parameters - - real(8) :: mu_k ! angular cosine in bin k - real(8) :: mu_k1 ! angular cosine in bin k+1 - real(8) :: p_k ! angular pdf in bin k - real(8) :: p_k1 ! angular pdf in bin k+1 - - real(8) :: r ! interpolation factor on incoming energy - real(8) :: frac ! interpolation factor on outgoing energy - real(8) :: U ! restriction energy - real(8) :: T ! nuclear temperature - - real(8) :: Ap ! total mass ratio for n-body dist - integer :: n_bodies ! number of bodies for n-body dist - real(8) :: E_max ! parameter for n-body dist - real(8) :: x, y, v ! intermediate variables for n-body dist - real(8) :: r1, r2, r3, r4, r5, r6 - logical :: histogram_interp ! use histogram interpolation on incoming energy - - ! ========================================================================== - ! SAMPLE ENERGY DISTRIBUTION IF THERE ARE MULTIPLE - - if (associated(edist % next)) then - p_valid = interpolate_tab1(edist % p_valid, E_in) - - if (prn() > p_valid) then - if (edist % law == 44 .or. edist % law == 61) then - call sample_energy(edist%next, E_in, E_out, mu_out) - elseif (edist % law == 66) then - call sample_energy(edist%next, E_in, E_out, A=A, Q=Q) - else - call sample_energy(edist%next, E_in, E_out) - end if - return - end if - end if - - ! Determine which secondary energy distribution law to use - select case (edist % law) - case (1) - ! ======================================================================= - ! TABULAR EQUIPROBABLE ENERGY BINS - - ! read number of interpolation regions, incoming energies, and outgoing - ! energies - NR = int(edist % data(1)) - NE = int(edist % data(2 + 2*NR)) - NET = int(edist % data(3 + 2*NR + NE)) - if (NR > 0) then - ! call write_particle_restart(p) - call fatal_error("Multiple interpolation regions not supported while & - &attempting to sample equiprobable energy bins.") - end if - - ! determine index on incoming energy grid and interpolation factor - lc = 2 + 2*NR - i = binary_search(edist % data(lc+1:lc+NE), NE, E_in) - r = (E_in - edist%data(lc+i)) / & - (edist%data(lc+i+1) - edist%data(lc+i)) - - ! Sample outgoing energy bin - r1 = prn() - k = 1 + int(NET * r1) - - ! Determine E_1 and E_K - lc = 3 + 3*NR + NE + (i-1)*NET - E_i_1 = edist % data(lc + 1) - E_i_K = edist % data(lc + NET) - - lc = 3 + 3*NR + NE + i*NET - E_i1_1 = edist % data(lc + 1) - E_i1_K = edist % data(lc + NET) - - E_1 = E_i_1 + r*(E_i1_1 - E_i_1) - E_K = E_i_K + r*(E_i1_K - E_i_K) - - ! Randomly select between the outgoing table for incoming energy E_i and - ! E_(i+1) - if (prn() < r) then - l = i + 1 - else - l = i - end if - - ! Determine E_l_k and E_l_k+1 - lc = 3 + 2*NR + NE + (l-1)*NET - E_l_k = edist % data(lc+k) - E_l_k1 = edist % data(lc+k+1) - - ! Determine E' (denoted here as E_out) - r2 = prn() - E_out = E_l_k + r2*(E_l_k1 - E_l_k) - - ! Now interpolate between incident energy bins i and i + 1 - if (l == i) then - E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1) - else - E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1) - end if - - case (3) - ! ======================================================================= - ! INELASTIC LEVEL SCATTERING - - E_out = edist%data(2) * (E_in - edist%data(1)) - - case (4) - ! ======================================================================= - ! CONTINUOUS TABULAR DISTRIBUTION - - ! read number of interpolation regions and incoming energies - NR = int(edist % data(1)) - NE = int(edist % data(2 + 2*NR)) - if (NR == 1) then - histogram_interp = (edist % data(3) == 1) - else if (NR > 1) then - ! call write_particle_restart(p) - call fatal_error("Multiple interpolation regions not supported while & - &attempting to sample continuous tabular distribution.") - else - histogram_interp = .false. - end if - - ! find energy bin and calculate interpolation factor -- if the energy is - ! outside the range of the tabulated energies, choose the first or last - ! bins - lc = 2 + 2*NR - if (E_in < edist % data(lc+1)) then - i = 1 - r = ZERO - elseif (E_in > edist % data(lc+NE)) then - i = NE - 1 - r = ONE - else - i = binary_search(edist % data(lc+1:lc+NE), NE, E_in) - r = (E_in - edist%data(lc+i)) / & - (edist%data(lc+i+1) - edist%data(lc+i)) - end if - - ! Sample between the ith and (i+1)th bin - if (histogram_interp) then - l = i - else - r2 = prn() - if (r > r2) then - l = i + 1 - else - l = i - end if - end if - - ! interpolation for energy E1 and EK - lc = int(edist%data(2 + 2*NR + NE + i)) - NP = int(edist%data(lc + 2)) - E_i_1 = edist%data(lc + 2 + 1) - E_i_K = edist%data(lc + 2 + NP) - - lc = int(edist%data(2 + 2*NR + NE + i + 1)) - NP = int(edist%data(lc + 2)) - E_i1_1 = edist%data(lc + 2 + 1) - E_i1_K = edist%data(lc + 2 + NP) - - E_1 = E_i_1 + r*(E_i1_1 - E_i_1) - E_K = E_i_K + r*(E_i1_K - E_i_K) - - ! determine location of outgoing energies, pdf, cdf for E(l) - lc = int(edist % data(2 + 2*NR + NE + l)) - - ! determine type of interpolation and number of discrete lines - INTTp = int(edist % data(lc + 1)) - NP = int(edist % data(lc + 2)) - if (INTTp > 10) then - INTT = mod(INTTp,10) - ND = (INTTp - INTT)/10 - else - INTT = INTTp - ND = 0 - end if - - if (ND > 0) then - ! discrete lines present - ! call write_particle_restart(p) - call fatal_error("Discrete lines in continuous tabular distributed not & - &yet supported") - end if - - ! determine outgoing energy bin - r1 = prn() - lc = lc + 2 ! start of EOUT - c_k = edist % data(lc + 2*NP + 1) - do k = 1, NP - 1 - c_k1 = edist % data(lc + 2*NP + k+1) - if (r1 < c_k1) exit - c_k = c_k1 - end do - - ! check to make sure k is <= NP - 1 - k = min(k, NP - 1) - - E_l_k = edist % data(lc+k) - p_l_k = edist % data(lc+NP+k) - if (INTT == HISTOGRAM) then - ! Histogram interpolation - if (p_l_k > ZERO) then - E_out = E_l_k + (r1 - c_k)/p_l_k - else - E_out = E_l_k - end if - - elseif (INTT == LINEAR_LINEAR) then - ! Linear-linear interpolation - E_l_k1 = edist % data(lc+k+1) - p_l_k1 = edist % data(lc+NP+k+1) - - frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k) - if (frac == ZERO) then - E_out = E_l_k + (r1 - c_k)/p_l_k - else - E_out = E_l_k + (sqrt(max(ZERO, p_l_k*p_l_k + & - 2*frac*(r1 - c_k))) - p_l_k)/frac - end if - else - ! call write_particle_restart(p) - call fatal_error("Unknown interpolation type: " // trim(to_str(INTT))) - end if - - ! Now interpolate between incident energy bins i and i + 1 - if (.not. histogram_interp) then - if (l == i) then - E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1) - else - E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1) - end if - end if - - case (5) - ! ======================================================================= - ! GENERAL EVAPORATION SPECTRUM - - case (7) - ! ======================================================================= - ! MAXWELL FISSION SPECTRUM - - ! read number of interpolation regions and incoming energies - NR = int(edist % data(1)) - NE = int(edist % data(2 + 2*NR)) - - ! determine nuclear temperature from tabulated function - T = interpolate_tab1(edist % data, E_in) - - ! determine restriction energy - lc = 2 + 2*NR + 2*NE - U = edist % data(lc + 1) - - n_sample = 0 - do - ! sample maxwell fission spectrum - E_out = maxwell_spectrum(T) - - ! accept energy based on restriction energy - if (E_out <= E_in - U) exit - - ! check for large number of rejections - n_sample = n_sample + 1 - if (n_sample == MAX_SAMPLE) then - ! call write_particle_restart(p) - call fatal_error("Too many rejections on Maxwell fission spectrum.") - end if - end do - - case (9) - ! ======================================================================= - ! EVAPORATION SPECTRUM - - ! read number of interpolation regions and incoming energies - NR = int(edist % data(1)) - NE = int(edist % data(2 + 2*NR)) - - ! determine nuclear temperature from tabulated function - T = interpolate_tab1(edist % data, E_in) - - ! determine restriction energy - lc = 2 + 2*NR + 2*NE - U = edist % data(lc + 1) - - y = (E_in - U)/T - v = 1 - exp(-y) - - ! sample outgoing energy based on evaporation spectrum probability - ! density function - n_sample = 0 - do - x = -log((1 - v*prn())*(1 - v*prn())) - if (x <= y) exit - - ! check for large number of rejections - n_sample = n_sample + 1 - if (n_sample == MAX_SAMPLE) then - ! call write_particle_restart(p) - call fatal_error("Too many rejections on evaporation spectrum.") - end if - end do - - E_out = x*T - - case (11) - ! ======================================================================= - ! ENERGY-DEPENDENT WATT SPECTRUM - - ! read number of interpolation regions and incoming energies for - ! parameter 'a' - NR = int(edist % data(1)) - NE = int(edist % data(2 + 2*NR)) - - ! determine Watt parameter 'a' from tabulated function - Watt_a = interpolate_tab1(edist % data, E_in) - - ! determine Watt parameter 'b' from tabulated function - lc = 2 + 2*(NR + NE) - Watt_b = interpolate_tab1(edist % data, E_in, lc + 1) - - ! read number of interpolation regions and incoming energies for - ! parameter 'a' - NR = int(edist % data(lc + 1)) - NE = int(edist % data(lc + 2 + 2*NR)) - - ! determine restriction energy - lc = lc + 2 + 2*(NR + NE) - U = edist % data(lc + 1) - - n_sample = 0 - do - ! Sample energy-dependent Watt fission spectrum - E_out = watt_spectrum(Watt_a, Watt_b) - - ! accept energy based on restriction energy - if (E_out <= E_in - U) exit - - ! check for large number of rejections - n_sample = n_sample + 1 - if (n_sample == MAX_SAMPLE) then - ! call write_particle_restart(p) - call fatal_error("Too many rejections on Watt spectrum.") - end if - end do - - case (44) - ! ======================================================================= - ! KALBACH-MANN CORRELATED SCATTERING - - if (.not. present(mu_out)) then - ! call write_particle_restart(p) - call fatal_error("Law 44 called without giving mu_out as argument.") - end if - - ! read number of interpolation regions and incoming energies - NR = int(edist % data(1)) - NE = int(edist % data(2 + 2*NR)) - if (NR > 0) then - ! call write_particle_restart(p) - call fatal_error("Multiple interpolation regions not supported while & - &attempting to sample Kalbach-Mann distribution.") - end if - - ! find energy bin and calculate interpolation factor -- if the energy is - ! outside the range of the tabulated energies, choose the first or last - ! bins - lc = 2 + 2*NR - if (E_in < edist % data(lc+1)) then - i = 1 - r = ZERO - elseif (E_in > edist % data(lc+NE)) then - i = NE - 1 - r = ONE - else - i = binary_search(edist % data(lc+1:lc+NE), NE, E_in) - r = (E_in - edist%data(lc+i)) / & - (edist%data(lc+i+1) - edist%data(lc+i)) - end if - - ! Sample between the ith and (i+1)th bin - r2 = prn() - if (r > r2) then - l = i + 1 - else - l = i - end if - - ! determine endpoints on grid i - lc = int(edist%data(2+2*NR+NE + i)) ! start of LDAT for i - NP = int(edist%data(lc + 2)) - E_i_1 = edist%data(lc + 2 + 1) - E_i_K = edist%data(lc + 2 + NP) - - ! determine endpoints on grid i+1 - lc = int(edist%data(2+2*NR+NE + i+1)) ! start of LDAT for i+1 - NP = int(edist%data(lc + 2)) - E_i1_1 = edist%data(lc + 2 + 1) - E_i1_K = edist%data(lc + 2 + NP) - - E_1 = E_i_1 + r*(E_i1_1 - E_i_1) - E_K = E_i_K + r*(E_i1_K - E_i_K) - - ! determine location of outgoing energies, pdf, cdf for E(l) - lc = int(edist % data(2 + 2*NR + NE + l)) - - ! determine type of interpolation and number of discrete lines - INTTp = int(edist % data(lc + 1)) - NP = int(edist % data(lc + 2)) - if (INTTp > 10) then - INTT = mod(INTTp,10) - ND = (INTTp - INTT)/10 - else - INTT = INTTp - ND = 0 - end if - - if (ND > 0) then - ! discrete lines present - ! call write_particle_restart(p) - call fatal_error("Discrete lines in continuous tabular distributed not & - &yet supported") - end if - - ! determine outgoing energy bin - r1 = prn() - lc = lc + 2 ! start of EOUT - c_k = edist % data(lc + 2*NP + 1) - do k = 1, NP - 1 - c_k1 = edist % data(lc + 2*NP + k+1) - if (r1 < c_k1) exit - c_k = c_k1 - end do - - ! check to make sure k is <= NP - 1 - k = min(k, NP - 1) - - E_l_k = edist % data(lc+k) - p_l_k = edist % data(lc+NP+k) - if (INTT == HISTOGRAM) then - ! Histogram interpolation - if (p_l_k > ZERO) then - E_out = E_l_k + (r1 - c_k)/p_l_k - else - E_out = E_l_k - end if - - ! Determine Kalbach-Mann parameters - KM_R = edist % data(lc + 3*NP + k) - KM_A = edist % data(lc + 4*NP + k) - - elseif (INTT == LINEAR_LINEAR) then - ! Linear-linear interpolation - E_l_k1 = edist % data(lc+k+1) - p_l_k1 = edist % data(lc+NP+k+1) - - ! Find E prime - frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k) - if (frac == ZERO) then - E_out = E_l_k + (r1 - c_k)/p_l_k - else - E_out = E_l_k + (sqrt(max(ZERO, p_l_k*p_l_k + & - 2*frac*(r1 - c_k))) - p_l_k)/frac - end if - - ! Determine Kalbach-Mann parameters - R_k = edist % data(lc + 3*NP + k) - R_k1 = edist % data(lc + 3*NP + k+1) - A_k = edist % data(lc + 4*NP + k) - A_k1 = edist % data(lc + 4*NP + k+1) - - KM_R = R_k + (R_k1 - R_k)*(E_out - E_l_k)/(E_l_k1 - E_l_k) - KM_A = A_k + (A_k1 - A_k)*(E_out - E_l_k)/(E_l_k1 - E_l_k) - else - ! call write_particle_restart() - call fatal_error("Unknown interpolation type: " // trim(to_str(INTT))) - end if - - ! Now interpolate between incident energy bins i and i + 1 - if (l == i) then - E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1) - else - E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1) - end if - - ! Sampled correlated angle from Kalbach-Mann parameters - r3 = prn() - r4 = prn() - if (r3 > KM_R) then - T = (TWO*r4 - ONE) * sinh(KM_A) - mu_out = log(T + sqrt(T*T + ONE))/KM_A - else - mu_out = log(r4*exp(KM_A) + (ONE - r4)*exp(-KM_A))/KM_A - end if - - case (61) - ! ======================================================================= - ! CORRELATED ENERGY AND ANGLE DISTRIBUTION - - if (.not. present(mu_out)) then - ! call write_particle_restart() - call fatal_error("Law 61 called without giving mu_out as argument.") - end if - - ! read number of interpolation regions and incoming energies - NR = int(edist % data(1)) - NE = int(edist % data(2 + 2*NR)) - if (NR > 0) then - ! call write_particle_restart() - call fatal_error("Multiple interpolation regions not supported while & - &attempting to sample correlated energy-angle distribution.") - end if - - ! find energy bin and calculate interpolation factor -- if the energy is - ! outside the range of the tabulated energies, choose the first or last - ! bins - lc = 2 + 2*NR - if (E_in < edist % data(lc+1)) then - i = 1 - r = ZERO - elseif (E_in > edist % data(lc+NE)) then - i = NE - 1 - r = ONE - else - i = binary_search(edist % data(lc+1:lc+NE), NE, E_in) - r = (E_in - edist%data(lc+i)) / & - (edist%data(lc+i+1) - edist%data(lc+i)) - end if - - ! Sample between the ith and (i+1)th bin - r2 = prn() - if (r > r2) then - l = i + 1 - else - l = i - end if - - ! determine endpoints on grid i - lc = int(edist%data(2+2*NR+NE + i)) ! start of LDAT for i - NP = int(edist%data(lc + 2)) - E_i_1 = edist%data(lc + 2 + 1) - E_i_K = edist%data(lc + 2 + NP) - - ! determine endpoints on grid i+1 - lc = int(edist%data(2+2*NR+NE + i+1)) ! start of LDAT for i+1 - NP = int(edist%data(lc + 2)) - E_i1_1 = edist%data(lc + 2 + 1) - E_i1_K = edist%data(lc + 2 + NP) - - E_1 = E_i_1 + r*(E_i1_1 - E_i_1) - E_K = E_i_K + r*(E_i1_K - E_i_K) - - ! determine location of outgoing energies, pdf, cdf for E(l) - lc = int(edist % data(2 + 2*NR + NE + l)) - - ! determine type of interpolation and number of discrete lines - INTTp = int(edist % data(lc + 1)) - NP = int(edist % data(lc + 2)) - if (INTTp > 10) then - INTT = mod(INTTp,10) - ND = (INTTp - INTT)/10 - else - INTT = INTTp - ND = 0 - end if - - if (ND > 0) then - ! discrete lines present - ! call write_particle_restart() - call fatal_error("Discrete lines in continuous tabular distributed not & - &yet supported") - end if - - ! determine outgoing energy bin - r1 = prn() - lc = lc + 2 ! start of EOUT - c_k = edist % data(lc + 2*NP + 1) - do k = 1, NP - 1 - c_k1 = edist % data(lc + 2*NP + k+1) - if (r1 < c_k1) exit - c_k = c_k1 - end do - - ! check to make sure k is <= NP - 1 - k = min(k, NP - 1) - - E_l_k = edist % data(lc+k) - p_l_k = edist % data(lc+NP+k) - if (INTT == HISTOGRAM) then - ! Histogram interpolation - if (p_l_k > ZERO) then - E_out = E_l_k + (r1 - c_k)/p_l_k - else - E_out = E_l_k - end if - - elseif (INTT == LINEAR_LINEAR) then - ! Linear-linear interpolation - E_l_k1 = edist % data(lc+k+1) - p_l_k1 = edist % data(lc+NP+k+1) - - ! Find E prime - frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k) - if (frac == ZERO) then - E_out = E_l_k + (r1 - c_k)/p_l_k - else - E_out = E_l_k + (sqrt(max(ZERO, p_l_k*p_l_k + & - 2*frac*(r1 - c_k))) - p_l_k)/frac - end if - else - ! call write_particle_restart() - call fatal_error("Unknown interpolation type: " // trim(to_str(INTT))) - end if - - ! Now interpolate between incident energy bins i and i + 1 - if (l == i) then - E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1) - else - E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1) - end if - - ! Find correlated angular distribution for closest outgoing energy bin - if (r1 - c_k < c_k1 - r1) then - lc = int(edist % data(lc + 3*NP + k)) - else - lc = int(edist % data(lc + 3*NP + k + 1)) - end if - - ! Check if angular distribution is isotropic - if (lc == 0) then - mu_out = TWO * prn() - ONE - return - end if - - ! interpolation type and number of points in angular distribution - JJ = int(edist % data(lc + 1)) - NP = int(edist % data(lc + 2)) - - ! determine outgoing cosine bin - r3 = prn() - lc = lc + 2 - c_k = edist % data(lc + 2*NP + 1) - do k = 1, NP - 1 - c_k1 = edist % data(lc + 2*NP + k+1) - if (r3 < c_k1) exit - c_k = c_k1 - end do - - ! check to make sure k is <= NP - 1 - k = min(k, NP - 1) - - p_k = edist % data(lc + NP + k) - mu_k = edist % data(lc + k) - if (JJ == HISTOGRAM) then - ! Histogram interpolation - if (p_k > ZERO) then - mu_out = mu_k + (r3 - c_k)/p_k - else - mu_out = mu_k - end if - - elseif (JJ == LINEAR_LINEAR) then - ! Linear-linear interpolation - p_k1 = edist % data(lc + NP + k+1) - mu_k1 = edist % data(lc + k+1) - - frac = (p_k1 - p_k)/(mu_k1 - mu_k) - if (frac == ZERO) then - mu_out = mu_k + (r3 - c_k)/p_k - else - mu_out = mu_k + (sqrt(p_k*p_k + 2*frac*(r3 - c_k))-p_k)/frac - end if - else - ! call write_particle_restart() - call fatal_error("Unknown interpolation type: " // trim(to_str(JJ))) - end if - - case (66) - ! ======================================================================= - ! N-BODY PHASE SPACE DISTRIBUTION - - ! read number of bodies in phase space and total mass ratio - n_bodies = int(edist % data(1)) - Ap = edist % data(2) - - ! determine E_max parameter - E_max = (Ap - ONE)/Ap * (A/(A+ONE)*E_in + Q) - - ! x is essentially a Maxwellian distribution - x = maxwell_spectrum(ONE) - - select case (n_bodies) - case (3) - y = maxwell_spectrum(ONE) - case (4) - r1 = prn() - r2 = prn() - r3 = prn() - y = -log(r1*r2*r3) - case (5) - r1 = prn() - r2 = prn() - r3 = prn() - r4 = prn() - r5 = prn() - r6 = prn() - y = -log(r1*r2*r3*r4) - log(r5) * cos(PI/TWO*r6)**2 - end select - - ! now determine v and E_out - v = x/(x+y) - E_out = E_max * v - - case (67) - ! ======================================================================= - ! LABORATORY ENERGY-ANGLE LAW - - end select - - end subroutine sample_energy - end module physics diff --git a/src/search.F90 b/src/search.F90 index 0c345471c2..f338105f4d 100644 --- a/src/search.F90 +++ b/src/search.F90 @@ -1,7 +1,6 @@ module search use constants - use error, only: fatal_error implicit none diff --git a/src/secondary_correlated.F90 b/src/secondary_correlated.F90 index c3c82b588d..39ce5a1a8d 100644 --- a/src/secondary_correlated.F90 +++ b/src/secondary_correlated.F90 @@ -1,18 +1,21 @@ module secondary_correlated use constants, only: ZERO, ONE, TWO, HISTOGRAM, LINEAR_LINEAR - use distribution_univariate, only: Tabular - use error, only: fatal_error - use secondary_header, only: SecondaryDistribution + use distribution_univariate, only: DistributionContainer + use secondary_header, only: AngleEnergy use random_lcg, only: prn use search, only: binary_search type AngleEnergyTable - type(Tabular) :: energy - type(Tabular), allocatable :: angle(:) + integer :: interpolation + integer :: n_discrete + real(8), allocatable :: e_out(:) + real(8), allocatable :: p(:) + real(8), allocatable :: c(:) + type(DistributionContainer), allocatable :: angle(:) end type AngleEnergyTable - type, extends(SecondaryDistribution) :: CorrelatedAngleEnergy + type, extends(AngleEnergy) :: CorrelatedAngleEnergy integer :: n_region integer, allocatable :: breakpoints(:) integer, allocatable :: interpolation(:) @@ -43,12 +46,6 @@ contains real(8) :: p_l_k, p_l_k1 ! adjacent p on outgoing grid l real(8) :: c_k, c_k1 ! cumulative probability - ! TODO: Write error during initialization - if (this%n_region > 1) then - call fatal_error("Multiple interpolation regions not supported while & - &attempting to sample Kalbach-Mann distribution.") - end if - ! find energy bin and calculate interpolation factor -- if the energy is ! outside the range of the tabulated energies, choose the first or last bins n_energy_in = size(this%energy_in) @@ -72,30 +69,23 @@ contains end if ! interpolation for energy E1 and EK - n_energy_out = size(this%table(i)%energy%x) - E_i_1 = this%table(i)%energy%x(1) - E_i_K = this%table(i)%energy%x(n_energy_out) + n_energy_out = size(this%table(i)%e_out) + E_i_1 = this%table(i)%e_out(1) + E_i_K = this%table(i)%e_out(n_energy_out) - n_energy_out = size(this%table(i+1)%energy%x) - E_i1_1 = this%table(i+1)%energy%x(1) - E_i1_K = this%table(i+1)%energy%x(n_energy_out) + n_energy_out = size(this%table(i+1)%e_out) + E_i1_1 = this%table(i+1)%e_out(1) + E_i1_K = this%table(i+1)%e_out(n_energy_out) E_1 = E_i_1 + r*(E_i1_1 - E_i_1) E_K = E_i_K + r*(E_i1_K - E_i_K) -!!$ ! TODO: Write error at initizliation -!!$ if (this%table(l)%n_discrete > 0) then -!!$ ! discrete lines present -!!$ call fatal_error("Discrete lines in continuous tabular distributed not & -!!$ &yet supported") -!!$ end if - ! determine outgoing energy bin - n_energy_out = size(this%table(l)%energy%x) + n_energy_out = size(this%table(l)%e_out) r1 = prn() - c_k = this%table(l)%energy%c(1) + c_k = this%table(l)%c(1) do k = 1, n_energy_out - 1 - c_k1 = this%table(l)%energy%c(k+1) + c_k1 = this%table(l)%c(k+1) if (r1 < c_k1) exit c_k = c_k1 end do @@ -103,9 +93,9 @@ contains ! check to make sure k is <= NP - 1 k = min(k, n_energy_out - 1) - E_l_k = this%table(l)%energy%x(k) - p_l_k = this%table(l)%energy%p(k) - if (this%table(l)%energy%interpolation == HISTOGRAM) then + E_l_k = this%table(l)%e_out(k) + p_l_k = this%table(l)%p(k) + if (this%table(l)%interpolation == HISTOGRAM) then ! Histogram interpolation if (p_l_k > ZERO) then E_out = E_l_k + (r1 - c_k)/p_l_k @@ -113,10 +103,10 @@ contains E_out = E_l_k end if - elseif (this%table(l)%energy%interpolation == LINEAR_LINEAR) then + elseif (this%table(l)%interpolation == LINEAR_LINEAR) then ! Linear-linear interpolation - E_l_k1 = this%table(l)%energy%x(k+1) - p_l_k1 = this%table(l)%energy%p(k+1) + E_l_k1 = this%table(l)%e_out(k+1) + p_l_k1 = this%table(l)%p(k+1) frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k) if (frac == ZERO) then @@ -136,9 +126,9 @@ contains ! Find correlated angular distribution for closest outgoing energy bin if (r1 - c_k < c_k1 - r1) then - mu = this%table(l)%angle(k)%sample() + mu = this%table(l)%angle(k)%obj%sample() else - mu = this%table(l)%angle(k + 1)%sample() + mu = this%table(l)%angle(k + 1)%obj%sample() end if end subroutine correlated_sample diff --git a/src/secondary_header.F90 b/src/secondary_header.F90 index 1f8ad3aed1..30cdab8b88 100644 --- a/src/secondary_header.F90 +++ b/src/secondary_header.F90 @@ -13,13 +13,6 @@ module secondary_header class(AngleEnergy), allocatable :: obj end type AngleEnergyContainer - type :: SecondaryDistribution - type(Tab1), allocatable :: applicability(:) - type(AngleEnergyContainer), allocatable :: distribution(:) - contains - procedure :: sample => secondary_sample - end type SecondaryDistribution - abstract interface subroutine iSampleAngleEnergy(this, E_in, E_out, mu) import AngleEnergy @@ -30,6 +23,18 @@ module secondary_header end subroutine iSampleAngleEnergy end interface +!=============================================================================== +! SECONDARYDISTRIBUTION stores a secondary distribution for angle and energy, +! whether correlated or uncorrelated. +!=============================================================================== + + type :: SecondaryDistribution + type(Tab1), allocatable :: applicability(:) + type(AngleEnergyContainer), allocatable :: distribution(:) + contains + procedure :: sample => secondary_sample + end type SecondaryDistribution + contains subroutine secondary_sample(this, E_in, E_out, mu) diff --git a/src/secondary_kalbach.F90 b/src/secondary_kalbach.F90 index d0d2f32bf2..9764502ef6 100644 --- a/src/secondary_kalbach.F90 +++ b/src/secondary_kalbach.F90 @@ -1,8 +1,7 @@ module secondary_kalbach use constants, only: ZERO, ONE, TWO, HISTOGRAM, LINEAR_LINEAR - use error, only: fatal_error - use secondary_header, only: SecondaryDistribution + use secondary_header, only: AngleEnergy use random_lcg, only: prn use search, only: binary_search @@ -16,7 +15,7 @@ module secondary_kalbach real(8), allocatable :: a(:) end type KalbachMannTable - type, extends(SecondaryDistribution) :: KalbachMann + type, extends(AngleEnergy) :: KalbachMann integer :: n_region integer, allocatable :: breakpoints(:) integer, allocatable :: interpolation(:) @@ -49,12 +48,6 @@ contains real(8) :: km_r, km_a ! Kalbach-Mann parameters real(8) :: T - ! TODO: Write error during initialization - if (this%n_region > 1) then - call fatal_error("Multiple interpolation regions not supported while & - &attempting to sample Kalbach-Mann distribution.") - end if - ! find energy bin and calculate interpolation factor -- if the energy is ! outside the range of the tabulated energies, choose the first or last bins n_energy_in = size(this%energy_in) @@ -89,13 +82,6 @@ contains E_1 = E_i_1 + r*(E_i1_1 - E_i_1) E_K = E_i_K + r*(E_i1_K - E_i_K) - ! TODO: Write error at initizliation - if (this%table(l)%n_discrete > 0) then - ! discrete lines present - call fatal_error("Discrete lines in continuous tabular distributed not & - &yet supported") - end if - ! determine outgoing energy bin n_energy_out = size(this%table(l)%e_out) r1 = prn() diff --git a/src/secondary_uncorrelated.F90 b/src/secondary_uncorrelated.F90 index 7d0046aa4f..bd10512bd1 100644 --- a/src/secondary_uncorrelated.F90 +++ b/src/secondary_uncorrelated.F90 @@ -1,10 +1,12 @@ module secondary_uncorrelated use angle_distribution, only: AngleDistribution + use constants, only: ONE, TWO use energy_distribution, only: EnergyDistribution - use secondary_header, only: SecondaryDistribution + use secondary_header, only: AngleEnergy + use random_lcg, only: prn - type, extends(SecondaryDistribution) :: UncorrelatedAngleEnergy + type, extends(AngleEnergy) :: UncorrelatedAngleEnergy type(AngleDistribution) :: angle class(EnergyDistribution), allocatable :: energy contains @@ -20,7 +22,12 @@ contains real(8), intent(out) :: mu ! Sample cosine of scattering angle - mu = this%angle%sample(E_in) + if (allocated(this%angle%energy)) then + mu = this%angle%sample(E_in) + else + ! no angle distribution given => assume isotropic for all energies + mu = TWO*prn() - ONE + end if ! Sample outgoing energy E_out = this%energy%sample(E_in) From a7ca67b33468b7973e5ac4036dbdb7bebf9f8a9e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 18 Nov 2015 15:10:05 -0600 Subject: [PATCH 03/11] Don't sample between energy distributions if there's only one --- src/secondary_header.F90 | 26 +++++++++++++++++--------- 1 file changed, 17 insertions(+), 9 deletions(-) diff --git a/src/secondary_header.F90 b/src/secondary_header.F90 index 30cdab8b88..e0c414826b 100644 --- a/src/secondary_header.F90 +++ b/src/secondary_header.F90 @@ -43,18 +43,26 @@ contains real(8), intent(out) :: E_out real(8), intent(out) :: mu + integer :: n real(8) :: p_valid - do i = 1, size(this%applicability) - ! Determine probability that i-th energy distribution is sampled - p_valid = interpolate_tab1(this%applicability(i), E_in) + n = size(this%applicability) + if (n > 1) then + do i = 1, n + ! Determine probability that i-th energy distribution is sampled + p_valid = interpolate_tab1(this%applicability(i), E_in) + + ! If i-th distribution is sampled, sample energy from the distribution + if (prn() <= p_valid) then + call this%distribution(i)%obj%sample(E_in, E_out, mu) + exit + end if + end do + else + ! If only one distribution is present, go ahead and sample it + call this%distribution(1)%obj%sample(E_in, E_out, mu) + end if - ! If i-th distribution is sampled, sample energy from the distribution - if (prn() <= p_valid) then - call this%distribution(i)%obj%sample(E_in, E_out, mu) - exit - end if - end do end subroutine secondary_sample end module secondary_header From 306f7016b63b25015a20fed0ed2394ae796a3ef5 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 18 Nov 2015 16:08:58 -0600 Subject: [PATCH 04/11] Use ACE-provided cumulative frequencies for angular distributions --- src/ace.F90 | 26 ++++++++------------------ 1 file changed, 8 insertions(+), 18 deletions(-) diff --git a/src/ace.F90 b/src/ace.F90 index 1407a4b132..6fb37c54ad 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -902,8 +902,6 @@ contains integer :: k ! index over energy distributions integer :: interp integer, allocatable :: LC(:) ! locator - real(8), allocatable :: x(:) - real(8), allocatable :: p(:) ! loop over all reactions with secondary neutrons -- NXS(5) does not include ! elastic scattering @@ -977,13 +975,10 @@ contains ! Get probability density data XSS_index = XSS_index + 2 - allocate(x(NP), p(NP)) - x(:) = get_real(NP) - p(:) = get_real(NP) - - ! initialize distribution - call adist%initialize(x, p, interp) - deallocate(x, p) + allocate(adist%x(NP), adist%p(NP), adist%c(NP)) + adist%x(:) = get_real(NP) + adist%p(:) = get_real(NP) + adist%c(:) = get_real(NP) end select end do deallocate(LC) @@ -1069,8 +1064,6 @@ contains integer :: interp integer, allocatable :: L(:) ! locations of distributions for each Ein integer, allocatable :: LC(:) ! locations of distributions for each Ein - real(8), allocatable :: x(:) - real(8), allocatable :: p(:) XSS_index = LDIS + IDAT - 1 @@ -1355,13 +1348,10 @@ contains ! Get probability density data XSS_index = XSS_index + 2 - allocate(x(NP), p(NP)) - x(:) = get_real(NP) - p(:) = get_real(NP) - - ! initialize distribution - call adist%initialize(x, p, interp) - deallocate(x, p) + allocate(adist%x(NP), adist%p(NP), adist%c(NP)) + adist%x(:) = get_real(NP) + adist%p(:) = get_real(NP) + adist%c(:) = get_real(NP) end select end do deallocate(LC) From 6262e07cbd83bb3eb379543f8729c898339b2961 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 22 Dec 2015 09:09:17 -0600 Subject: [PATCH 05/11] A few patches so that RNG streams still match --- src/ace.F90 | 16 ++++++++++++++++ src/secondary_correlated.F90 | 8 ++++++++ src/secondary_kalbach.F90 | 9 +++++++++ src/secondary_uncorrelated.F90 | 8 +++++++- 4 files changed, 40 insertions(+), 1 deletion(-) diff --git a/src/ace.F90 b/src/ace.F90 index 6fb37c54ad..5b003d4067 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -1035,6 +1035,22 @@ contains call get_energy_dist(secondary%distribution(n)%obj, LAW, & JXS(11), IDAT, nuc%awr, nuc%reactions(i)%Q_value) + ! <<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<< + ! Before the secondary distribution refactor, when the angle/energy + ! distribution was uncorrelated, no angle was actually sampled. With + ! the refactor, an angle is always sampled for an uncorrelated + ! distribution even when no angle distribution exists in the ACE file + ! (isotropic is assumed). To preserve the RNG stream, we explicitly + ! mark fission reactions so that we avoid the angle sampling. + if (any(nuc%reactions(i + 1)%MT == & + [N_FISSION, N_F, N_NF, N_2NF, N_3NF])) then + select type (aedist => secondary%distribution(n)%obj) + type is (UncorrelatedAngleEnergy) + aedist%fission = .true. + end select + end if + ! <<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<< + ! Get locator for next distribution LNW = nint(XSS(JXS(11) + LNW - 1)) end do diff --git a/src/secondary_correlated.F90 b/src/secondary_correlated.F90 index 39ce5a1a8d..5b1d49e28d 100644 --- a/src/secondary_correlated.F90 +++ b/src/secondary_correlated.F90 @@ -46,6 +46,14 @@ contains real(8) :: p_l_k, p_l_k1 ! adjacent p on outgoing grid l real(8) :: c_k, c_k1 ! cumulative probability + ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<< + ! Before the secondary distribution refactor, an isotropic polar cosine was + ! always sampled but then overwritten with the polar cosine sampled from the + ! correlated distribution. To preserve the random number stream, we keep + ! this dummy sampling here but can remove it later (will change answers) + mu = TWO*prn() - ONE + ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<< + ! find energy bin and calculate interpolation factor -- if the energy is ! outside the range of the tabulated energies, choose the first or last bins n_energy_in = size(this%energy_in) diff --git a/src/secondary_kalbach.F90 b/src/secondary_kalbach.F90 index 9764502ef6..d38541edd5 100644 --- a/src/secondary_kalbach.F90 +++ b/src/secondary_kalbach.F90 @@ -48,6 +48,14 @@ contains real(8) :: km_r, km_a ! Kalbach-Mann parameters real(8) :: T + ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<< + ! Before the secondary distribution refactor, an isotropic polar cosine was + ! always sampled but then overwritten with the polar cosine sampled from the + ! correlated distribution. To preserve the random number stream, we keep + ! this dummy sampling here but can remove it later (will change answers) + mu = TWO*prn() - ONE + ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<< + ! find energy bin and calculate interpolation factor -- if the energy is ! outside the range of the tabulated energies, choose the first or last bins n_energy_in = size(this%energy_in) @@ -144,6 +152,7 @@ contains r1 = prn() mu = log(r1*exp(km_a) + (ONE - r1)*exp(-km_a))/km_a end if + end subroutine kalbachmann_sample end module secondary_kalbach diff --git a/src/secondary_uncorrelated.F90 b/src/secondary_uncorrelated.F90 index bd10512bd1..40fa53a382 100644 --- a/src/secondary_uncorrelated.F90 +++ b/src/secondary_uncorrelated.F90 @@ -7,6 +7,7 @@ module secondary_uncorrelated use random_lcg, only: prn type, extends(AngleEnergy) :: UncorrelatedAngleEnergy + logical :: fission = .false. type(AngleDistribution) :: angle class(EnergyDistribution), allocatable :: energy contains @@ -22,7 +23,12 @@ contains real(8), intent(out) :: mu ! Sample cosine of scattering angle - if (allocated(this%angle%energy)) then + if (this%fission) then + ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<< + ! For fission, the angle is not used, so just assign a dummy value + mu = ONE + ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<< + elseif (allocated(this%angle%energy)) then mu = this%angle%sample(E_in) else ! no angle distribution given => assume isotropic for all energies From f4d655605a9c3ccde109d764f758a09962ad8e64 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 4 Jan 2016 14:58:08 -0600 Subject: [PATCH 06/11] Remove SAVE attribute (implicit) in global module. This was needed to remove a weird ICE in gfortran 4.6. --- src/global.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/global.F90 b/src/global.F90 index a5943b6b74..4c243ee41f 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -22,7 +22,6 @@ module global #endif implicit none - save ! ============================================================================ ! GEOMETRY-RELATED VARIABLES From 615d733cc8c1df038e8ac916239648e20ef298c2 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 22 Jan 2016 08:29:42 -0600 Subject: [PATCH 07/11] Add lots of comments --- src/angle_distribution.F90 | 22 +++-- src/energy_distribution.F90 | 153 +++++++++++++++++++++------------ src/physics.F90 | 4 +- src/secondary_correlated.F90 | 27 +++--- src/secondary_header.F90 | 32 ++++--- src/secondary_kalbach.F90 | 28 +++--- src/secondary_uncorrelated.F90 | 12 ++- 7 files changed, 177 insertions(+), 101 deletions(-) diff --git a/src/angle_distribution.F90 b/src/angle_distribution.F90 index 448bafd8f7..282d51b73b 100644 --- a/src/angle_distribution.F90 +++ b/src/angle_distribution.F90 @@ -8,6 +8,14 @@ module angle_distribution implicit none private +!=============================================================================== +! ANGLEDISTRIBUTION represents an angular distribution that is to be used in an +! uncorrelated angle-energy distribution. This occurs whenever the angle +! distrbution is given in File 4 in an ENDF file. The distribution of angles +! depends on the incoming energy of the neutron, so this type stores a +! distribution for each of a set of incoming energies. +!=============================================================================== + type, public :: AngleDistribution real(8), allocatable :: energy(:) type(DistributionContainer), allocatable :: distribution(:) @@ -19,17 +27,17 @@ contains function angle_sample(this, E) result(mu) class(AngleDistribution), intent(in) :: this - real(8), intent(in) :: E - real(8) :: mu + real(8), intent(in) :: E ! incoming energy + real(8) :: mu ! sampled cosine of scattering angle - integer :: i - integer :: n - real(8) :: r + integer :: i ! index on incoming energy grid + integer :: n ! number of incoming energies + real(8) :: r ! interpolation factor on incoming energy grid - ! determine number of incoming energies + ! Determine number of incoming energies n = size(this%energy) - ! find energy bin and calculate interpolation factor -- if the energy is + ! Find energy bin and calculate interpolation factor -- if the energy is ! outside the range of the tabulated energies, choose the first or last bins if (E < this%energy(1)) then i = 1 diff --git a/src/energy_distribution.F90 b/src/energy_distribution.F90 index a8b975879f..1dea1bfece 100644 --- a/src/energy_distribution.F90 +++ b/src/energy_distribution.F90 @@ -9,7 +9,9 @@ module energy_distribution !=============================================================================== ! ENERGYDISTRIBUTION (abstract) defines an energy distribution that is a -! function of the incident energy of a projectile +! function of the incident energy of a projectile. Each derived type must +! implement a sample() function that returns a sampled outgoing energy given an +! incoming energy !=============================================================================== type, abstract :: EnergyDistribution @@ -31,19 +33,31 @@ module energy_distribution end type EnergyDistributionContainer !=============================================================================== -! Derived classes +! Derived classes +!=============================================================================== + +!=============================================================================== +! TABULAREQUIPROBABLE represents an energy distribution with tabular +! equiprobable energy bins as given in ACE law 1. This is an older +! representation that has largely been replaced with ACE laws 4, 44, and 61. !=============================================================================== type, extends(EnergyDistribution) :: TabularEquiprobable - integer :: n_region - integer, allocatable :: breakpoints(:) - integer, allocatable :: interpolation(:) - real(8), allocatable :: energy_in(:) - real(8), allocatable :: energy_out(:,:) + integer :: n_region ! number of interpolation regions + integer, allocatable :: breakpoints(:) ! breakpoints of interpolation regions + integer, allocatable :: interpolation(:) ! interpolation region codes + real(8), allocatable :: energy_in(:) ! incoming energies + real(8), allocatable :: energy_out(:,:) ! table of outgoing energies for + ! each incoming energy contains procedure :: sample => equiprobable_sample end type TabularEquiprobable +!=============================================================================== +! LEVELINELASTIC gives the energy distribution for level inelastic scattering by +! neutrons as in ENDF MT=51--90. +!=============================================================================== + type, extends(EnergyDistribution) :: LevelInelastic real(8) :: threshold real(8) :: mass_ratio @@ -51,6 +65,12 @@ module energy_distribution procedure :: sample => level_inelastic_sample end type LevelInelastic +!=============================================================================== +! CONTINUOUSTABULAR gives an energy distribution represented as a tabular +! distribution with histogram or linear-linear interpolation. This corresponds +! to ACE law 4, which NJOY produces for a number of ENDF energy distributions. +!=============================================================================== + type CTTable integer :: interpolation integer :: n_discrete @@ -69,13 +89,23 @@ module energy_distribution procedure :: sample => continuous_sample end type ContinuousTabular +!=============================================================================== +! MAXWELLENERGY gives the energy distribution of neutrons emitted from a Maxwell +! fission spectrum. This corresponds to ACE law 7 and ENDF File 5, LF=7. +!=============================================================================== + type, extends(EnergyDistribution) :: MaxwellEnergy - type(Tab1) :: theta - real(8) :: u + type(Tab1) :: theta ! incoming-energy-dependent parameter + real(8) :: u ! restriction energy contains procedure :: sample => maxwellenergy_sample end type MaxwellEnergy +!=============================================================================== +! EVAPORATION represents an evaporation spectrum corresponding to ACE law 9 and +! ENDF File 5, LF=9. +!=============================================================================== + type, extends(EnergyDistribution) :: Evaporation type(Tab1) :: theta real(8) :: u @@ -83,6 +113,11 @@ module energy_distribution procedure :: sample => evaporation_sample end type Evaporation +!=============================================================================== +! EVAPORATION gives the energy distribution of neutrons emitted from a Watt +! fission spectrum. This corresponds to ACE law 11 and ENDF File 5, LF=11. +!=============================================================================== + type, extends(EnergyDistribution) :: WattEnergy type(Tab1) :: a type(Tab1) :: b @@ -91,6 +126,12 @@ module energy_distribution procedure :: sample => watt_sample end type WattEnergy +!=============================================================================== +! NBODYPHASESPACE gives the energy distribution for particles emitted from +! neutron and charged-particle reactions. This corresponds to ACE law 66 and +! ENDF File 6, LAW=6. +!=============================================================================== + type, extends(EnergyDistribution) :: NBodyPhaseSpace integer :: n_bodies real(8) :: mass_ratio @@ -104,12 +145,12 @@ contains function equiprobable_sample(this, E_in) result(E_out) class(TabularEquiprobable), intent(in) :: this - real(8), intent(in) :: E_in - real(8) :: E_out + real(8), intent(in) :: E_in ! incoming energy + real(8) :: E_out ! sampled outgoing energy - integer :: i, k, l - integer :: n_energy_in - integer :: n_energy_out + integer :: i, k, l ! indices + integer :: n_energy_in ! number of incoming energies + integer :: n_energy_out ! number of outgoing energies real(8) :: r ! interpolation factor on incoming energy real(8) :: E_i_1, E_i_K ! endpoints on outgoing grid i real(8) :: E_i1_1, E_i1_K ! endpoints on outgoing grid i+1 @@ -120,7 +161,7 @@ contains n_energy_in = size(this%energy_in) n_energy_out = size(this%energy_out, 1) - ! determine index on incoming energy grid and interpolation factor + ! Determine index on incoming energy grid and interpolation factor i = binary_search(this%energy_in, size(this%energy_in), E_in) r = (E_in - this%energy_in(i)) / & (this%energy_in(i+1) - this%energy_in(i)) @@ -171,31 +212,31 @@ contains function continuous_sample(this, E_in) result(E_out) class(ContinuousTabular), intent(in) :: this - real(8), intent(in) :: E_in - real(8) :: E_out + real(8), intent(in) :: E_in ! incoming energy + real(8) :: E_out ! sampled outgoing energy - integer :: i, k, l - integer :: n_energy_in - integer :: n_energy_out - real(8) :: r ! interpolation factor on incoming energy - real(8) :: r1 ! random number on [0,1) - real(8) :: frac ! interpolation factor on outgoing energy - real(8) :: E_i_1, E_i_K ! endpoints on outgoing grid i - real(8) :: E_i1_1, E_i1_K ! endpoints on outgoing grid i+1 - real(8) :: E_1, E_K ! endpoints interpolated between i and i+1 - real(8) :: E_l_k, E_l_k1 ! adjacent E on outgoing grid l - real(8) :: p_l_k, p_l_k1 ! adjacent p on outgoing grid l - real(8) :: c_k, c_k1 ! cumulative probability - logical :: histogram_interp + integer :: i, k, l ! indices + integer :: n_energy_in ! number of incoming energies + integer :: n_energy_out ! number of outgoing energies + real(8) :: r ! interpolation factor on incoming energy + real(8) :: r1 ! random number on [0,1) + real(8) :: frac ! interpolation factor on outgoing energy + real(8) :: E_i_1, E_i_K ! endpoints on outgoing grid i + real(8) :: E_i1_1, E_i1_K ! endpoints on outgoing grid i+1 + real(8) :: E_1, E_K ! endpoints interpolated between i and i+1 + real(8) :: E_l_k, E_l_k1 ! adjacent E on outgoing grid l + real(8) :: p_l_k, p_l_k1 ! adjacent p on outgoing grid l + real(8) :: c_k, c_k1 ! cumulative probability + logical :: histogram_interp ! whether histogram interpolation is used - ! read number of interpolation regions and incoming energies + ! Read number of interpolation regions and incoming energies if (this%n_region == 1) then histogram_interp = (this%interpolation(1) == 1) else histogram_interp = .false. end if - ! find energy bin and calculate interpolation factor -- if the energy is + ! Find energy bin and calculate interpolation factor -- if the energy is ! outside the range of the tabulated energies, choose the first or last bins n_energy_in = size(this%energy_in) if (E_in < this%energy_in(1)) then @@ -221,7 +262,7 @@ contains end if end if - ! interpolation for energy E1 and EK + ! Interpolation for energy E1 and EK n_energy_out = size(this%energy_out(i)%e_out) E_i_1 = this%energy_out(i)%e_out(1) E_i_K = this%energy_out(i)%e_out(n_energy_out) @@ -233,7 +274,7 @@ contains E_1 = E_i_1 + r*(E_i1_1 - E_i_1) E_K = E_i_K + r*(E_i1_K - E_i_K) - ! determine outgoing energy bin + ! Determine outgoing energy bin n_energy_out = size(this%energy_out(l)%e_out) r1 = prn() c_k = this%energy_out(l)%c(1) @@ -243,7 +284,7 @@ contains c_k = c_k1 end do - ! check to make sure k is <= NP - 1 + ! Check to make sure k is <= NP - 1 k = min(k, n_energy_out - 1) E_l_k = this%energy_out(l)%e_out(k) @@ -282,29 +323,29 @@ contains function maxwellenergy_sample(this, E_in) result(E_out) class(MaxwellEnergy), intent(in) :: this - real(8), intent(in) :: E_in - real(8) :: E_out + real(8), intent(in) :: E_in ! incoming energy + real(8) :: E_out ! sampled outgoing energy - real(8) :: theta + real(8) :: theta ! Maxwell distribution parameter ! Get temperature corresponding to incoming energy theta = interpolate_tab1(this%theta, E_in) do - ! sample maxwell fission spectrum + ! Sample maxwell fission spectrum E_out = maxwell_spectrum(theta) - ! accept energy based on restriction energy + ! Accept energy based on restriction energy if (E_out <= E_in - this%u) exit end do end function maxwellenergy_sample function evaporation_sample(this, E_in) result(E_out) class(Evaporation), intent(in) :: this - real(8), intent(in) :: E_in - real(8) :: E_out + real(8), intent(in) :: E_in ! incoming energy + real(8) :: E_out ! sampled outgoing energy - real(8) :: theta + real(8) :: theta ! evaporation spectrum parameter real(8) :: x, y, v ! Get temperature corresponding to incoming energy @@ -313,7 +354,7 @@ contains y = (E_in - this%U)/theta v = 1 - exp(-y) - ! sample outgoing energy based on evaporation spectrum probability + ! Sample outgoing energy based on evaporation spectrum probability ! density function do x = -log((ONE - v*prn())*(ONE - v*prn())) @@ -325,37 +366,37 @@ contains function watt_sample(this, E_in) result(E_out) class(WattEnergy), intent(in) :: this - real(8), intent(in) :: E_in - real(8) :: E_out + real(8), intent(in) :: E_in ! incoming energy + real(8) :: E_out ! sampled outgoing energy - real(8) :: a, b + real(8) :: a, b ! Watt spectrum parameters - ! determine Watt parameter 'a' from tabulated function + ! Determine Watt parameter 'a' from tabulated function a = interpolate_tab1(this%a, E_in) - ! determine Watt parameter 'b' from tabulated function + ! Determine Watt parameter 'b' from tabulated function b = interpolate_tab1(this%b, E_in) do ! Sample energy-dependent Watt fission spectrum E_out = watt_spectrum(a, b) - ! accept energy based on restriction energy + ! Accept energy based on restriction energy if (E_out <= E_in - this%u) exit end do end function watt_sample function nbody_sample(this, E_in) result(E_out) class(NBodyPhaseSpace), intent(in) :: this - real(8), intent(in) :: E_in - real(8) :: E_out + real(8), intent(in) :: E_in ! incoming energy + real(8) :: E_out ! sampled outgoing energy - real(8) :: Ap - real(8) :: E_max + real(8) :: Ap ! total mass of particles in neutron masses + real(8) :: E_max ! maximum possible COM energy real(8) :: x, y, v real(8) :: r1, r2, r3, r4, r5, r6 - ! determine E_max parameter + ! Determine E_max parameter Ap = this%mass_ratio E_max = (Ap - ONE)/Ap * (this%A/(this%A + ONE)*E_in + this%Q) @@ -380,7 +421,7 @@ contains y = -log(r1*r2*r3*r4) - log(r5) * cos(PI/TWO*r6)**2 end select - ! now determine v and E_out + ! Now determine v and E_out v = x/(x+y) E_out = E_max * v end function nbody_sample diff --git a/src/physics.F90 b/src/physics.F90 index 8519f4d2b0..31eb981abb 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -1315,7 +1315,7 @@ contains ! copy energy of neutron E_in = p % E - ! sample outgoing energy + ! sample outgoing energy and scattering cosine call rxn%secondary%sample(E_in, E, mu) ! if scattering system is in center-of-mass, transfer cosine of scattering @@ -1337,7 +1337,7 @@ contains ! or 1 if (abs(mu) > ONE) mu = sign(ONE,mu) - ! Set outgoing energy and scattering angle + ! Set outgoing energy and scattering angle p % E = E p % mu = mu diff --git a/src/secondary_correlated.F90 b/src/secondary_correlated.F90 index 5b1d49e28d..c0289d55eb 100644 --- a/src/secondary_correlated.F90 +++ b/src/secondary_correlated.F90 @@ -6,6 +6,11 @@ module secondary_correlated use random_lcg, only: prn use search, only: binary_search +!=============================================================================== +! CORRELATEDANGLEENERGY represents a correlated angle-energy distribution. This +! corresponds to ACE law 61 and ENDF File 6, LAW=1, LANG/=2. +!=============================================================================== + type AngleEnergyTable integer :: interpolation integer :: n_discrete @@ -16,11 +21,11 @@ module secondary_correlated end type AngleEnergyTable type, extends(AngleEnergy) :: CorrelatedAngleEnergy - integer :: n_region - integer, allocatable :: breakpoints(:) - integer, allocatable :: interpolation(:) - real(8), allocatable :: energy_in(:) - type(AngleEnergyTable), allocatable :: table(:) + integer :: n_region ! number of interpolation regions + integer, allocatable :: breakpoints(:) ! breakpoints of interpolation regions + integer, allocatable :: interpolation(:) ! interpolation region codes + real(8), allocatable :: energy_in(:) ! incoming energies + type(AngleEnergyTable), allocatable :: table(:) ! outgoing E/mu distributions contains procedure :: sample => correlated_sample end type CorrelatedAngleEnergy @@ -29,13 +34,13 @@ contains subroutine correlated_sample(this, E_in, E_out, mu) class(CorrelatedAngleEnergy), intent(in) :: this - real(8), intent(in) :: E_in - real(8), intent(out) :: E_out - real(8), intent(out) :: mu + real(8), intent(in) :: E_in ! incoming energy + real(8), intent(out) :: E_out ! sampled outgoing energy + real(8), intent(out) :: mu ! sapmled scattering cosine - integer :: i, k, l - integer :: n_energy_in - integer :: n_energy_out + integer :: i, k, l ! indices + integer :: n_energy_in ! number of incoming energies + integer :: n_energy_out ! number of outgoing energies real(8) :: r ! interpolation factor on incoming energy real(8) :: r1 ! random number on [0,1) real(8) :: frac ! interpolation factor on outgoing energy diff --git a/src/secondary_header.F90 b/src/secondary_header.F90 index e0c414826b..d858074129 100644 --- a/src/secondary_header.F90 +++ b/src/secondary_header.F90 @@ -4,15 +4,18 @@ module secondary_header use interpolation, only: interpolate_tab1 use random_lcg, only: prn +!=============================================================================== +! ANGLEENERGY (abstract) defines a correlated or uncorrelated angle-energy +! distribution that is a function of incoming energy. Each derived type must +! implement a sample() subroutine that returns an outgoing energy and scattering +! cosine given an incoming energy. +!=============================================================================== + type, abstract :: AngleEnergy contains procedure(iSampleAngleEnergy), deferred :: sample end type AngleEnergy - type :: AngleEnergyContainer - class(AngleEnergy), allocatable :: obj - end type AngleEnergyContainer - abstract interface subroutine iSampleAngleEnergy(this, E_in, E_out, mu) import AngleEnergy @@ -23,9 +26,16 @@ module secondary_header end subroutine iSampleAngleEnergy end interface + type :: AngleEnergyContainer + class(AngleEnergy), allocatable :: obj + end type AngleEnergyContainer + !=============================================================================== -! SECONDARYDISTRIBUTION stores a secondary distribution for angle and energy, -! whether correlated or uncorrelated. +! SECONDARYDISTRIBUTION stores multiple angle-energy distributions, each of +! which has a given probability of occurring for a given incoming energy. In +! general, most secondary distributions only have one angle-energy distribution, +! but for some cases (e.g., (n,2n) in certain nuclides) multiple distinct +! distributions exist. !=============================================================================== type :: SecondaryDistribution @@ -39,12 +49,12 @@ contains subroutine secondary_sample(this, E_in, E_out, mu) class(SecondaryDistribution), intent(in) :: this - real(8), intent(in) :: E_in - real(8), intent(out) :: E_out - real(8), intent(out) :: mu + real(8), intent(in) :: E_in ! incoming energy + real(8), intent(out) :: E_out ! sampled outgoing energy + real(8), intent(out) :: mu ! sampled scattering cosine - integer :: n - real(8) :: p_valid + integer :: n ! number of angle-energy distributions + real(8) :: p_valid ! probability that given distribution is valid n = size(this%applicability) if (n > 1) then diff --git a/src/secondary_kalbach.F90 b/src/secondary_kalbach.F90 index d38541edd5..5e6949206e 100644 --- a/src/secondary_kalbach.F90 +++ b/src/secondary_kalbach.F90 @@ -5,6 +5,12 @@ module secondary_kalbach use random_lcg, only: prn use search, only: binary_search +!=============================================================================== +! KalbachMann represents a correlated angle-energy distribution with the angular +! distribution represented using Kalbach-Mann systematics. This corresponds to +! ACE law 44 and ENDF File 6, LAW=1, LANG=2. +!=============================================================================== + type KalbachMannTable integer :: n_discrete integer :: interpolation @@ -16,11 +22,11 @@ module secondary_kalbach end type KalbachMannTable type, extends(AngleEnergy) :: KalbachMann - integer :: n_region - integer, allocatable :: breakpoints(:) - integer, allocatable :: interpolation(:) - real(8), allocatable :: energy_in(:) - type(KalbachMannTable), allocatable :: table(:) + integer :: n_region ! number of interpolation regions + integer, allocatable :: breakpoints(:) ! breakpoints of interpolation regions + integer, allocatable :: interpolation(:) ! interpolation region codes + real(8), allocatable :: energy_in(:) ! incoming energies + type(KalbachMannTable), allocatable :: table(:) ! outgoing E/mu parameters contains procedure :: sample => kalbachmann_sample end type KalbachMann @@ -29,13 +35,13 @@ contains subroutine kalbachmann_sample(this, E_in, E_out, mu) class(KalbachMann), intent(in) :: this - real(8), intent(in) :: E_in - real(8), intent(out) :: E_out - real(8), intent(out) :: mu + real(8), intent(in) :: E_in ! incoming energy + real(8), intent(out) :: E_out ! sampled outgoing energy + real(8), intent(out) :: mu ! sampled scattering cosine - integer :: i, k, l - integer :: n_energy_in - integer :: n_energy_out + integer :: i, k, l ! indices + integer :: n_energy_in ! number of incoming energies + integer :: n_energy_out ! number of outgoing energies real(8) :: r ! interpolation factor on incoming energy real(8) :: r1 ! random number on [0,1) real(8) :: frac ! interpolation factor on outgoing energy diff --git a/src/secondary_uncorrelated.F90 b/src/secondary_uncorrelated.F90 index 40fa53a382..22a56aa127 100644 --- a/src/secondary_uncorrelated.F90 +++ b/src/secondary_uncorrelated.F90 @@ -6,6 +6,12 @@ module secondary_uncorrelated use secondary_header, only: AngleEnergy use random_lcg, only: prn +!=============================================================================== +! UNCORRELATEDANGLEENERGY represents an uncorrelated angle-energy +! distribution. This corresponds to when an energy distribution is given in ENDF +! File 5/6 and an angular distribution is given in ENDF File 4. +!=============================================================================== + type, extends(AngleEnergy) :: UncorrelatedAngleEnergy logical :: fission = .false. type(AngleDistribution) :: angle @@ -18,9 +24,9 @@ contains subroutine uncorrelated_sample(this, E_in, E_out, mu) class(UncorrelatedAngleEnergy), intent(in) :: this - real(8), intent(in) :: E_in - real(8), intent(out) :: E_out - real(8), intent(out) :: mu + real(8), intent(in) :: E_in ! incoming energy + real(8), intent(out) :: E_out ! sampled outgoing energy + real(8), intent(out) :: mu ! sampled scattering cosine ! Sample cosine of scattering angle if (this%fission) then From 7bd0b2fe9749f3c6d6828588664b90a53ad965ff Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 22 Jan 2016 08:59:54 -0600 Subject: [PATCH 08/11] Cut down on print_nuclide routine --- src/output.F90 | 43 +++---------------------------------------- 1 file changed, 3 insertions(+), 40 deletions(-) diff --git a/src/output.F90 b/src/output.F90 index 3a48de02da..f133bad4bb 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -322,14 +322,8 @@ contains integer :: i ! loop index over nuclides integer :: unit_ ! unit to write to - integer :: size_total ! memory used by nuclide (bytes) - integer :: size_angle_total ! total memory used for angle dist. (bytes) - integer :: size_energy_total ! total memory used for energy dist. (bytes) integer :: size_xs ! memory used for cross-sections (bytes) - integer :: size_angle ! memory used for an angle distribution (bytes) - integer :: size_energy ! memory used for a energy distributions (bytes) integer :: size_urr ! memory used for probability tables (bytes) - character(11) :: law ! secondary energy distribution law type(UrrData), pointer :: urr ! set default unit for writing information @@ -340,8 +334,6 @@ contains end if ! Initialize totals - size_angle_total = 0 - size_energy_total = 0 size_urr = 0 size_xs = 0 @@ -356,36 +348,15 @@ contains write(unit_,*) ' # of reactions = ' // trim(to_str(nuc % n_reaction)) ! Information on each reaction - write(unit_,*) ' Reaction Q-value COM Law IE size(angle) size(energy)' + write(unit_,*) ' Reaction Q-value COM IE' do i = 1, nuc % n_reaction associate (rxn => nuc % reactions(i)) -!!$ ! Determine size of angle distribution -!!$ if (rxn % has_angle_dist) then -!!$ size_angle = rxn % adist % n_energy * 16 + size(rxn % adist % data) * 8 -!!$ else -!!$ size_angle = 0 -!!$ end if - size_angle = 0 - -!!$ ! Determine size of energy distribution and law -!!$ if (rxn % has_energy_dist) then -!!$ size_energy = size(rxn % edist % data) * 8 -!!$ law = to_str(rxn % edist % law) -!!$ else -!!$ size_energy = 0 -!!$ law = 'None' -!!$ end if - size_energy = 0 - law = 'None' - - write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,A4,1X,I6,1X,I11,1X,I11)') & + write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,I6)') & reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, & - law(1:4), rxn % threshold, size_angle, size_energy + rxn % threshold ! Accumulate data size size_xs = size_xs + (nuc % n_grid - rxn%threshold + 1) * 8 - size_angle_total = size_angle_total + size_angle - size_energy_total = size_energy_total + size_energy end associate end do @@ -411,19 +382,11 @@ contains size_urr = urr % n_energy * (urr % n_prob * 6 + 1) * 8 end if - ! Calculate total memory - size_total = size_xs + size_angle_total + size_energy_total + size_urr - ! Write memory used write(unit_,*) ' Memory Requirements' write(unit_,*) ' Cross sections = ' // trim(to_str(size_xs)) // ' bytes' - write(unit_,*) ' Secondary angle distributions = ' // & - trim(to_str(size_angle_total)) // ' bytes' - write(unit_,*) ' Secondary energy distributions = ' // & - trim(to_str(size_energy_total)) // ' bytes' write(unit_,*) ' Probability Tables = ' // & trim(to_str(size_urr)) // ' bytes' - write(unit_,*) ' Total = ' // trim(to_str(size_total)) // ' bytes' ! Blank line at end of nuclide write(unit_,*) From 869456ddcffbdd2a3546aeb756fab9ddb2563164 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 22 Jan 2016 09:45:43 -0600 Subject: [PATCH 09/11] Fix comment Evaporation -> WattEnergy --- src/energy_distribution.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/energy_distribution.F90 b/src/energy_distribution.F90 index 1dea1bfece..db3dcc4411 100644 --- a/src/energy_distribution.F90 +++ b/src/energy_distribution.F90 @@ -114,7 +114,7 @@ module energy_distribution end type Evaporation !=============================================================================== -! EVAPORATION gives the energy distribution of neutrons emitted from a Watt +! WATTENERGY gives the energy distribution of neutrons emitted from a Watt ! fission spectrum. This corresponds to ACE law 11 and ENDF File 5, LF=11. !=============================================================================== From c635084221860ee0347ae7699e955e644953aee9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 22 Jan 2016 11:38:47 -0600 Subject: [PATCH 10/11] Add test to cover Watt, N-body phase space, and evaporation spectra. Fix bug with N-body phase space. --- src/ace.F90 | 2 +- tests/test_energy_laws/geometry.xml | 7 +++++++ tests/test_energy_laws/materials.xml | 11 +++++++++++ tests/test_energy_laws/results_true.dat | 2 ++ tests/test_energy_laws/settings.xml | 14 ++++++++++++++ tests/test_energy_laws/test_energy_laws.py | 21 +++++++++++++++++++++ 6 files changed, 56 insertions(+), 1 deletion(-) create mode 100644 tests/test_energy_laws/geometry.xml create mode 100644 tests/test_energy_laws/materials.xml create mode 100644 tests/test_energy_laws/results_true.dat create mode 100644 tests/test_energy_laws/settings.xml create mode 100644 tests/test_energy_laws/test_energy_laws.py diff --git a/src/ace.F90 b/src/ace.F90 index 5b003d4067..1c49262985 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -1033,7 +1033,7 @@ contains ! Read energy law data call get_energy_dist(secondary%distribution(n)%obj, LAW, & - JXS(11), IDAT, nuc%awr, nuc%reactions(i)%Q_value) + JXS(11), IDAT, nuc%awr, nuc%reactions(i + 1)%Q_value) ! <<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<< ! Before the secondary distribution refactor, when the angle/energy diff --git a/tests/test_energy_laws/geometry.xml b/tests/test_energy_laws/geometry.xml new file mode 100644 index 0000000000..bb0e15d47f --- /dev/null +++ b/tests/test_energy_laws/geometry.xml @@ -0,0 +1,7 @@ + + + + + + + diff --git a/tests/test_energy_laws/materials.xml b/tests/test_energy_laws/materials.xml new file mode 100644 index 0000000000..1a40aeb1db --- /dev/null +++ b/tests/test_energy_laws/materials.xml @@ -0,0 +1,11 @@ + + + + + + + + + + + diff --git a/tests/test_energy_laws/results_true.dat b/tests/test_energy_laws/results_true.dat new file mode 100644 index 0000000000..777d0790d4 --- /dev/null +++ b/tests/test_energy_laws/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +2.421637E+00 5.324089E-03 diff --git a/tests/test_energy_laws/settings.xml b/tests/test_energy_laws/settings.xml new file mode 100644 index 0000000000..9f0e8ed05f --- /dev/null +++ b/tests/test_energy_laws/settings.xml @@ -0,0 +1,14 @@ + + + + + 10 + 5 + 1000 + + + + + + + diff --git a/tests/test_energy_laws/test_energy_laws.py b/tests/test_energy_laws/test_energy_laws.py new file mode 100644 index 0000000000..5f436e4694 --- /dev/null +++ b/tests/test_energy_laws/test_energy_laws.py @@ -0,0 +1,21 @@ +#!/usr/bin/env python + +"""The purpose of this test is to provide coverage of Watt, N-body phase space, +and evaporation energy distributions. The only nuclide that uses a Watt fission +spectrum in ENDF/B-VII.1 is U-233. The only nuclide that has a reaction using +the N-body phase space distribution is H-2(n.2n). Several nuclides have +reactions with evaporation spectra. In this test, the material is composed of +U-233, H-2, and Na-23 (which has several reactions with evaporation spectra. + +""" + +import glob +import os +import sys +sys.path.insert(0, os.pardir) +from testing_harness import TestHarness + + +if __name__ == '__main__': + harness = TestHarness('statepoint.10.*') + harness.main() From 2ff33e55111a6ad7b9720c558f5918b05d52afd1 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 22 Jan 2016 13:29:59 -0600 Subject: [PATCH 11/11] Modify test_energy_laws to include Ta-181 Ta-181 is among a handful of nuclides that uses Kalbach-Mann systematics with linear-linear interpolation. --- tests/test_energy_laws/geometry.xml | 2 -- tests/test_energy_laws/materials.xml | 10 +++++----- tests/test_energy_laws/results_true.dat | 2 +- tests/test_energy_laws/settings.xml | 3 --- tests/test_energy_laws/test_energy_laws.py | 21 +++++++++++++++------ 5 files changed, 21 insertions(+), 17 deletions(-) diff --git a/tests/test_energy_laws/geometry.xml b/tests/test_energy_laws/geometry.xml index bb0e15d47f..c42f45597d 100644 --- a/tests/test_energy_laws/geometry.xml +++ b/tests/test_energy_laws/geometry.xml @@ -1,7 +1,5 @@ - - diff --git a/tests/test_energy_laws/materials.xml b/tests/test_energy_laws/materials.xml index 1a40aeb1db..97f9b84be1 100644 --- a/tests/test_energy_laws/materials.xml +++ b/tests/test_energy_laws/materials.xml @@ -1,11 +1,11 @@ - + 71c - - - + + + + - diff --git a/tests/test_energy_laws/results_true.dat b/tests/test_energy_laws/results_true.dat index 777d0790d4..48eb6bc81b 100644 --- a/tests/test_energy_laws/results_true.dat +++ b/tests/test_energy_laws/results_true.dat @@ -1,2 +1,2 @@ k-combined: -2.421637E+00 5.324089E-03 +2.130076E+00 1.938907E-03 diff --git a/tests/test_energy_laws/settings.xml b/tests/test_energy_laws/settings.xml index 9f0e8ed05f..1c3f444f0f 100644 --- a/tests/test_energy_laws/settings.xml +++ b/tests/test_energy_laws/settings.xml @@ -1,14 +1,11 @@ - 10 5 1000 - - diff --git a/tests/test_energy_laws/test_energy_laws.py b/tests/test_energy_laws/test_energy_laws.py index 5f436e4694..0e593876fc 100644 --- a/tests/test_energy_laws/test_energy_laws.py +++ b/tests/test_energy_laws/test_energy_laws.py @@ -1,11 +1,20 @@ #!/usr/bin/env python -"""The purpose of this test is to provide coverage of Watt, N-body phase space, -and evaporation energy distributions. The only nuclide that uses a Watt fission -spectrum in ENDF/B-VII.1 is U-233. The only nuclide that has a reaction using -the N-body phase space distribution is H-2(n.2n). Several nuclides have -reactions with evaporation spectra. In this test, the material is composed of -U-233, H-2, and Na-23 (which has several reactions with evaporation spectra. +"""The purpose of this test is to provide coverage of energy distributions that +are not covered in other tests. It has a single material with the following +nuclides: + +U-233: Only nuclide that has a Watt fission spectrum + +H-2: Only nuclide that has an N-body phase space distribution, in this case for +(n,2n) + +Na-23: Has an evaporation spectrum and also has reactions that have multiple +angle-energy distributions, so it provides coverage for both of those +situations. + +Ta-181: One of a few nuclides that has reactions with Kalbach-Mann distributions +that use linear-linear interpolation. """