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

Merge pull request #557 from paulromano/secondary-refactor

Browse source 
Refactor of secondary angle/energy distributions
This commit is contained in:
Sterling Harper 2016-01-23 19:41:45 -05:00
commit c4e5f4380f
20 changed files with 1508 additions and 1396 deletions
Download patch file Download diff file Expand all files Collapse all files

770
src/ace.F90
View file

@ -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,96 @@ 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
! 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(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)
end select
end do
end associate
end do
@ -981,25 +998,62 @@ 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 + 1)%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
end associate
end do
@ -1011,281 +1065,319 @@ 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
! 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(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)
! 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.

77
src/ace_header.F90
View file

@ -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)

63
src/angle_distribution.F90 Normal file
View file

@ -0,0 +1,63 @@
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
!===============================================================================
! 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(:)
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 ! incoming energy
real(8) :: mu ! sampled cosine of scattering angle
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
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

29
src/distribution_univariate.F90
View file

@ -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

34
src/endf_header.F90
View file

@ -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

429
src/energy_distribution.F90 Normal file
View file

@ -0,0 +1,429 @@
module energy_distribution
use constants, only: ZERO, ONE, TWO, PI, HISTOGRAM, LINEAR_LINEAR
use endf_header, only: Tab1
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. Each derived type must
! implement a sample() function that returns a sampled outgoing energy given an
! incoming energy
!===============================================================================
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
type :: EnergyDistributionContainer
class(EnergyDistribution), allocatable :: obj
end type EnergyDistributionContainer
!===============================================================================
! 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 ! 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
contains
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
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
!===============================================================================
! 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 ! 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
contains
procedure :: sample => evaporation_sample
end type Evaporation
!===============================================================================
! 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.
!===============================================================================
type, extends(EnergyDistribution) :: WattEnergy
type(Tab1) :: a
type(Tab1) :: b
real(8) :: u
contains
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
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 ! incoming energy
real(8) :: E_out ! sampled outgoing energy
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
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)
! 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 ! incoming energy
real(8) :: E_out ! sampled outgoing energy
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
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
! 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)
! 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 ! incoming energy
real(8) :: E_out ! sampled outgoing energy
real(8) :: theta ! Maxwell distribution parameter
! 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 ! incoming energy
real(8) :: E_out ! sampled outgoing energy
real(8) :: theta ! evaporation spectrum parameter
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 ! incoming energy
real(8) :: E_out ! sampled outgoing energy
real(8) :: a, b ! Watt spectrum parameters
! 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 ! incoming energy
real(8) :: E_out ! sampled outgoing energy
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
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

1
src/global.F90
View file

@ -22,7 +22,6 @@ module global
#endif
implicit none
save
! ============================================================================
! GEOMETRY-RELATED VARIABLES

1
src/interpolation.F90
View file

@ -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

40
src/output.F90
View file

@ -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,33 +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
! 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
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
@ -408,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_,*)

962
src/physics.F90
View file

File diff suppressed because it is too large Load diff

1
src/search.F90
View file

@ -1,7 +1,6 @@
module search
use constants
use error, only: fatal_error
implicit none

148
src/secondary_correlated.F90 Normal file
View file

@ -0,0 +1,148 @@
module secondary_correlated
use constants, only: ZERO, ONE, TWO, HISTOGRAM, LINEAR_LINEAR
use distribution_univariate, only: DistributionContainer
use secondary_header, only: AngleEnergy
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
real(8), allocatable :: e_out(:)
real(8), allocatable :: p(:)
real(8), allocatable :: c(:)
type(DistributionContainer), allocatable :: angle(:)
end type AngleEnergyTable
type, extends(AngleEnergy) :: CorrelatedAngleEnergy
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
contains
subroutine correlated_sample(this, E_in, E_out, mu)
class(CorrelatedAngleEnergy), intent(in) :: this
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 ! 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
! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 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)
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)
! 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
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
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)%obj%sample()
else
mu = this%table(l)%angle(k + 1)%obj%sample()
end if
end subroutine correlated_sample
end module secondary_correlated

78
src/secondary_header.F90 Normal file
View file

@ -0,0 +1,78 @@
module secondary_header
use endf_header, only: Tab1
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
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
type :: AngleEnergyContainer
class(AngleEnergy), allocatable :: obj
end type AngleEnergyContainer
!===============================================================================
! 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
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)
class(SecondaryDistribution), intent(in) :: this
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 ! number of angle-energy distributions
real(8) :: p_valid ! probability that given distribution is valid
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
end subroutine secondary_sample
end module secondary_header

164
src/secondary_kalbach.F90 Normal file
View file

@ -0,0 +1,164 @@
module secondary_kalbach
use constants, only: ZERO, ONE, TWO, HISTOGRAM, LINEAR_LINEAR
use secondary_header, only: AngleEnergy
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
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(AngleEnergy) :: KalbachMann
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
contains
subroutine kalbachmann_sample(this, E_in, E_out, mu)
class(KalbachMann), intent(in) :: this
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 ! 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
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)
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)
! 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

48
src/secondary_uncorrelated.F90 Normal file
View file

@ -0,0 +1,48 @@
module secondary_uncorrelated
use angle_distribution, only: AngleDistribution
use constants, only: ONE, TWO
use energy_distribution, only: EnergyDistribution
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
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 ! 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
! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 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
mu = TWO*prn() - ONE
end if
! Sample outgoing energy
E_out = this%energy%sample(E_in)
end subroutine uncorrelated_sample
end module secondary_uncorrelated

5
tests/test_energy_laws/geometry.xml Normal file
View file

@ -0,0 +1,5 @@
<?xml version="1.0"?>
<geometry>
<surface id="1" type="sphere" coeffs="0 0 0 100" boundary="vacuum"/>
<cell id="1" material="1" region="-1" />
</geometry>

11
tests/test_energy_laws/materials.xml Normal file
View file

@ -0,0 +1,11 @@
<?xml version="1.0"?>
<materials>
<default_xs>71c</default_xs>
<material id="1">
<density value="20" units="g/cc" />
<nuclide name="U-233" ao="1.0" />
<nuclide name="H-2" ao="1.0" />
<nuclide name="Na-23" ao="1.0" />
<nuclide name="Ta-181" ao="1.0" />
</material>
</materials>

2
tests/test_energy_laws/results_true.dat Normal file
View file

@ -0,0 +1,2 @@
k-combined:
2.130076E+00 1.938907E-03

11
tests/test_energy_laws/settings.xml Normal file
View file

@ -0,0 +1,11 @@
<?xml version="1.0"?>
<settings>
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
</eigenvalue>
<source>
<space type="point" parameters="0. 0. 0." />
</source>
</settings>

30
tests/test_energy_laws/test_energy_laws.py Normal file
View file

@ -0,0 +1,30 @@
#!/usr/bin/env python
"""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.
"""
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()