mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-29 06:35:48 -04:00
Remove unused angle/energy distribution Fortran files
This commit is contained in:
Paul Romano
parent
4c30c61dc0
commit
56b5e733ac
9 changed files with 0 additions and 1333 deletions
|
|
@ -288,7 +288,6 @@ set_target_properties(faddeeva PROPERTIES
|
|||
|
||||
add_library(libopenmc SHARED
|
||||
src/algorithm.F90
|
||||
src/angle_distribution.F90
|
||||
src/angleenergy_header.F90
|
||||
src/bank_header.F90
|
||||
src/api.F90
|
||||
|
|
@ -306,7 +305,6 @@ add_library(libopenmc SHARED
|
|||
src/eigenvalue.F90
|
||||
src/endf.F90
|
||||
src/endf_header.F90
|
||||
src/energy_distribution.F90
|
||||
src/error.F90
|
||||
src/geometry.F90
|
||||
src/geometry_header.F90
|
||||
|
|
@ -334,7 +332,6 @@ add_library(libopenmc SHARED
|
|||
src/physics_mg.F90
|
||||
src/plot.F90
|
||||
src/plot_header.F90
|
||||
src/product_header.F90
|
||||
src/progress_header.F90
|
||||
src/pugixml/pugixml_f.F90
|
||||
src/random_lcg.F90
|
||||
|
|
@ -342,9 +339,6 @@ add_library(libopenmc SHARED
|
|||
src/relaxng
|
||||
src/sab_header.F90
|
||||
src/secondary_correlated.F90
|
||||
src/secondary_kalbach.F90
|
||||
src/secondary_nbody.F90
|
||||
src/secondary_uncorrelated.F90
|
||||
src/set_header.F90
|
||||
src/settings.F90
|
||||
src/simulation_header.F90
|
||||
|
|
|
|||
|
|
@ -1,130 +0,0 @@
|
|||
module angle_distribution
|
||||
|
||||
use algorithm, only: binary_search
|
||||
use constants, only: ZERO, ONE, HISTOGRAM, LINEAR_LINEAR
|
||||
use distribution_univariate, only: DistributionContainer, Tabular
|
||||
use hdf5_interface, only: read_attribute, get_shape, read_dataset, &
|
||||
open_dataset, close_dataset, HID_T, HSIZE_T
|
||||
use random_lcg, only: prn
|
||||
|
||||
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
|
||||
procedure :: from_hdf5 => angle_from_hdf5
|
||||
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
|
||||
|
||||
subroutine angle_from_hdf5(this, group_id)
|
||||
class(AngleDistribution), intent(inout) :: this
|
||||
integer(HID_T), intent(in) :: group_id
|
||||
|
||||
integer :: i, j
|
||||
integer :: n
|
||||
integer :: n_energy
|
||||
integer(HID_T) :: dset_id
|
||||
integer(HSIZE_T) :: dims(1), dims2(2)
|
||||
integer, allocatable :: offsets(:)
|
||||
integer, allocatable :: interp(:)
|
||||
real(8), allocatable :: temp(:,:)
|
||||
|
||||
! Get incoming energies
|
||||
dset_id = open_dataset(group_id, 'energy')
|
||||
call get_shape(dset_id, dims)
|
||||
n_energy = int(dims(1), 4)
|
||||
allocate(this % energy(n_energy))
|
||||
allocate(this % distribution(n_energy))
|
||||
call read_dataset(this % energy, dset_id)
|
||||
call close_dataset(dset_id)
|
||||
|
||||
! Get outgoing energy distribution data
|
||||
dset_id = open_dataset(group_id, 'mu')
|
||||
call read_attribute(offsets, dset_id, 'offsets')
|
||||
call read_attribute(interp, dset_id, 'interpolation')
|
||||
call get_shape(dset_id, dims2)
|
||||
allocate(temp(dims2(1), dims2(2)))
|
||||
call read_dataset(temp, dset_id)
|
||||
call close_dataset(dset_id)
|
||||
|
||||
do i = 1, n_energy
|
||||
! Determine number of outgoing energies
|
||||
j = offsets(i)
|
||||
if (i < n_energy) then
|
||||
n = offsets(i+1) - j
|
||||
else
|
||||
n = size(temp, 1) - j
|
||||
end if
|
||||
|
||||
! Create and initialize tabular distribution
|
||||
allocate(Tabular :: this % distribution(i) % obj)
|
||||
select type (mudist => this % distribution(i) % obj)
|
||||
type is (Tabular)
|
||||
mudist % interpolation = interp(i)
|
||||
allocate(mudist % x(n), mudist % p(n), mudist % c(n))
|
||||
mudist % x(:) = temp(j+1:j+n, 1)
|
||||
mudist % p(:) = temp(j+1:j+n, 2)
|
||||
|
||||
! To get answers that match ACE data, for now we still use the tabulated
|
||||
! CDF values that were passed through to the HDF5 library. At a later
|
||||
! time, we can remove the CDF values from the HDF5 library and
|
||||
! reconstruct them using the PDF
|
||||
if (.true.) then
|
||||
mudist % c(:) = temp(j+1:j+n, 3)
|
||||
else
|
||||
call mudist % initialize(temp(j+1:j+n, 1), temp(j+1:j+n, 2), interp(i))
|
||||
end if
|
||||
end select
|
||||
|
||||
j = j + n
|
||||
end do
|
||||
end subroutine angle_from_hdf5
|
||||
|
||||
end module angle_distribution
|
||||
|
|
@ -1,583 +0,0 @@
|
|||
module energy_distribution
|
||||
|
||||
use algorithm, only: binary_search
|
||||
use constants, only: ZERO, ONE, HALF, TWO, PI, HISTOGRAM, LINEAR_LINEAR
|
||||
use endf_header, only: Tabulated1D
|
||||
use hdf5_interface
|
||||
use math, only: maxwell_spectrum, watt_spectrum
|
||||
use random_lcg, only: prn
|
||||
|
||||
!===============================================================================
|
||||
! 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(energy_distribution_sample_), deferred :: sample
|
||||
procedure(energy_distribution_from_hdf5_), deferred :: from_hdf5
|
||||
end type EnergyDistribution
|
||||
|
||||
abstract interface
|
||||
function energy_distribution_sample_(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 energy_distribution_sample_
|
||||
|
||||
subroutine energy_distribution_from_hdf5_(this, group_id)
|
||||
import EnergyDistribution
|
||||
import HID_T
|
||||
class(EnergyDistribution), intent(inout) :: this
|
||||
integer(HID_T), intent(in) :: group_id
|
||||
end subroutine energy_distribution_from_hdf5_
|
||||
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
|
||||
procedure :: from_hdf5 => equiprobable_from_hdf5
|
||||
end type TabularEquiprobable
|
||||
|
||||
!===============================================================================
|
||||
! DISCRETEPHOTON gives the energy distribution for a discrete photon (usually
|
||||
! used for photon production from an incident-neutron reaction)
|
||||
!===============================================================================
|
||||
|
||||
type, extends(EnergyDistribution) :: DiscretePhoton
|
||||
integer :: primary_flag
|
||||
real(8) :: energy
|
||||
real(8) :: A
|
||||
contains
|
||||
procedure :: sample => discrete_photon_sample
|
||||
procedure :: from_hdf5 => discrete_photon_from_hdf5
|
||||
end type DiscretePhoton
|
||||
|
||||
!===============================================================================
|
||||
! 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
|
||||
procedure :: from_hdf5 => level_inelastic_from_hdf5
|
||||
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(:)
|
||||
type(CTTable), allocatable :: distribution(:)
|
||||
contains
|
||||
procedure :: sample => continuous_sample
|
||||
procedure :: from_hdf5 => continuous_from_hdf5
|
||||
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(Tabulated1D) :: theta ! incoming-energy-dependent parameter
|
||||
real(8) :: u ! restriction energy
|
||||
contains
|
||||
procedure :: sample => maxwellenergy_sample
|
||||
procedure :: from_hdf5 => maxwellenergy_from_hdf5
|
||||
end type MaxwellEnergy
|
||||
|
||||
!===============================================================================
|
||||
! EVAPORATION represents an evaporation spectrum corresponding to ACE law 9 and
|
||||
! ENDF File 5, LF=9.
|
||||
!===============================================================================
|
||||
|
||||
type, extends(EnergyDistribution) :: Evaporation
|
||||
type(Tabulated1D) :: theta
|
||||
real(8) :: u
|
||||
contains
|
||||
procedure :: sample => evaporation_sample
|
||||
procedure :: from_hdf5 => evaporation_from_hdf5
|
||||
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(Tabulated1D) :: a
|
||||
type(Tabulated1D) :: b
|
||||
real(8) :: u
|
||||
contains
|
||||
procedure :: sample => watt_sample
|
||||
procedure :: from_hdf5 => watt_from_hdf5
|
||||
end type WattEnergy
|
||||
|
||||
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
|
||||
|
||||
subroutine equiprobable_from_hdf5(this, group_id)
|
||||
class(TabularEquiprobable), intent(inout) :: this
|
||||
integer(HID_T), intent(in) :: group_id
|
||||
end subroutine equiprobable_from_hdf5
|
||||
|
||||
function discrete_photon_sample(this, E_in) result(E_out)
|
||||
class(DiscretePhoton), intent(in) :: this
|
||||
real(8), intent(in) :: E_in
|
||||
real(8) :: E_out
|
||||
|
||||
if (this % primary_flag == 2) then
|
||||
E_out = this % energy + this % A/(this % A + 1)*E_in
|
||||
else
|
||||
E_out = this % energy
|
||||
end if
|
||||
end function discrete_photon_sample
|
||||
|
||||
subroutine discrete_photon_from_hdf5(this, group_id)
|
||||
class(DiscretePhoton), intent(inout) :: this
|
||||
integer(HID_T), intent(in) :: group_id
|
||||
|
||||
call read_attribute(this % primary_flag, group_id, 'primary_flag')
|
||||
call read_attribute(this % energy, group_id, 'energy')
|
||||
call read_attribute(this % A, group_id, 'atomic_weight_ratio')
|
||||
end subroutine discrete_photon_from_hdf5
|
||||
|
||||
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
|
||||
|
||||
subroutine level_inelastic_from_hdf5(this, group_id)
|
||||
class(LevelInelastic), intent(inout) :: this
|
||||
integer(HID_T), intent(in) :: group_id
|
||||
|
||||
call read_attribute(this%threshold, group_id, 'threshold')
|
||||
call read_attribute(this%mass_ratio, group_id, 'mass_ratio')
|
||||
end subroutine level_inelastic_from_hdf5
|
||||
|
||||
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)
|
||||
if (E_in < this%energy(1)) then
|
||||
i = 1
|
||||
r = ZERO
|
||||
elseif (E_in > this%energy(n_energy_in)) then
|
||||
i = n_energy_in - 1
|
||||
r = ONE
|
||||
else
|
||||
i = binary_search(this%energy, n_energy_in, E_in)
|
||||
r = (E_in - this%energy(i)) / &
|
||||
(this%energy(i+1) - this%energy(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%distribution(i)%e_out)
|
||||
E_i_1 = this%distribution(i)%e_out(1)
|
||||
E_i_K = this%distribution(i)%e_out(n_energy_out)
|
||||
|
||||
n_energy_out = size(this%distribution(i+1)%e_out)
|
||||
E_i1_1 = this%distribution(i+1)%e_out(1)
|
||||
E_i1_K = this%distribution(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%distribution(l)%e_out)
|
||||
r1 = prn()
|
||||
c_k = this%distribution(l)%c(1)
|
||||
do k = 1, n_energy_out - 1
|
||||
c_k1 = this%distribution(l)%c(k+1)
|
||||
if (r1 < c_k1) exit
|
||||
c_k = c_k1
|
||||
end do
|
||||
|
||||
! Check to make sure 1 <= k <= NP - 1
|
||||
k = max(1, min(k, n_energy_out - 1))
|
||||
|
||||
E_l_k = this%distribution(l)%e_out(k)
|
||||
p_l_k = this%distribution(l)%p(k)
|
||||
if (this%distribution(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%distribution(l)%interpolation == LINEAR_LINEAR) then
|
||||
! Linear-linear interpolation
|
||||
E_l_k1 = this%distribution(l)%e_out(k+1)
|
||||
p_l_k1 = this%distribution(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 .and. n_energy_out > 1) 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
|
||||
|
||||
subroutine continuous_from_hdf5(this, group_id)
|
||||
class(ContinuousTabular), intent(inout) :: this
|
||||
integer(HID_T), intent(in) :: group_id
|
||||
|
||||
integer :: i, j, k
|
||||
integer :: n
|
||||
integer :: n_energy
|
||||
integer(HID_T) :: dset_id
|
||||
integer(HSIZE_T) :: dims(1), dims2(2)
|
||||
integer, allocatable :: temp(:,:)
|
||||
integer, allocatable :: offsets(:)
|
||||
integer, allocatable :: interp(:)
|
||||
integer, allocatable :: n_discrete(:)
|
||||
real(8), allocatable :: eout(:,:)
|
||||
|
||||
! Open incoming energy dataset
|
||||
dset_id = open_dataset(group_id, 'energy')
|
||||
|
||||
! Get interpolation parameters
|
||||
call read_attribute(temp, dset_id, 'interpolation')
|
||||
allocate(this%breakpoints(size(temp, 1)))
|
||||
allocate(this%interpolation(size(temp, 1)))
|
||||
this%breakpoints(:) = temp(:, 1)
|
||||
this%interpolation(:) = temp(:, 2)
|
||||
this%n_region = size(this%breakpoints)
|
||||
|
||||
! Get incoming energies
|
||||
call get_shape(dset_id, dims)
|
||||
n_energy = int(dims(1), 4)
|
||||
allocate(this%energy(n_energy))
|
||||
allocate(this%distribution(n_energy))
|
||||
call read_dataset(this%energy, dset_id)
|
||||
call close_dataset(dset_id)
|
||||
|
||||
! Get outgoing energy distribution data
|
||||
dset_id = open_dataset(group_id, 'distribution')
|
||||
call read_attribute(offsets, dset_id, 'offsets')
|
||||
call read_attribute(interp, dset_id, 'interpolation')
|
||||
call read_attribute(n_discrete, dset_id, 'n_discrete_lines')
|
||||
call get_shape(dset_id, dims2)
|
||||
allocate(eout(dims2(1), dims2(2)))
|
||||
call read_dataset(eout, dset_id)
|
||||
call close_dataset(dset_id)
|
||||
|
||||
do i = 1, n_energy
|
||||
! Determine number of outgoing energies
|
||||
j = offsets(i)
|
||||
if (i < n_energy) then
|
||||
n = offsets(i+1) - j
|
||||
else
|
||||
n = size(eout, 1) - j
|
||||
end if
|
||||
|
||||
associate (d => this % distribution(i))
|
||||
! Assign interpolation scheme and number of discrete lines
|
||||
d % interpolation = interp(i)
|
||||
d % n_discrete = n_discrete(i)
|
||||
|
||||
! Allocate arrays for energies and PDF/CDF
|
||||
allocate(d % e_out(n))
|
||||
allocate(d % p(n))
|
||||
allocate(d % c(n))
|
||||
|
||||
! Copy data
|
||||
d % e_out(:) = eout(j+1:j+n, 1)
|
||||
d % p(:) = eout(j+1:j+n, 2)
|
||||
|
||||
! To get answers that match ACE data, for now we still use the tabulated
|
||||
! CDF values that were passed through to the HDF5 library. At a later
|
||||
! time, we can remove the CDF values from the HDF5 library and
|
||||
! reconstruct them using the PDF
|
||||
if (.true.) then
|
||||
d % c(:) = eout(j+1:j+n, 3)
|
||||
else
|
||||
! Calculate cumulative distribution function -- discrete portion
|
||||
do k = 1, n_discrete(i)
|
||||
if (k == 1) then
|
||||
d % c(k) = d % p(k)
|
||||
else
|
||||
d % c(k) = d % c(k-1) + d % p(k)
|
||||
end if
|
||||
end do
|
||||
|
||||
! Continuous portion
|
||||
do k = d % n_discrete + 1, n
|
||||
if (k == d % n_discrete + 1) then
|
||||
d % c(k) = sum(d % p(1:d % n_discrete))
|
||||
else
|
||||
if (d % interpolation == HISTOGRAM) then
|
||||
d % c(k) = d % c(k-1) + d % p(k-1) * &
|
||||
(d % e_out(k) - d % e_out(k-1))
|
||||
elseif (d % interpolation == LINEAR_LINEAR) then
|
||||
d % c(k) = d % c(k-1) + HALF*(d % p(k-1) + d % p(k)) * &
|
||||
(d % e_out(k) - d % e_out(k-1))
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
|
||||
! Normalize density and distribution functions
|
||||
d % p(:) = d % p(:)/d % c(n)
|
||||
d % c(:) = d % c(:)/d % c(n)
|
||||
end if
|
||||
end associate
|
||||
end do
|
||||
end subroutine continuous_from_hdf5
|
||||
|
||||
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 = this % theta % evaluate(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
|
||||
|
||||
subroutine maxwellenergy_from_hdf5(this, group_id)
|
||||
class(MaxwellEnergy), intent(inout) :: this
|
||||
integer(HID_T), intent(in) :: group_id
|
||||
|
||||
integer(HID_T) :: dset_id
|
||||
|
||||
call read_attribute(this%u, group_id, 'u')
|
||||
dset_id = open_dataset(group_id, 'theta')
|
||||
call this%theta%from_hdf5(dset_id)
|
||||
call close_dataset(dset_id)
|
||||
end subroutine maxwellenergy_from_hdf5
|
||||
|
||||
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 = this % theta % evaluate(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
|
||||
|
||||
subroutine evaporation_from_hdf5(this, group_id)
|
||||
class(Evaporation), intent(inout) :: this
|
||||
integer(HID_T), intent(in) :: group_id
|
||||
|
||||
integer(HID_T) :: dset_id
|
||||
|
||||
call read_attribute(this%u, group_id, 'u')
|
||||
dset_id = open_dataset(group_id, 'theta')
|
||||
call this%theta%from_hdf5(dset_id)
|
||||
call close_dataset(dset_id)
|
||||
end subroutine evaporation_from_hdf5
|
||||
|
||||
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 = this % a % evaluate(E_in)
|
||||
|
||||
! Determine Watt parameter 'b' from tabulated function
|
||||
b = this % b % evaluate(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
|
||||
|
||||
subroutine watt_from_hdf5(this, group_id)
|
||||
class(WattEnergy), intent(inout) :: this
|
||||
integer(HID_T), intent(in) :: group_id
|
||||
|
||||
integer(HID_T) :: dset_id
|
||||
|
||||
call read_attribute(this%u, group_id, 'u')
|
||||
|
||||
dset_id = open_dataset(group_id, 'a')
|
||||
call this%a%from_hdf5(dset_id)
|
||||
call close_dataset(dset_id)
|
||||
|
||||
dset_id = open_dataset(group_id, 'b')
|
||||
call this%b%from_hdf5(dset_id)
|
||||
call close_dataset(dset_id)
|
||||
end subroutine watt_from_hdf5
|
||||
|
||||
end module energy_distribution
|
||||
|
|
@ -17,11 +17,9 @@ module nuclide_header
|
|||
FIT_T, FIT_A, FIT_F, MultipoleArray
|
||||
use message_passing
|
||||
use multipole_header, only: MultipoleArray
|
||||
use product_header, only: AngleEnergyContainer
|
||||
use random_lcg, only: prn, future_prn, prn_set_stream
|
||||
use reaction_header, only: Reaction
|
||||
use sab_header, only: SAlphaBeta, sab_tables
|
||||
use secondary_uncorrelated, only: UncorrelatedAngleEnergy
|
||||
use settings
|
||||
use stl_vector, only: VectorInt, VectorReal
|
||||
use string
|
||||
|
|
|
|||
|
|
@ -18,7 +18,6 @@ module physics
|
|||
use random_lcg, only: prn, advance_prn_seed, prn_set_stream
|
||||
use reaction_header, only: Reaction
|
||||
use sab_header, only: sab_tables
|
||||
use secondary_uncorrelated, only: UncorrelatedAngleEnergy
|
||||
use settings
|
||||
use simulation_header
|
||||
use string, only: to_str
|
||||
|
|
|
|||
|
|
@ -1,153 +0,0 @@
|
|||
module product_header
|
||||
|
||||
use angleenergy_header, only: AngleEnergyContainer
|
||||
use constants, only: ZERO, MAX_WORD_LEN, EMISSION_PROMPT, EMISSION_DELAYED, &
|
||||
EMISSION_TOTAL, NEUTRON, PHOTON
|
||||
use endf_header, only: Tabulated1D, Function1D, Polynomial
|
||||
use hdf5_interface, only: read_attribute, open_group, close_group, &
|
||||
open_dataset, close_dataset, read_dataset, HID_T
|
||||
use random_lcg, only: prn
|
||||
use secondary_correlated, only: CorrelatedAngleEnergy
|
||||
use secondary_kalbach, only: KalbachMann
|
||||
use secondary_nbody, only: NBodyPhaseSpace
|
||||
use secondary_uncorrelated, only: UncorrelatedAngleEnergy
|
||||
use string, only: to_str
|
||||
|
||||
!===============================================================================
|
||||
! REACTIONPRODUCT stores a data for a reaction product including its yield and
|
||||
! angle-energy distributions, each of which has a given probability of occurring
|
||||
! for a given incoming energy. In general, most products only have one
|
||||
! angle-energy distribution, but for some cases (e.g., (n,2n) in certain
|
||||
! nuclides) multiple distinct distributions exist.
|
||||
!===============================================================================
|
||||
|
||||
type :: ReactionProduct
|
||||
integer :: particle
|
||||
integer :: emission_mode ! prompt, delayed, or total emission
|
||||
real(8) :: decay_rate ! Decay rate for delayed neutron precursors
|
||||
class(Function1D), pointer :: yield => null() ! Energy-dependent neutron yield
|
||||
type(Tabulated1D), allocatable :: applicability(:)
|
||||
type(AngleEnergyContainer), allocatable :: distribution(:)
|
||||
contains
|
||||
procedure :: sample => reactionproduct_sample
|
||||
procedure :: from_hdf5 => reactionproduct_from_hdf5
|
||||
end type ReactionProduct
|
||||
|
||||
contains
|
||||
|
||||
subroutine reactionproduct_sample(this, E_in, E_out, mu)
|
||||
class(ReactionProduct), 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 ! loop counter
|
||||
integer :: n ! number of angle-energy distributions
|
||||
real(8) :: prob ! cumulative probability
|
||||
real(8) :: c ! sampled cumulative probability
|
||||
|
||||
n = size(this%applicability)
|
||||
if (n > 1) then
|
||||
prob = ZERO
|
||||
c = prn()
|
||||
do i = 1, n
|
||||
! Determine probability that i-th energy distribution is sampled
|
||||
prob = prob + this % applicability(i) % evaluate(E_in)
|
||||
|
||||
! If i-th distribution is sampled, sample energy from the distribution
|
||||
if (c <= prob) 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 reactionproduct_sample
|
||||
|
||||
subroutine reactionproduct_from_hdf5(this, group_id)
|
||||
class(ReactionProduct), intent(inout) :: this
|
||||
integer(HID_T), intent(in) :: group_id
|
||||
|
||||
integer :: i
|
||||
integer :: n
|
||||
integer(HID_T) :: dgroup
|
||||
integer(HID_T) :: app
|
||||
integer(HID_T) :: yield
|
||||
character(MAX_WORD_LEN) :: temp
|
||||
|
||||
! Read particle type
|
||||
call read_attribute(temp, group_id, 'particle')
|
||||
select case (temp)
|
||||
case ('neutron')
|
||||
this % particle = NEUTRON
|
||||
case ('photon')
|
||||
this % particle = PHOTON
|
||||
end select
|
||||
|
||||
! Read emission mode and decay rate
|
||||
call read_attribute(temp, group_id, 'emission_mode')
|
||||
select case (temp)
|
||||
case ('prompt')
|
||||
this % emission_mode = EMISSION_PROMPT
|
||||
case ('delayed')
|
||||
this % emission_mode = EMISSION_DELAYED
|
||||
case ('total')
|
||||
this % emission_mode = EMISSION_TOTAL
|
||||
end select
|
||||
|
||||
! Read decay rate for delayed emission
|
||||
if (this % emission_mode == EMISSION_DELAYED) then
|
||||
call read_attribute(this % decay_rate, group_id, 'decay_rate')
|
||||
end if
|
||||
|
||||
! Read secondary particle yield
|
||||
yield = open_dataset(group_id, 'yield')
|
||||
call read_attribute(temp, yield, 'type')
|
||||
select case (temp)
|
||||
case ('Tabulated1D')
|
||||
allocate(Tabulated1D :: this % yield)
|
||||
case ('Polynomial')
|
||||
allocate(Polynomial :: this % yield)
|
||||
end select
|
||||
call this % yield % from_hdf5(yield)
|
||||
call close_dataset(yield)
|
||||
|
||||
call read_attribute(n, group_id, 'n_distribution')
|
||||
allocate(this%applicability(n))
|
||||
allocate(this%distribution(n))
|
||||
|
||||
do i = 1, n
|
||||
dgroup = open_group(group_id, trim('distribution_' // to_str(i - 1)))
|
||||
|
||||
! Read applicability
|
||||
if (n > 1) then
|
||||
app = open_dataset(dgroup, 'applicability')
|
||||
call this%applicability(i)%from_hdf5(app)
|
||||
call close_dataset(app)
|
||||
end if
|
||||
|
||||
! Read type of distribution and allocate accordingly
|
||||
call read_attribute(temp, dgroup, 'type')
|
||||
select case (temp)
|
||||
case ('uncorrelated')
|
||||
allocate(UncorrelatedAngleEnergy :: this%distribution(i)%obj)
|
||||
case ('correlated')
|
||||
allocate(CorrelatedAngleEnergy :: this%distribution(i)%obj)
|
||||
case ('nbody')
|
||||
allocate(NBodyPhaseSpace :: this%distribution(i)%obj)
|
||||
case ('kalbach-mann')
|
||||
allocate(KalbachMann :: this%distribution(i)%obj)
|
||||
end select
|
||||
|
||||
! Read distribution data
|
||||
call this%distribution(i)%obj%from_hdf5(dgroup)
|
||||
|
||||
call close_group(dgroup)
|
||||
end do
|
||||
|
||||
end subroutine reactionproduct_from_hdf5
|
||||
|
||||
end module product_header
|
||||
|
|
@ -1,278 +0,0 @@
|
|||
module secondary_kalbach
|
||||
|
||||
use algorithm, only: binary_search
|
||||
use angleenergy_header, only: AngleEnergy
|
||||
use constants, only: ZERO, HALF, ONE, TWO, HISTOGRAM, LINEAR_LINEAR
|
||||
use hdf5_interface
|
||||
use random_lcg, only: prn
|
||||
|
||||
!===============================================================================
|
||||
! 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(:) ! incoming energies
|
||||
type(KalbachMannTable), allocatable :: distribution(:) ! outgoing E/mu parameters
|
||||
contains
|
||||
procedure :: sample => kalbachmann_sample
|
||||
procedure :: from_hdf5 => kalbachmann_from_hdf5
|
||||
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)
|
||||
if (E_in < this%energy(1)) then
|
||||
i = 1
|
||||
r = ZERO
|
||||
elseif (E_in > this%energy(n_energy_in)) then
|
||||
i = n_energy_in - 1
|
||||
r = ONE
|
||||
else
|
||||
i = binary_search(this%energy, n_energy_in, E_in)
|
||||
r = (E_in - this%energy(i)) / &
|
||||
(this%energy(i+1) - this%energy(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%distribution(i)%e_out)
|
||||
E_i_1 = this%distribution(i)%e_out(1)
|
||||
E_i_K = this%distribution(i)%e_out(n_energy_out)
|
||||
|
||||
n_energy_out = size(this%distribution(i+1)%e_out)
|
||||
E_i1_1 = this%distribution(i+1)%e_out(1)
|
||||
E_i1_K = this%distribution(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%distribution(l)%e_out)
|
||||
r1 = prn()
|
||||
c_k = this%distribution(l)%c(1)
|
||||
do k = 1, n_energy_out - 1
|
||||
c_k1 = this%distribution(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%distribution(l)%e_out(k)
|
||||
p_l_k = this%distribution(l)%p(k)
|
||||
if (this%distribution(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%distribution(l)%r(k)
|
||||
km_a = this%distribution(l)%a(k)
|
||||
|
||||
elseif (this%distribution(l)%interpolation == LINEAR_LINEAR) then
|
||||
! Linear-linear interpolation
|
||||
E_l_k1 = this%distribution(l)%e_out(k+1)
|
||||
p_l_k1 = this%distribution(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%distribution(l)%r(k) + (E_out - E_l_k)/(E_l_k1 - E_l_k) * &
|
||||
(this%distribution(l)%r(k+1) - this%distribution(l)%r(k))
|
||||
km_a = this%distribution(l)%a(k) + (E_out - E_l_k)/(E_l_k1 - E_l_k) * &
|
||||
(this%distribution(l)%a(k+1) - this%distribution(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
|
||||
|
||||
subroutine kalbachmann_from_hdf5(this, group_id)
|
||||
class(KalbachMann), intent(inout) :: this
|
||||
integer(HID_T), intent(in) :: group_id
|
||||
|
||||
integer :: i, j, k
|
||||
integer :: n
|
||||
integer :: n_energy
|
||||
integer(HID_T) :: dset_id
|
||||
integer(HSIZE_T) :: dims(1), dims2(2)
|
||||
integer, allocatable :: temp(:,:)
|
||||
integer, allocatable :: offsets(:)
|
||||
integer, allocatable :: interp(:)
|
||||
integer, allocatable :: n_discrete(:)
|
||||
real(8), allocatable :: eout(:,:)
|
||||
|
||||
! Open incoming energy dataset
|
||||
dset_id = open_dataset(group_id, 'energy')
|
||||
|
||||
! Get interpolation parameters
|
||||
call read_attribute(temp, dset_id, 'interpolation')
|
||||
allocate(this%breakpoints(size(temp, 1)))
|
||||
allocate(this%interpolation(size(temp, 1)))
|
||||
this%breakpoints(:) = temp(:, 1)
|
||||
this%interpolation(:) = temp(:, 2)
|
||||
this%n_region = size(this%breakpoints)
|
||||
|
||||
! Get incoming energies
|
||||
call get_shape(dset_id, dims)
|
||||
n_energy = int(dims(1), 4)
|
||||
allocate(this%energy(n_energy))
|
||||
allocate(this%distribution(n_energy))
|
||||
call read_dataset(this%energy, dset_id)
|
||||
call close_dataset(dset_id)
|
||||
|
||||
! Get outgoing energy distribution data
|
||||
dset_id = open_dataset(group_id, 'distribution')
|
||||
call read_attribute(offsets, dset_id, 'offsets')
|
||||
call read_attribute(interp, dset_id, 'interpolation')
|
||||
call read_attribute(n_discrete, dset_id, 'n_discrete_lines')
|
||||
call get_shape(dset_id, dims2)
|
||||
allocate(eout(dims2(1), dims2(2)))
|
||||
call read_dataset(eout, dset_id)
|
||||
call close_dataset(dset_id)
|
||||
|
||||
do i = 1, n_energy
|
||||
! Determine number of outgoing energies
|
||||
j = offsets(i)
|
||||
if (i < n_energy) then
|
||||
n = offsets(i+1) - j
|
||||
else
|
||||
n = size(eout, 1) - j
|
||||
end if
|
||||
|
||||
associate (d => this%distribution(i))
|
||||
! Assign interpolation scheme and number of discrete lines
|
||||
d % interpolation = interp(i)
|
||||
d % n_discrete = n_discrete(i)
|
||||
|
||||
! Allocate arrays for energies and PDF/CDF
|
||||
allocate(d % e_out(n))
|
||||
allocate(d % p(n))
|
||||
allocate(d % c(n))
|
||||
allocate(d % r(n))
|
||||
allocate(d % a(n))
|
||||
|
||||
! Copy data
|
||||
d % e_out(:) = eout(j+1:j+n, 1)
|
||||
d % p(:) = eout(j+1:j+n, 2)
|
||||
d % c(:) = eout(j+1:j+n, 3)
|
||||
d % r(:) = eout(j+1:j+n, 4)
|
||||
d % a(:) = eout(j+1:j+n, 5)
|
||||
|
||||
|
||||
! To get answers that match ACE data, for now we still use the tabulated
|
||||
! CDF values that were passed through to the HDF5 library. At a later
|
||||
! time, we can remove the CDF values from the HDF5 library and
|
||||
! reconstruct them using the PDF
|
||||
if (.false.) then
|
||||
! Calculate cumulative distribution function -- discrete portion
|
||||
do k = 1, d % n_discrete
|
||||
if (k == 1) then
|
||||
d % c(k) = d % p(k)
|
||||
else
|
||||
d % c(k) = d % c(k-1) + d % p(k)
|
||||
end if
|
||||
end do
|
||||
|
||||
! Continuous portion
|
||||
do k = d % n_discrete + 1, n
|
||||
if (k == d % n_discrete + 1) then
|
||||
d % c(k) = sum(d % p(1:d % n_discrete))
|
||||
else
|
||||
if (d % interpolation == HISTOGRAM) then
|
||||
d % c(k) = d % c(k-1) + d % p(k-1) * &
|
||||
(d % e_out(k) - d % e_out(k-1))
|
||||
elseif (d % interpolation == LINEAR_LINEAR) then
|
||||
d % c(k) = d % c(k-1) + HALF*(d % p(k-1) + d % p(k)) * &
|
||||
(d % e_out(k) - d % e_out(k-1))
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
|
||||
! Normalize density and distribution functions
|
||||
d % p(:) = d % p(:)/d % c(n)
|
||||
d % c(:) = d % c(:)/d % c(n)
|
||||
end if
|
||||
end associate
|
||||
|
||||
j = j + n
|
||||
end do
|
||||
end subroutine kalbachmann_from_hdf5
|
||||
|
||||
end module secondary_kalbach
|
||||
|
|
@ -1,82 +0,0 @@
|
|||
module secondary_nbody
|
||||
|
||||
use angleenergy_header, only: AngleEnergy
|
||||
use constants, only: ONE, TWO, PI
|
||||
use hdf5_interface, only: read_attribute, HID_T
|
||||
use math, only: maxwell_spectrum
|
||||
use random_lcg, only: prn
|
||||
|
||||
!===============================================================================
|
||||
! 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(AngleEnergy) :: NBodyPhaseSpace
|
||||
integer :: n_bodies
|
||||
real(8) :: mass_ratio
|
||||
real(8) :: A
|
||||
real(8) :: Q
|
||||
contains
|
||||
procedure :: sample => nbody_sample
|
||||
procedure :: from_hdf5 => nbody_from_hdf5
|
||||
end type NBodyPhaseSpace
|
||||
|
||||
contains
|
||||
|
||||
subroutine nbody_sample(this, E_in, E_out, mu)
|
||||
class(NBodyPhaseSpace), 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 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
|
||||
|
||||
! By definition, the distribution of the angle is isotropic for an N-body
|
||||
! phase space distribution
|
||||
mu = TWO*prn() - ONE
|
||||
|
||||
! 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 subroutine nbody_sample
|
||||
|
||||
subroutine nbody_from_hdf5(this, group_id)
|
||||
class(NBodyPhaseSpace), intent(inout) :: this
|
||||
integer(HID_T), intent(in) :: group_id
|
||||
|
||||
call read_attribute(this%mass_ratio, group_id, 'total_mass')
|
||||
call read_attribute(this%n_bodies, group_id, 'n_particles')
|
||||
call read_attribute(this%A, group_id, 'atomic_weight_ratio')
|
||||
call read_attribute(this%Q, group_id, 'q_value')
|
||||
end subroutine nbody_from_hdf5
|
||||
|
||||
end module secondary_nbody
|
||||
|
|
@ -1,98 +0,0 @@
|
|||
module secondary_uncorrelated
|
||||
|
||||
use angle_distribution, only: AngleDistribution
|
||||
use angleenergy_header, only: AngleEnergy
|
||||
use constants, only: ONE, TWO, MAX_WORD_LEN
|
||||
use energy_distribution, only: EnergyDistribution, LevelInelastic, &
|
||||
ContinuousTabular, MaxwellEnergy, Evaporation, WattEnergy, DiscretePhoton
|
||||
use error, only: warning
|
||||
use hdf5_interface, only: read_attribute, open_group, close_group, &
|
||||
object_exists, HID_T
|
||||
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
|
||||
procedure :: from_hdf5 => uncorrelated_from_hdf5
|
||||
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
|
||||
|
||||
subroutine uncorrelated_from_hdf5(this, group_id)
|
||||
class(UncorrelatedAngleEnergy), intent(inout) :: this
|
||||
integer(HID_T), intent(in) :: group_id
|
||||
|
||||
integer(HID_T) :: energy_group
|
||||
integer(HID_T) :: angle_group
|
||||
character(MAX_WORD_LEN) :: type
|
||||
|
||||
! Check if angle group is present & read
|
||||
if (object_exists(group_id, 'angle')) then
|
||||
angle_group = open_group(group_id, 'angle')
|
||||
call this%angle%from_hdf5(angle_group)
|
||||
call close_group(angle_group)
|
||||
end if
|
||||
|
||||
! Check if energy group is present & read
|
||||
if (object_exists(group_id, 'energy')) then
|
||||
energy_group = open_group(group_id, 'energy')
|
||||
call read_attribute(type, energy_group, 'type')
|
||||
select case (type)
|
||||
case ('discrete_photon')
|
||||
allocate(DiscretePhoton :: this%energy)
|
||||
case ('level')
|
||||
allocate(LevelInelastic :: this%energy)
|
||||
case ('continuous')
|
||||
allocate(ContinuousTabular :: this%energy)
|
||||
case ('maxwell')
|
||||
allocate(MaxwellEnergy :: this%energy)
|
||||
case ('evaporation')
|
||||
allocate(Evaporation :: this%energy)
|
||||
case ('watt')
|
||||
allocate(WattEnergy :: this%energy)
|
||||
case default
|
||||
call warning("Energy distribution type '" // trim(type) &
|
||||
// "' not implemented.")
|
||||
end select
|
||||
|
||||
if (allocated(this % energy)) then
|
||||
call this%energy%from_hdf5(energy_group)
|
||||
end if
|
||||
|
||||
call close_group(energy_group)
|
||||
end if
|
||||
end subroutine uncorrelated_from_hdf5
|
||||
|
||||
end module secondary_uncorrelated
|
||||
Loading…
Add table
Add a link
Reference in a new issue