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 #610 from paulromano/refactor-reaction

Browse source 
Refactor of reaction product data
This commit is contained in:
Adam Nelson 2016-03-22 05:15:57 -04:00
commit 9c020cfb62
27 changed files with 2546 additions and 2708 deletions
Download patch file Download diff file Expand all files Collapse all files

721
src/ace.F90
View file

File diff suppressed because it is too large Load diff

58
src/ace_header.F90
View file

@ -1,58 +0,0 @@
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
!===============================================================================
! REACTION contains the cross-section and secondary energy and angle
! distributions for a single reaction in a continuous-energy ACE-format table
!===============================================================================
type Reaction
integer :: MT ! ENDF MT value
real(8) :: Q_value ! Reaction Q value
integer :: multiplicity ! Number of secondary particles released
type(Tab1), pointer :: multiplicity_E => null() ! Energy-dependent neutron yield
integer :: threshold ! Energy grid index of threshold
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
type(SecondaryDistribution) :: secondary
contains
procedure :: clear => reaction_clear ! Deallocates Reaction
end type Reaction
!===============================================================================
! URRDATA contains probability tables for the unresolved resonance range.
!===============================================================================
type UrrData
integer :: n_energy ! # of incident neutron energies
integer :: n_prob ! # of probabilities
integer :: interp ! inteprolation (2=lin-lin, 5=log-log)
integer :: inelastic_flag ! inelastic competition flag
integer :: absorption_flag ! other absorption flag
logical :: multiply_smooth ! multiply by smooth cross section?
real(8), allocatable :: energy(:) ! incident energies
real(8), allocatable :: prob(:,:,:) ! actual probabibility tables
end type UrrData
contains
!===============================================================================
! 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)
end subroutine reaction_clear
end module ace_header

29
src/angleenergy_header.F90 Normal file
View file

@ -0,0 +1,29 @@
module angleenergy_header
!===============================================================================
! 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(angleenergy_sample_), deferred :: sample
end type AngleEnergy
abstract interface
subroutine angleenergy_sample_(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 angleenergy_sample_
end interface
type :: AngleEnergyContainer
class(AngleEnergy), allocatable :: obj
end type AngleEnergyContainer
end module angleenergy_header

6
src/constants.F90
View file

@ -223,6 +223,12 @@ module constants
NU_POLYNOMIAL = 1, & ! Nu values given by polynomial
NU_TABULAR = 2 ! Nu values given by tabular distribution
! Secondary particle emission type
integer, parameter :: &
EMISSION_PROMPT = 1, & ! Prompt emission of secondary particle
EMISSION_DELAYED = 2, & ! Delayed emission of secondary particle
EMISSION_TOTAL = 3 ! Yield represents total emission (prompt + delayed)
! Cross section filetypes
integer, parameter :: &
ASCII = 1, & ! ASCII cross section file

220
src/cross_section.F90
View file

@ -1,10 +1,8 @@
module cross_section
use ace_header, only: Reaction, UrrData
use constants
use energy_grid, only: grid_method, log_spacing
use error, only: fatal_error
use fission, only: nu_total
use global
use list_header, only: ListElemInt
use material_header, only: Material
@ -363,131 +361,127 @@ contains
real(8) :: capture ! (n,gamma) cross section
real(8) :: fission ! fission cross section
real(8) :: inelastic ! inelastic cross section
type(UrrData), pointer :: urr
type(NuclideCE), pointer :: nuc
micro_xs(i_nuclide) % use_ptable = .true.
! get pointer to probability table
nuc => nuclides(i_nuclide)
urr => nuc % urr_data
associate (nuc => nuclides(i_nuclide), urr => nuclides(i_nuclide) % urr_data)
! determine energy table
i_energy = 1
do
if (E < urr % energy(i_energy + 1)) exit
i_energy = i_energy + 1
end do
! determine energy table
i_energy = 1
do
if (E < urr % energy(i_energy + 1)) exit
i_energy = i_energy + 1
end do
! determine interpolation factor on table
f = (E - urr % energy(i_energy)) / &
(urr % energy(i_energy + 1) - urr % energy(i_energy))
! determine interpolation factor on table
f = (E - urr % energy(i_energy)) / &
(urr % energy(i_energy + 1) - urr % energy(i_energy))
! sample probability table using the cumulative distribution
! sample probability table using the cumulative distribution
! Random numbers for xs calculation are sampled from a separated stream.
! This guarantees the randomness and, at the same time, makes sure we reuse
! random number for the same nuclide at different temperatures, therefore
! preserving correlation of temperature in probability tables.
call prn_set_stream(STREAM_URR_PTABLE)
r = future_prn(int(nuc_zaid_dict % get_key(nuc % zaid), 8))
call prn_set_stream(STREAM_TRACKING)
! Random numbers for xs calculation are sampled from a separated stream.
! This guarantees the randomness and, at the same time, makes sure we reuse
! random number for the same nuclide at different temperatures, therefore
! preserving correlation of temperature in probability tables.
call prn_set_stream(STREAM_URR_PTABLE)
r = future_prn(int(nuc_zaid_dict % get_key(nuc % zaid), 8))
call prn_set_stream(STREAM_TRACKING)
i_low = 1
do
if (urr % prob(i_energy, URR_CUM_PROB, i_low) > r) exit
i_low = i_low + 1
end do
i_up = 1
do
if (urr % prob(i_energy + 1, URR_CUM_PROB, i_up) > r) exit
i_up = i_up + 1
end do
i_low = 1
do
if (urr % prob(i_energy, URR_CUM_PROB, i_low) > r) exit
i_low = i_low + 1
end do
i_up = 1
do
if (urr % prob(i_energy + 1, URR_CUM_PROB, i_up) > r) exit
i_up = i_up + 1
end do
! determine elastic, fission, and capture cross sections from probability
! table
if (urr % interp == LINEAR_LINEAR) then
elastic = (ONE - f) * urr % prob(i_energy, URR_ELASTIC, i_low) + &
f * urr % prob(i_energy + 1, URR_ELASTIC, i_up)
fission = (ONE - f) * urr % prob(i_energy, URR_FISSION, i_low) + &
f * urr % prob(i_energy + 1, URR_FISSION, i_up)
capture = (ONE - f) * urr % prob(i_energy, URR_N_GAMMA, i_low) + &
f * urr % prob(i_energy + 1, URR_N_GAMMA, i_up)
elseif (urr % interp == LOG_LOG) then
! Get logarithmic interpolation factor
f = log(E / urr % energy(i_energy)) / &
log(urr % energy(i_energy + 1) / urr % energy(i_energy))
! determine elastic, fission, and capture cross sections from probability
! table
if (urr % interp == LINEAR_LINEAR) then
elastic = (ONE - f) * urr % prob(i_energy, URR_ELASTIC, i_low) + &
f * urr % prob(i_energy + 1, URR_ELASTIC, i_up)
fission = (ONE - f) * urr % prob(i_energy, URR_FISSION, i_low) + &
f * urr % prob(i_energy + 1, URR_FISSION, i_up)
capture = (ONE - f) * urr % prob(i_energy, URR_N_GAMMA, i_low) + &
f * urr % prob(i_energy + 1, URR_N_GAMMA, i_up)
elseif (urr % interp == LOG_LOG) then
! Get logarithmic interpolation factor
f = log(E / urr % energy(i_energy)) / &
log(urr % energy(i_energy + 1) / urr % energy(i_energy))
! Calculate elastic cross section/factor
elastic = ZERO
if (urr % prob(i_energy, URR_ELASTIC, i_low) > ZERO .and. &
urr % prob(i_energy + 1, URR_ELASTIC, i_up) > ZERO) then
elastic = exp((ONE - f) * log(urr % prob(i_energy, URR_ELASTIC, &
i_low)) + f * log(urr % prob(i_energy + 1, URR_ELASTIC, &
i_up)))
end if
! Calculate fission cross section/factor
fission = ZERO
if (urr % prob(i_energy, URR_FISSION, i_low) > ZERO .and. &
urr % prob(i_energy + 1, URR_FISSION, i_up) > ZERO) then
fission = exp((ONE - f) * log(urr % prob(i_energy, URR_FISSION, &
i_low)) + f * log(urr % prob(i_energy + 1, URR_FISSION, &
i_up)))
end if
! Calculate capture cross section/factor
capture = ZERO
if (urr % prob(i_energy, URR_N_GAMMA, i_low) > ZERO .and. &
urr % prob(i_energy + 1, URR_N_GAMMA, i_up) > ZERO) then
capture = exp((ONE - f) * log(urr % prob(i_energy, URR_N_GAMMA, &
i_low)) + f * log(urr % prob(i_energy + 1, URR_N_GAMMA, &
i_up)))
end if
end if
! Determine treatment of inelastic scattering
inelastic = ZERO
if (urr % inelastic_flag > 0) then
! Get index on energy grid and interpolation factor
i_energy = micro_xs(i_nuclide) % index_grid
f = micro_xs(i_nuclide) % interp_factor
! Determine inelastic scattering cross section
associate (rxn => nuc % reactions(nuc % urr_inelastic))
if (i_energy >= rxn % threshold) then
inelastic = (ONE - f) * rxn % sigma(i_energy - rxn%threshold + 1) + &
f * rxn % sigma(i_energy - rxn%threshold + 2)
! Calculate elastic cross section/factor
elastic = ZERO
if (urr % prob(i_energy, URR_ELASTIC, i_low) > ZERO .and. &
urr % prob(i_energy + 1, URR_ELASTIC, i_up) > ZERO) then
elastic = exp((ONE - f) * log(urr % prob(i_energy, URR_ELASTIC, &
i_low)) + f * log(urr % prob(i_energy + 1, URR_ELASTIC, &
i_up)))
end if
end associate
end if
! Multiply by smooth cross-section if needed
if (urr % multiply_smooth) then
elastic = elastic * micro_xs(i_nuclide) % elastic
capture = capture * (micro_xs(i_nuclide) % absorption - &
micro_xs(i_nuclide) % fission)
fission = fission * micro_xs(i_nuclide) % fission
end if
! Calculate fission cross section/factor
fission = ZERO
if (urr % prob(i_energy, URR_FISSION, i_low) > ZERO .and. &
urr % prob(i_energy + 1, URR_FISSION, i_up) > ZERO) then
fission = exp((ONE - f) * log(urr % prob(i_energy, URR_FISSION, &
i_low)) + f * log(urr % prob(i_energy + 1, URR_FISSION, &
i_up)))
end if
! Check for negative values
if (elastic < ZERO) elastic = ZERO
if (fission < ZERO) fission = ZERO
if (capture < ZERO) capture = ZERO
! Calculate capture cross section/factor
capture = ZERO
if (urr % prob(i_energy, URR_N_GAMMA, i_low) > ZERO .and. &
urr % prob(i_energy + 1, URR_N_GAMMA, i_up) > ZERO) then
capture = exp((ONE - f) * log(urr % prob(i_energy, URR_N_GAMMA, &
i_low)) + f * log(urr % prob(i_energy + 1, URR_N_GAMMA, &
i_up)))
end if
end if
! Set elastic, absorption, fission, and total cross sections. Note that the
! total cross section is calculated as sum of partials rather than using the
! table-provided value
micro_xs(i_nuclide) % elastic = elastic
micro_xs(i_nuclide) % absorption = capture + fission
micro_xs(i_nuclide) % fission = fission
micro_xs(i_nuclide) % total = elastic + inelastic + capture + fission
! Determine treatment of inelastic scattering
inelastic = ZERO
if (urr % inelastic_flag > 0) then
! Get index on energy grid and interpolation factor
i_energy = micro_xs(i_nuclide) % index_grid
f = micro_xs(i_nuclide) % interp_factor
! Determine nu-fission cross section
if (nuc % fissionable) then
micro_xs(i_nuclide) % nu_fission = nu_total(nuc, E) * &
micro_xs(i_nuclide) % fission
end if
! Determine inelastic scattering cross section
associate (rxn => nuc % reactions(nuc % urr_inelastic))
if (i_energy >= rxn % threshold) then
inelastic = (ONE - f) * rxn % sigma(i_energy - rxn%threshold + 1) + &
f * rxn % sigma(i_energy - rxn%threshold + 2)
end if
end associate
end if
! Multiply by smooth cross-section if needed
if (urr % multiply_smooth) then
elastic = elastic * micro_xs(i_nuclide) % elastic
capture = capture * (micro_xs(i_nuclide) % absorption - &
micro_xs(i_nuclide) % fission)
fission = fission * micro_xs(i_nuclide) % fission
end if
! Check for negative values
if (elastic < ZERO) elastic = ZERO
if (fission < ZERO) fission = ZERO
if (capture < ZERO) capture = ZERO
! Set elastic, absorption, fission, and total cross sections. Note that the
! total cross section is calculated as sum of partials rather than using the
! table-provided value
micro_xs(i_nuclide) % elastic = elastic
micro_xs(i_nuclide) % absorption = capture + fission
micro_xs(i_nuclide) % fission = fission
micro_xs(i_nuclide) % total = elastic + inelastic + capture + fission
! Determine nu-fission cross section
if (nuc % fissionable) then
micro_xs(i_nuclide) % nu_fission = nuc % nu(E, EMISSION_TOTAL) * &
micro_xs(i_nuclide) % fission
end if
end associate
end subroutine calculate_urr_xs

190
src/endf_header.F90
View file

@ -1,12 +1,51 @@
module endf_header
implicit none
use constants, only: ZERO, HISTOGRAM, LINEAR_LINEAR, LINEAR_LOG, &
LOG_LINEAR, LOG_LOG
use search, only: binary_search
implicit none
type, abstract :: Function1D
contains
procedure(function1d_evaluate_), deferred :: evaluate
end type Function1D
abstract interface
pure function function1d_evaluate_(this, x) result(y)
import Function1D
class(Function1D), intent(in) :: this
real(8), intent(in) :: x
real(8) :: y
end function function1d_evaluate_
end interface
!===============================================================================
! TAB1 represents a one-dimensional interpolable function
! CONSTANT1D represents a constant one-dimensional function
!===============================================================================
type Tab1
type, extends(Function1D) :: Constant1D
real(8) :: y
contains
procedure :: evaluate => constant1d_evaluate
end type Constant1D
!===============================================================================
! POLYNOMIAL represents a one-dimensional function expressed as a polynomial
!===============================================================================
type, extends(Function1D) :: Polynomial
real(8), allocatable :: coef(:) ! coefficients
contains
procedure :: evaluate => polynomial_evaluate
procedure :: from_ace => polynomial_from_ace
end type Polynomial
!===============================================================================
! TABULATED1D represents a one-dimensional interpolable function
!===============================================================================
type, extends(Function1D) :: Tabulated1D
integer :: n_regions = 0 ! # of interpolation regions
integer, allocatable :: nbt(:) ! values separating interpolation regions
integer, allocatable :: int(:) ! interpolation scheme
@ -14,18 +53,78 @@ module endf_header
real(8), allocatable :: x(:) ! values of abscissa
real(8), allocatable :: y(:) ! values of ordinate
contains
procedure :: from_ace
end type Tab1
procedure :: from_ace => tabulated1d_from_ace
procedure :: evaluate => tabulated1d_evaluate
end type Tabulated1D
contains
subroutine from_ace(this, xss, idx)
class(Tab1), intent(inout) :: this
!===============================================================================
! Constant1D implementation
!===============================================================================
pure function constant1d_evaluate(this, x) result(y)
class(Constant1D), intent(in) :: this
real(8), intent(in) :: x
real(8) :: y
y = this % y
end function constant1d_evaluate
!===============================================================================
! Polynomial implementation
!===============================================================================
subroutine polynomial_from_ace(this, xss, idx)
class(Polynomial), intent(inout) :: this
real(8), intent(in) :: xss(:)
integer, intent(in) :: idx
integer :: nc ! number of coefficients (order - 1)
! Clear space
if (allocated(this % coef)) deallocate(this % coef)
! Determine number of coefficients
nc = nint(xss(idx))
! Allocate space for and read coefficients
allocate(this % coef(nc))
this % coef(:) = xss(idx + 1 : idx + nc)
end subroutine polynomial_from_ace
pure function polynomial_evaluate(this, x) result(y)
class(Polynomial), intent(in) :: this
real(8), intent(in) :: x
real(8) :: y
integer :: i
! Use Horner's rule to evaluate polynomial. Note that coefficients are
! ordered in increasing powers of x.
y = ZERO
do i = size(this % coef), 1, -1
y = y*x + this % coef(i)
end do
end function polynomial_evaluate
!===============================================================================
! Tabulated1D implementation
!===============================================================================
subroutine tabulated1d_from_ace(this, xss, idx)
class(Tabulated1D), intent(inout) :: this
real(8), intent(in) :: xss(:)
integer, intent(in) :: idx
integer :: nr, ne
! Clear space
if (allocated(this % nbt)) deallocate(this % nbt)
if (allocated(this % int)) deallocate(this % int)
if (allocated(this % x)) deallocate(this % x)
if (allocated(this % y)) deallocate(this % y)
! Determine number of regions
nr = nint(xss(idx))
this%n_regions = nr
@ -47,6 +146,81 @@ contains
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 subroutine tabulated1d_from_ace
pure function tabulated1d_evaluate(this, x) result(y)
class(Tabulated1D), intent(in) :: this
real(8), intent(in) :: x ! x value to find y at
real(8) :: y ! y(x)
integer :: i ! bin in which to interpolate
integer :: j ! index for interpolation region
integer :: n_regions ! number of interpolation regions
integer :: n_pairs ! number of tabulated values
integer :: interp ! ENDF interpolation scheme
real(8) :: r ! interpolation factor
real(8) :: x0, x1 ! bounding x values
real(8) :: y0, y1 ! bounding y values
! determine number of interpolation regions and pairs
n_regions = this % n_regions
n_pairs = this % n_pairs
! find which bin the abscissa is in -- if the abscissa is outside the
! tabulated range, the first or last point is chosen, i.e. no interpolation
! is done outside the energy range
if (x < this % x(1)) then
y = this % y(1)
return
elseif (x > this % x(n_pairs)) then
y = this % y(n_pairs)
return
else
i = binary_search(this % x, n_pairs, x)
end if
! determine interpolation scheme
if (n_regions == 0) then
interp = LINEAR_LINEAR
elseif (n_regions == 1) then
interp = this % int(1)
elseif (n_regions > 1) then
do j = 1, n_regions
if (i < this % nbt(j)) then
interp = this % int(j)
exit
end if
end do
end if
! handle special case of histogram interpolation
if (interp == HISTOGRAM) then
y = this % y(i)
return
end if
! determine bounding values
x0 = this % x(i)
x1 = this % x(i + 1)
y0 = this % y(i)
y1 = this % y(i + 1)
! determine interpolation factor and interpolated value
select case (interp)
case (LINEAR_LINEAR)
r = (x - x0)/(x1 - x0)
y = y0 + r*(y1 - y0)
case (LINEAR_LOG)
r = log(x/x0)/log(x1/x0)
y = y0 + r*(y1 - y0)
case (LOG_LINEAR)
r = (x - x0)/(x1 - x0)
y = y0*exp(r*log(y1/y0))
case (LOG_LOG)
r = log(x/x0)/log(x1/x0)
y = y0*exp(r*log(y1/y0))
end select
end function tabulated1d_evaluate
end module endf_header

125
src/energy_distribution.F90
View file

@ -1,8 +1,7 @@
module energy_distribution
use constants, only: ZERO, ONE, TWO, PI, HISTOGRAM, LINEAR_LINEAR
use endf_header, only: Tab1
use interpolation, only: interpolate_tab1
use endf_header, only: Tabulated1D
use math, only: maxwell_spectrum, watt_spectrum
use random_lcg, only: prn
use search, only: binary_search
@ -83,8 +82,8 @@ module energy_distribution
integer :: n_region
integer, allocatable :: breakpoints(:)
integer, allocatable :: interpolation(:)
real(8), allocatable :: energy_in(:)
type(CTTable), allocatable :: energy_out(:)
real(8), allocatable :: energy(:)
type(CTTable), allocatable :: distribution(:)
contains
procedure :: sample => continuous_sample
end type ContinuousTabular
@ -95,7 +94,7 @@ module energy_distribution
!===============================================================================
type, extends(EnergyDistribution) :: MaxwellEnergy
type(Tab1) :: theta ! incoming-energy-dependent parameter
type(Tabulated1D) :: theta ! incoming-energy-dependent parameter
real(8) :: u ! restriction energy
contains
procedure :: sample => maxwellenergy_sample
@ -107,7 +106,7 @@ module energy_distribution
!===============================================================================
type, extends(EnergyDistribution) :: Evaporation
type(Tab1) :: theta
type(Tabulated1D) :: theta
real(8) :: u
contains
procedure :: sample => evaporation_sample
@ -119,28 +118,13 @@ module energy_distribution
!===============================================================================
type, extends(EnergyDistribution) :: WattEnergy
type(Tab1) :: a
type(Tab1) :: b
type(Tabulated1D) :: a
type(Tabulated1D) :: 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)
@ -202,6 +186,7 @@ contains
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
@ -210,6 +195,7 @@ contains
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
@ -238,17 +224,17 @@ contains
! 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
n_energy_in = size(this%energy)
if (E_in < this%energy(1)) then
i = 1
r = ZERO
elseif (E_in > this%energy_in(n_energy_in)) then
elseif (E_in > this%energy(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))
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
@ -263,23 +249,23 @@ contains
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%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%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)
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%energy_out(l)%e_out)
n_energy_out = size(this%distribution(l)%e_out)
r1 = prn()
c_k = this%energy_out(l)%c(1)
c_k = this%distribution(l)%c(1)
do k = 1, n_energy_out - 1
c_k1 = this%energy_out(l)%c(k+1)
c_k1 = this%distribution(l)%c(k+1)
if (r1 < c_k1) exit
c_k = c_k1
end do
@ -287,9 +273,9 @@ contains
! 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
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
@ -297,10 +283,10 @@ contains
E_out = E_l_k
end if
elseif (this%energy_out(l)%interpolation == LINEAR_LINEAR) then
elseif (this%distribution(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)
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
@ -321,6 +307,7 @@ contains
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
@ -329,7 +316,7 @@ contains
real(8) :: theta ! Maxwell distribution parameter
! Get temperature corresponding to incoming energy
theta = interpolate_tab1(this%theta, E_in)
theta = this % theta % evaluate(E_in)
do
! Sample maxwell fission spectrum
@ -349,9 +336,9 @@ contains
real(8) :: x, y, v
! Get temperature corresponding to incoming energy
theta = interpolate_tab1(this%theta, E_in)
theta = this % theta % evaluate(E_in)
y = (E_in - this%U)/theta
y = (E_in - this%u)/theta
v = 1 - exp(-y)
! Sample outgoing energy based on evaporation spectrum probability
@ -372,10 +359,10 @@ contains
real(8) :: a, b ! Watt spectrum parameters
! Determine Watt parameter 'a' from tabulated function
a = interpolate_tab1(this%a, E_in)
a = this % a % evaluate(E_in)
! Determine Watt parameter 'b' from tabulated function
b = interpolate_tab1(this%b, E_in)
b = this % b % evaluate(E_in)
do
! Sample energy-dependent Watt fission spectrum
@ -386,44 +373,4 @@ contains
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

161
src/fission.F90
View file

@ -1,161 +0,0 @@
module fission
use nuclide_header, only: NuclideCE
use constants
use error, only: fatal_error
use interpolation, only: interpolate_tab1
use search, only: binary_search
implicit none
contains
!===============================================================================
! NU_TOTAL calculates the total number of neutrons emitted per fission for a
! given nuclide and incoming neutron energy
!===============================================================================
pure function nu_total(nuc, E) result(nu)
type(NuclideCE), intent(in) :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: nu ! number of total neutrons emitted per fission
integer :: i ! loop index
integer :: NC ! number of polynomial coefficients
real(8) :: c ! polynomial coefficient
if (nuc % nu_t_type == NU_NONE) then
nu = ERROR_REAL
elseif (nuc % nu_t_type == NU_POLYNOMIAL) then
! determine number of coefficients
NC = int(nuc % nu_t_data(1))
! sum up polynomial in energy
nu = ZERO
do i = 0, NC - 1
c = nuc % nu_t_data(i+2)
nu = nu + c * E**i
end do
elseif (nuc % nu_t_type == NU_TABULAR) then
! use ENDF interpolation laws to determine nu
nu = interpolate_tab1(nuc % nu_t_data, E)
end if
end function nu_total
!===============================================================================
! NU_PROMPT calculates the total number of prompt neutrons emitted per fission
! for a given nuclide and incoming neutron energy
!===============================================================================
pure function nu_prompt(nuc, E) result(nu)
type(NuclideCE), intent(in) :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: nu ! number of prompt neutrons emitted per fission
integer :: i ! loop index
integer :: NC ! number of polynomial coefficients
real(8) :: c ! polynomial coefficient
if (nuc % nu_p_type == NU_NONE) then
! since no prompt or delayed data is present, this means all neutron
! emission is prompt -- WARNING: This currently returns zero. The calling
! routine needs to know this situation is occurring since we don't want
! to call nu_total unnecessarily if it has already been called.
nu = ZERO
elseif (nuc % nu_p_type == NU_POLYNOMIAL) then
! determine number of coefficients
NC = int(nuc % nu_p_data(1))
! sum up polynomial in energy
nu = ZERO
do i = 0, NC - 1
c = nuc % nu_p_data(i+2)
nu = nu + c * E**i
end do
elseif (nuc % nu_p_type == NU_TABULAR) then
! use ENDF interpolation laws to determine nu
nu = interpolate_tab1(nuc % nu_p_data, E)
end if
end function nu_prompt
!===============================================================================
! NU_DELAYED calculates the total number of delayed neutrons emitted per fission
! for a given nuclide and incoming neutron energy
!===============================================================================
pure function nu_delayed(nuc, E) result(nu)
type(NuclideCE), intent(in) :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: nu ! number of delayed neutrons emitted per fission
if (nuc % nu_d_type == NU_NONE) then
! since no prompt or delayed data is present, this means all neutron
! emission is prompt -- WARNING: This currently returns zero. The calling
! routine needs to know this situation is occurring since we don't want
! to call nu_delayed unnecessarily if it has already been called.
nu = ZERO
elseif (nuc % nu_d_type == NU_TABULAR) then
! use ENDF interpolation laws to determine nu
nu = interpolate_tab1(nuc % nu_d_data, E)
end if
end function nu_delayed
!===============================================================================
! YIELD_DELAYED calculates the fractional yield of delayed neutrons emitted for
! a given nuclide and incoming neutron energy in a given delayed group.
!===============================================================================
pure function yield_delayed(nuc, E, g) result(yield)
type(NuclideCE), intent(in) :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: yield ! delayed neutron precursor yield
integer, intent(in) :: g ! the delayed neutron precursor group
integer :: d ! precursor group
integer :: lc ! index before start of energies/nu values
integer :: NR ! number of interpolation regions
integer :: NE ! number of energies tabulated
yield = ZERO
if (g > nuc % n_precursor .or. g < 1) then
! if the precursor group is outside the range of precursor groups for
! the input nuclide, return ZERO.
yield = ZERO
else if (nuc % nu_d_type == NU_NONE) then
! since no prompt or delayed data is present, this means all neutron
! emission is prompt -- WARNING: This currently returns zero. The calling
! routine needs to know this situation is occurring since we don't want
! to call yield_delayed unnecessarily if it has already been called.
yield = ZERO
else if (nuc % nu_d_type == NU_TABULAR) then
lc = 1
! loop over delayed groups and determine the yield for the desired group
do d = 1, nuc % n_precursor
! determine number of interpolation regions and energies
NR = int(nuc % nu_d_precursor_data(lc + 1))
NE = int(nuc % nu_d_precursor_data(lc + 2 + 2*NR))
! check if this is the desired group
if (d == g) then
! determine delayed neutron precursor yield for group g
yield = interpolate_tab1(nuc % nu_d_precursor_data( &
lc+1:lc+2+2*NR+2*NE), E)
exit
end if
! advance pointer
lc = lc + 2 + 2*NR + 2*NE + 1
end do
end if
end function yield_delayed
end module fission

205
src/interpolation.F90
View file

@ -1,205 +0,0 @@
module interpolation
use constants
use endf_header, only: Tab1
use search, only: binary_search
use string, only: to_str
implicit none
interface interpolate_tab1
module procedure interpolate_tab1_array, interpolate_tab1_object
end interface interpolate_tab1
contains
!===============================================================================
! INTERPOLATE_TAB1_ARRAY interpolates a function between two points based on
! particular interpolation scheme. The data needs to be organized as a ENDF TAB1
! type function containing the interpolation regions, break points, and
! tabulated x's and y's.
!===============================================================================
pure function interpolate_tab1_array(data, x, loc_start) result(y)
real(8), intent(in) :: data(:) ! array of data
real(8), intent(in) :: x ! x value to find y at
integer, intent(in), optional :: loc_start ! starting location in data
real(8) :: y ! y(x)
integer :: i ! bin in which to interpolate
integer :: j ! index for interpolation region
integer :: loc_0 ! starting location
integer :: n_regions ! number of interpolation regions
integer :: n_points ! number of tabulated values
integer :: interp ! ENDF interpolation scheme
integer :: loc_breakpoints ! location of breakpoints in data
integer :: loc_interp ! location of interpolation schemes in data
integer :: loc_x ! location of x's in data
integer :: loc_y ! location of y's in data
real(8) :: r ! interpolation factor
real(8) :: x0, x1 ! bounding x values
real(8) :: y0, y1 ! bounding y values
! determine starting location
if (present(loc_start)) then
loc_0 = loc_start - 1
else
loc_0 = 0
end if
! determine number of interpolation regions
n_regions = int(data(loc_0 + 1))
! set locations for breakpoints and interpolation schemes
loc_breakpoints = loc_0 + 1
loc_interp = loc_breakpoints + n_regions
! determine number of tabulated values
n_points = int(data(loc_interp + n_regions + 1))
! set locations for x's and y's
loc_x = loc_interp + n_regions + 1
loc_y = loc_x + n_points
! find which bin the abscissa is in -- if the abscissa is outside the
! tabulated range, the first or last point is chosen, i.e. no interpolation
! is done outside the energy range
if (x < data(loc_x + 1)) then
y = data(loc_y + 1)
return
elseif (x > data(loc_x + n_points)) then
y = data(loc_y + n_points)
return
else
i = binary_search(data(loc_x + 1:loc_x + n_points), n_points, x)
end if
! determine interpolation scheme
if (n_regions == 0) then
interp = LINEAR_LINEAR
elseif (n_regions == 1) then
interp = int(data(loc_interp + 1))
elseif (n_regions > 1) then
do j = 1, n_regions
if (i < data(loc_breakpoints + j)) then
interp = int(data(loc_interp + j))
exit
end if
end do
end if
! handle special case of histogram interpolation
if (interp == HISTOGRAM) then
y = data(loc_y + i)
return
end if
! determine bounding values
x0 = data(loc_x + i)
x1 = data(loc_x + i + 1)
y0 = data(loc_y + i)
y1 = data(loc_y + i + 1)
! determine interpolation factor and interpolated value
select case (interp)
case (LINEAR_LINEAR)
r = (x - x0)/(x1 - x0)
y = y0 + r*(y1 - y0)
case (LINEAR_LOG)
r = log(x/x0)/log(x1/x0)
y = y0 + r*(y1 - y0)
case (LOG_LINEAR)
r = (x - x0)/(x1 - x0)
y = y0*exp(r*log(y1/y0))
case (LOG_LOG)
r = log(x/x0)/log(x1/x0)
y = y0*exp(r*log(y1/y0))
end select
end function interpolate_tab1_array
!===============================================================================
! INTERPOLATE_TAB1_OBJECT interpolates a function between two points based on
! particular interpolation scheme. The data needs to be organized as a ENDF TAB1
! type function containing the interpolation regions, break points, and
! tabulated x's and y's.
!===============================================================================
pure function interpolate_tab1_object(obj, x) result(y)
type(Tab1), intent(in) :: obj ! ENDF Tab1 interpolable function
real(8), intent(in) :: x ! x value to find y at
real(8) :: y ! y(x)
integer :: i ! bin in which to interpolate
integer :: j ! index for interpolation region
integer :: n_regions ! number of interpolation regions
integer :: n_pairs ! number of tabulated values
integer :: interp ! ENDF interpolation scheme
real(8) :: r ! interpolation factor
real(8) :: x0, x1 ! bounding x values
real(8) :: y0, y1 ! bounding y values
! determine number of interpolation regions and pairs
n_regions = obj % n_regions
n_pairs = obj % n_pairs
! find which bin the abscissa is in -- if the abscissa is outside the
! tabulated range, the first or last point is chosen, i.e. no interpolation
! is done outside the energy range
if (x < obj % x(1)) then
y = obj % y(1)
return
elseif (x > obj % x(n_pairs)) then
y = obj % y(n_pairs)
return
else
i = binary_search(obj % x, n_pairs, x)
end if
! determine interpolation scheme
if (n_regions == 0) then
interp = LINEAR_LINEAR
elseif (n_regions == 1) then
interp = obj % int(1)
elseif (n_regions > 1) then
do j = 1, n_regions
if (i < obj % nbt(j)) then
interp = obj % int(j)
exit
end if
end do
end if
! handle special case of histogram interpolation
if (interp == HISTOGRAM) then
y = obj % y(i)
return
end if
! determine bounding values
x0 = obj % x(i)
x1 = obj % x(i + 1)
y0 = obj % y(i)
y1 = obj % y(i + 1)
! determine interpolation factor and interpolated value
select case (interp)
case (LINEAR_LINEAR)
r = (x - x0)/(x1 - x0)
y = y0 + r*(y1 - y0)
case (LINEAR_LOG)
r = log(x/x0)/log(x1/x0)
y = y0 + r*(y1 - y0)
case (LOG_LINEAR)
r = (x - x0)/(x1 - x0)
y = y0*exp(r*log(y1/y0))
case (LOG_LOG)
r = log(x/x0)/log(x1/x0)
y = y0*exp(r*log(y1/y0))
end select
end function interpolate_tab1_object
end module interpolation

212
src/nuclide_header.F90
View file

@ -2,13 +2,18 @@ module nuclide_header
use, intrinsic :: ISO_FORTRAN_ENV
use ace_header
use constants
use endf, only: reaction_name
use error, only: fatal_error
use dict_header, only: DictIntInt
use endf, only: reaction_name, is_fission, is_disappearance
use endf_header, only: Function1D
use error, only: fatal_error, warning
use list_header, only: ListInt
use math, only: evaluate_legendre, find_angle
use product_header, only: AngleEnergyContainer
use reaction_header, only: Reaction
use stl_vector, only: VectorInt
use string
use urr_header, only: UrrData
use xml_interface
implicit none
@ -67,24 +72,11 @@ module nuclide_header
real(8) :: E_max ! upper cutoff energy for res scattering
! Fission information
logical :: has_partial_fission ! nuclide has partial fission reactions?
integer :: n_fission ! # of fission reactions
logical :: has_partial_fission = .false. ! nuclide has partial fission reactions?
integer :: n_fission ! # of fission reactions
integer :: n_precursor = 0 ! # of delayed neutron precursors
integer, allocatable :: index_fission(:) ! indices in reactions
! Total fission neutron emission
integer :: nu_t_type
real(8), allocatable :: nu_t_data(:)
! Prompt fission neutron emission
integer :: nu_p_type
real(8), allocatable :: nu_p_data(:)
! Delayed fission neutron emission
integer :: nu_d_type
integer :: n_precursor ! # of delayed neutron precursors
real(8), allocatable :: nu_d_data(:)
real(8), allocatable :: nu_d_precursor_data(:)
type(AngleEnergyContainer), allocatable :: nu_d_edist(:)
class(Function1D), allocatable :: total_nu
! Unresolved resonance data
logical :: urr_present
@ -100,6 +92,7 @@ module nuclide_header
contains
procedure :: clear => nuclidece_clear
procedure :: print => nuclidece_print
procedure :: nu => nuclidece_nu
end type NuclideCE
type, abstract, extends(Nuclide) :: NuclideMG
@ -680,81 +673,133 @@ module nuclide_header
! or NuclideAngle
!===============================================================================
subroutine nuclidece_clear(this)
subroutine nuclidece_clear(this)
class(NuclideCE), intent(inout) :: this ! The Nuclide object to clear
class(NuclideCE), intent(inout) :: this ! The Nuclide object to clear
integer :: i ! Loop counter
if (associated(this % urr_data)) deallocate(this % urr_data)
if (associated(this % urr_data)) deallocate(this % urr_data)
call this % reaction_index % clear()
if (allocated(this % reactions)) then
do i = 1, size(this % reactions)
call this % reactions(i) % clear()
end do
end subroutine nuclidece_clear
function nuclidece_nu(this, E, emission_mode, group) result(nu)
class(NuclideCE), intent(in) :: this
real(8), intent(in) :: E
integer, intent(in) :: emission_mode
integer, optional, intent(in) :: group
real(8) :: nu
integer :: i
if (.not. this % fissionable) then
nu = ZERO
return
end if
select case (emission_mode)
case (EMISSION_PROMPT)
associate (product => this % reactions(this % index_fission(1)) % products(1))
nu = product % yield % evaluate(E)
end associate
case (EMISSION_DELAYED)
if (this % n_precursor > 0) then
if (present(group)) then
! If delayed group specified, determine yield immediately
associate(p => this % reactions(this % index_fission(1)) % products(1 + group))
nu = p % yield % evaluate(E)
end associate
else
nu = ZERO
associate (rx => this % reactions(this % index_fission(1)))
do i = 2, size(rx % products)
associate (product => rx % products(i))
! Skip any non-neutron products
if (product % particle /= NEUTRON) exit
! Evaluate yield
if (product % emission_mode == EMISSION_DELAYED) then
nu = nu + product % yield % evaluate(E)
end if
end associate
end do
end associate
end if
else
nu = ZERO
end if
call this % reaction_index % clear()
case (EMISSION_TOTAL)
if (allocated(this % total_nu)) then
nu = this % total_nu % evaluate(E)
else
associate (rx => this % reactions(this % index_fission(1)))
nu = rx % products(1) % yield % evaluate(E)
end associate
end if
end select
end subroutine nuclidece_clear
end function nuclidece_nu
!===============================================================================
! NUCLIDE*_PRINT displays information about a continuous-energy neutron
! cross_section table and its reactions and secondary angle/energy distributions
!===============================================================================
subroutine nuclidece_print(this, unit)
class(NuclideCE), intent(in) :: this
integer, intent(in), optional :: unit
subroutine nuclidece_print(this, unit)
class(NuclideCE), intent(in) :: this
integer, intent(in), optional :: unit
integer :: i ! loop index over nuclides
integer :: unit_ ! unit to write to
integer :: size_xs ! memory used for cross-sections (bytes)
integer :: size_urr ! memory used for probability tables (bytes)
type(UrrData), pointer :: urr
integer :: i ! loop index over nuclides
integer :: unit_ ! unit to write to
integer :: size_xs ! memory used for cross-sections (bytes)
integer :: size_urr ! memory used for probability tables (bytes)
! set default unit for writing information
if (present(unit)) then
unit_ = unit
else
unit_ = OUTPUT_UNIT
end if
! set default unit for writing information
if (present(unit)) then
unit_ = unit
else
unit_ = OUTPUT_UNIT
end if
! Initialize totals
size_urr = 0
size_xs = 0
! Initialize totals
size_urr = 0
size_xs = 0
! Basic nuclide information
write(unit_,*) 'Nuclide ' // trim(this % name)
write(unit_,*) ' zaid = ' // trim(to_str(this % zaid))
write(unit_,*) ' awr = ' // trim(to_str(this % awr))
write(unit_,*) ' kT = ' // trim(to_str(this % kT))
write(unit_,*) ' # of grid points = ' // trim(to_str(this % n_grid))
write(unit_,*) ' Fissionable = ', this % fissionable
write(unit_,*) ' # of fission reactions = ' // trim(to_str(this % n_fission))
write(unit_,*) ' # of reactions = ' // trim(to_str(this % n_reaction))
! Basic nuclide information
write(unit_,*) 'Nuclide ' // trim(this % name)
write(unit_,*) ' zaid = ' // trim(to_str(this % zaid))
write(unit_,*) ' awr = ' // trim(to_str(this % awr))
write(unit_,*) ' kT = ' // trim(to_str(this % kT))
write(unit_,*) ' # of grid points = ' // trim(to_str(this % n_grid))
write(unit_,*) ' Fissionable = ', this % fissionable
write(unit_,*) ' # of fission reactions = ' // trim(to_str(this % n_fission))
write(unit_,*) ' # of reactions = ' // trim(to_str(this % n_reaction))
! Information on each reaction
write(unit_,*) ' Reaction Q-value COM IE'
do i = 1, this % n_reaction
associate (rxn => this % reactions(i))
write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,I6)') &
reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, &
rxn % threshold
! Information on each reaction
write(unit_,*) ' Reaction Q-value COM IE'
do i = 1, this % n_reaction
associate (rxn => this % reactions(i))
write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,I6)') &
reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, &
rxn % threshold
! Accumulate data size
size_xs = size_xs + (this % n_grid - rxn%threshold + 1) * 8
end associate
end do
! Accumulate data size
size_xs = size_xs + (this % n_grid - rxn%threshold + 1) * 8
end associate
end do
! Add memory required for summary reactions (total, absorption, fission,
! nu-fission)
size_xs = 8 * this % n_grid * 4
! Add memory required for summary reactions (total, absorption, fission,
! nu-fission)
size_xs = 8 * this % n_grid * 4
! Write information about URR probability tables
size_urr = 0
if (this % urr_present) then
urr => this % urr_data
! Write information about URR probability tables
size_urr = 0
if (this % urr_present) then
associate(urr => this % urr_data)
write(unit_,*) ' Unresolved resonance probability table:'
write(unit_,*) ' # of energies = ' // trim(to_str(urr % n_energy))
write(unit_,*) ' # of probabilities = ' // trim(to_str(urr % n_prob))
@ -767,17 +812,18 @@ module nuclide_header
! Calculate memory used by probability tables and add to total
size_urr = urr % n_energy * (urr % n_prob * 6 + 1) * 8
end if
end associate
end if
! Write memory used
write(unit_,*) ' Memory Requirements'
write(unit_,*) ' Cross sections = ' // trim(to_str(size_xs)) // ' bytes'
write(unit_,*) ' Probability Tables = ' // &
trim(to_str(size_urr)) // ' bytes'
! Write memory used
write(unit_,*) ' Memory Requirements'
write(unit_,*) ' Cross sections = ' // trim(to_str(size_xs)) // ' bytes'
write(unit_,*) ' Probability Tables = ' // &
trim(to_str(size_urr)) // ' bytes'
! Blank line at end of nuclide
write(unit_,*)
end subroutine nuclidece_print
! Blank line at end of nuclide
write(unit_,*)
end subroutine nuclidece_print
subroutine nuclidemg_print(this, unit_)
class(NuclideMG), intent(in) :: this

1
src/output.F90
View file

@ -2,7 +2,6 @@ module output
use, intrinsic :: ISO_FORTRAN_ENV
use ace_header, only: Reaction, UrrData
use constants
use endf, only: reaction_name
use error, only: fatal_error, warning

132
src/physics.F90
View file

@ -1,13 +1,10 @@
module physics
use ace_header, only: Reaction
use constants
use cross_section, only: elastic_xs_0K
use endf, only: reaction_name
use error, only: fatal_error, warning
use fission, only: nu_total, nu_delayed
use global
use interpolation, only: interpolate_tab1
use material_header, only: Material
use math
use mesh, only: get_mesh_indices
@ -17,6 +14,7 @@ module physics
use particle_restart_write, only: write_particle_restart
use physics_common
use random_lcg, only: prn, advance_prn_seed, prn_set_stream
use reaction_header, only: Reaction
use search, only: binary_search
use secondary_uncorrelated, only: UncorrelatedAngleEnergy
use string, only: to_str
@ -447,9 +445,9 @@ contains
vel = sqrt(dot_product(v_n, v_n))
! Sample scattering angle
select type (dist => rxn%secondary%distribution(1)%obj)
select type (dist => rxn % products(1) % distribution(1) % obj)
type is (UncorrelatedAngleEnergy)
mu_cm = dist%angle%sample(E)
mu_cm = dist % angle % sample(E)
end select
! Determine direction cosines in CM
@ -1073,8 +1071,6 @@ contains
integer :: nu ! actual number of neutrons produced
integer :: ijk(3) ! indices in ufs mesh
real(8) :: nu_t ! total nu
real(8) :: mu ! fission neutron angular cosine
real(8) :: phi ! fission neutron azimuthal angle
real(8) :: weight ! weight adjustment for ufs method
logical :: in_mesh ! source site in ufs mesh?
type(NuclideCE), pointer :: nuc
@ -1145,25 +1141,12 @@ contains
! Set weight of fission bank site
bank_array(i) % wgt = ONE/weight
! Sample cosine of angle -- fission neutrons are always emitted
! isotropically. Sometimes in ACE data, fission reactions actually have
! an angular distribution listed, but for those that do, it's simply just
! a uniform distribution in mu
mu = TWO * prn() - ONE
! Sample delayed group and angle/energy for fission reaction
call sample_fission_neutron(nuc, nuc % reactions(i_reaction), &
p % E, bank_array(i))
! Sample azimuthal angle uniformly in [0,2*pi)
phi = TWO*PI*prn()
bank_array(i) % uvw(1) = mu
bank_array(i) % uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
bank_array(i) % uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
! Sample secondary energy distribution for fission reaction and set energy
! in fission bank
bank_array(i) % E = sample_fission_energy(nuc, &
nuc % reactions(i_reaction), p)
! Set the delayed group of the neutron
bank_array(i) % delayed_group = p % delayed_group
! Set delayed group on particle too
p % delayed_group = bank_array(i) % delayed_group
! Increment the number of neutrons born delayed
if (p % delayed_group > 0) then
@ -1182,35 +1165,41 @@ contains
end subroutine create_fission_sites
!===============================================================================
! SAMPLE_FISSION_ENERGY
! SAMPLE_FISSION_NEUTRON
!===============================================================================
function sample_fission_energy(nuc, rxn, p) result(E_out)
subroutine sample_fission_neutron(nuc, rxn, E_in, site)
type(NuclideCE), intent(in) :: nuc
type(Reaction), intent(in) :: rxn
real(8), intent(in) :: E_in
type(Bank), intent(inout) :: site
type(NuclideCE), intent(in) :: nuc
type(Reaction), intent(in) :: rxn
type(Particle), intent(inout) :: p ! Particle causing fission
real(8) :: E_out ! outgoing energy of fission neutron
integer :: j ! index on nu energy grid / precursor group
integer :: lc ! index before start of energies/nu values
integer :: NR ! number of interpolation regions
integer :: NE ! number of energies tabulated
integer :: n_sample ! number of times resampling
integer :: group ! index on nu energy grid / precursor group
integer :: n_sample ! number of resamples
real(8) :: nu_t ! total nu
real(8) :: nu_d ! delayed nu
real(8) :: beta ! delayed neutron fraction
real(8) :: xi ! random number
real(8) :: yield ! delayed neutron precursor yield
real(8) :: prob ! cumulative probability
real(8) :: mu ! cosine of scattering angle
real(8) :: phi ! azimuthal angle
! Determine total nu
nu_t = nu_total(nuc, p % E)
! Sample cosine of angle -- fission neutrons are always emitted
! isotropically. Sometimes in ACE data, fission reactions actually have
! an angular distribution listed, but for those that do, it's simply just
! a uniform distribution in mu
mu = TWO * prn() - ONE
! Determine delayed nu
nu_d = nu_delayed(nuc, p % E)
! Sample azimuthal angle uniformly in [0,2*pi)
phi = TWO*PI*prn()
site % uvw(1) = mu
site % uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
site % uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
! Determine delayed neutron fraction
! Determine total nu, delayed nu, and delayed neutron fraction
nu_t = nuc % nu(E_in, EMISSION_TOTAL)
nu_d = nuc % nu(E_in, EMISSION_DELAYED)
beta = nu_d / nu_t
if (prn() < beta) then
@ -1218,51 +1207,41 @@ contains
! DELAYED NEUTRON SAMPLED
! sampled delayed precursor group
xi = prn()
lc = 1
xi = prn()*nu_d
prob = ZERO
do j = 1, nuc % n_precursor
! determine number of interpolation regions and energies
NR = int(nuc % nu_d_precursor_data(lc + 1))
NE = int(nuc % nu_d_precursor_data(lc + 2 + 2*NR))
do group = 1, nuc % n_precursor
! determine delayed neutron precursor yield for group j
yield = interpolate_tab1(nuc % nu_d_precursor_data( &
lc+1:lc+2+2*NR+2*NE), p % E)
yield = rxn % products(1 + group) % yield % evaluate(E_in)
! Check if this group is sampled
prob = prob + yield
if (xi < prob) exit
! advance pointer
lc = lc + 2 + 2*NR + 2*NE + 1
end do
! if the sum of the probabilities is slightly less than one and the
! random number is greater, j will be greater than nuc %
! n_precursor -- check for this condition
j = min(j, nuc % n_precursor)
group = min(group, nuc % n_precursor)
! set the delayed group for the particle born from fission
p % delayed_group = j
site % delayed_group = group
! sample from energy distribution
n_sample = 0
do
select type (aedist => nuc%nu_d_edist(j)%obj)
type is (UncorrelatedAngleEnergy)
E_out = aedist%energy%sample(p%E)
end select
! sample from energy/angle distribution -- note that mu has already been
! sampled above and doesn't need to be resampled
call rxn % products(1 + group) % sample(E_in, site % E, mu)
! resample if energy is greater than maximum neutron energy
if (E_out < energy_max_neutron) exit
if (site % E < energy_max_neutron) exit
! check for large number of resamples
n_sample = n_sample + 1
if (n_sample == MAX_SAMPLE) then
! call write_particle_restart(p)
call fatal_error("Resampled energy distribution maximum number of " &
&// "times for nuclide " // nuc % name)
// "times for nuclide " // nuc % name)
end if
end do
@ -1271,28 +1250,27 @@ contains
! PROMPT NEUTRON SAMPLED
! set the delayed group for the particle born from fission to 0
p % delayed_group = 0
site % delayed_group = 0
! sample from prompt neutron energy distribution
n_sample = 0
do
call rxn%secondary%sample(p%E, E_out, prob)
call rxn % products(1) % sample(E_in, site % E, mu)
! resample if energy is greater than maximum neutron energy
if (E_out < energy_max_neutron) exit
if (site % E < energy_max_neutron) exit
! check for large number of resamples
n_sample = n_sample + 1
if (n_sample == MAX_SAMPLE) then
! call write_particle_restart(p)
call fatal_error("Resampled energy distribution maximum number of " &
&// "times for nuclide " // nuc % name)
// "times for nuclide " // nuc % name)
end if
end do
end if
end function sample_fission_energy
end subroutine sample_fission_neutron
!===============================================================================
! INELASTIC_SCATTER handles all reactions with a single secondary neutron (other
@ -1316,7 +1294,7 @@ contains
E_in = p % E
! sample outgoing energy and scattering cosine
call rxn%secondary%sample(E_in, E, mu)
call rxn % products(1) % sample(E_in, E, mu)
! if scattering system is in center-of-mass, transfer cosine of scattering
! angle and outgoing energy from CM to LAB
@ -1344,14 +1322,16 @@ contains
! change direction of particle
p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu)
! change weight of particle based on yield
if (rxn % multiplicity_with_E) then
yield = interpolate_tab1(rxn % multiplicity_E, E_in)
p % wgt = yield * p % wgt
else
do i = 1, rxn % multiplicity - 1
call p % create_secondary(p % coord(1) % uvw, NEUTRON, run_CE=.True.)
! evaluate yield
yield = rxn % products(1) % yield % evaluate(E_in)
if (mod(yield, ONE) == ZERO) then
! If yield is integral, create exactly that many secondary particles
do i = 1, nint(yield) - 1
call p % create_secondary(p % coord(1) % uvw, NEUTRON, run_CE=.true.)
end do
else
! Otherwise, change weight of particle based on yield
p % wgt = yield * p % wgt
end if
end subroutine inelastic_scatter

62
src/product_header.F90 Normal file
View file

@ -0,0 +1,62 @@
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, Constant1D, Polynomial
use random_lcg, only: prn
!===============================================================================
! 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
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
end module product_header

21
src/reaction_header.F90 Normal file
View file

@ -0,0 +1,21 @@
module reaction_header
use product_header, only: ReactionProduct
implicit none
!===============================================================================
! REACTION contains the cross-section and secondary energy and angle
! distributions for a single reaction in a continuous-energy ACE-format table
!===============================================================================
type Reaction
integer :: MT ! ENDF MT value
real(8) :: Q_value ! Reaction Q value
integer :: threshold ! Energy grid index of threshold
logical :: scatter_in_cm ! scattering system in center-of-mass?
real(8), allocatable :: sigma(:) ! Cross section values
type(ReactionProduct), allocatable :: products(:)
end type Reaction
end module reaction_header

52
src/secondary_correlated.F90
View file

@ -1,8 +1,8 @@
module secondary_correlated
use angleenergy_header, only: AngleEnergy
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
@ -24,8 +24,8 @@ module secondary_correlated
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
real(8), allocatable :: energy(:) ! incoming energies
type(AngleEnergyTable), allocatable :: distribution(:) ! outgoing E/mu distributions
contains
procedure :: sample => correlated_sample
end type CorrelatedAngleEnergy
@ -61,17 +61,17 @@ contains
! 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
n_energy_in = size(this%energy)
if (E_in < this%energy(1)) then
i = 1
r = ZERO
elseif (E_in > this%energy_in(n_energy_in)) then
elseif (E_in > this%energy(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))
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
@ -82,23 +82,23 @@ contains
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%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%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)
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%table(l)%e_out)
n_energy_out = size(this%distribution(l)%e_out)
r1 = prn()
c_k = this%table(l)%c(1)
c_k = this%distribution(l)%c(1)
do k = 1, n_energy_out - 1
c_k1 = this%table(l)%c(k+1)
c_k1 = this%distribution(l)%c(k+1)
if (r1 < c_k1) exit
c_k = c_k1
end do
@ -106,9 +106,9 @@ contains
! 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
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
@ -116,10 +116,10 @@ contains
E_out = E_l_k
end if
elseif (this%table(l)%interpolation == LINEAR_LINEAR) then
elseif (this%distribution(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)
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
@ -139,9 +139,9 @@ contains
! Find correlated angular distribution for closest outgoing energy bin
if (r1 - c_k < c_k1 - r1) then
mu = this%table(l)%angle(k)%obj%sample()
mu = this%distribution(l)%angle(k)%obj%sample()
else
mu = this%table(l)%angle(k + 1)%obj%sample()
mu = this%distribution(l)%angle(k + 1)%obj%sample()
end if
end subroutine correlated_sample

83
src/secondary_header.F90
View file

@ -1,83 +0,0 @@
module secondary_header
use constants, only: ZERO
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(angleenergy_sample_), deferred :: sample
end type AngleEnergy
abstract interface
subroutine angleenergy_sample_(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 angleenergy_sample_
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 :: 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 + interpolate_tab1(this%applicability(i), 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 secondary_sample
end module secondary_header

60
src/secondary_kalbach.F90
View file

@ -1,7 +1,7 @@
module secondary_kalbach
use angleenergy_header, only: AngleEnergy
use constants, only: ZERO, ONE, TWO, HISTOGRAM, LINEAR_LINEAR
use secondary_header, only: AngleEnergy
use random_lcg, only: prn
use search, only: binary_search
@ -25,8 +25,8 @@ module secondary_kalbach
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
real(8), allocatable :: energy(:) ! incoming energies
type(KalbachMannTable), allocatable :: distribution(:) ! outgoing E/mu parameters
contains
procedure :: sample => kalbachmann_sample
end type KalbachMann
@ -64,17 +64,17 @@ contains
! 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
n_energy_in = size(this%energy)
if (E_in < this%energy(1)) then
i = 1
r = ZERO
elseif (E_in > this%energy_in(n_energy_in)) then
elseif (E_in > this%energy(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))
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
@ -85,23 +85,23 @@ contains
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%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%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)
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%table(l)%e_out)
n_energy_out = size(this%distribution(l)%e_out)
r1 = prn()
c_k = this%table(l)%c(1)
c_k = this%distribution(l)%c(1)
do k = 1, n_energy_out - 1
c_k1 = this%table(l)%c(k+1)
c_k1 = this%distribution(l)%c(k+1)
if (r1 < c_k1) exit
c_k = c_k1
end do
@ -109,9 +109,9 @@ contains
! 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
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
@ -120,13 +120,13 @@ contains
end if
! Determine Kalbach-Mann parameters
km_r = this%table(l)%r(k)
km_a = this%table(l)%a(k)
km_r = this%distribution(l)%r(k)
km_a = this%distribution(l)%a(k)
elseif (this%table(l)%interpolation == LINEAR_LINEAR) then
elseif (this%distribution(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)
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
@ -137,10 +137,10 @@ contains
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))
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

70
src/secondary_nbody.F90 Normal file
View file

@ -0,0 +1,70 @@
module secondary_nbody
use angleenergy_header, only: AngleEnergy
use constants, only: ONE, TWO, PI
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
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
end module secondary_nbody

2
src/secondary_uncorrelated.F90
View file

@ -1,9 +1,9 @@
module secondary_uncorrelated
use angle_distribution, only: AngleDistribution
use angleenergy_header, only: AngleEnergy
use constants, only: ONE, TWO
use energy_distribution, only: EnergyDistribution
use secondary_header, only: AngleEnergy
use random_lcg, only: prn
!===============================================================================

1
src/summary.F90
View file

@ -1,6 +1,5 @@
module summary
use ace_header, only: Reaction, UrrData
use constants
use endf, only: reaction_name
use geometry_header, only: Cell, Universe, Lattice, RectLattice, &

95
src/tally.F90
View file

@ -1,7 +1,7 @@
module tally
use ace_header, only: Reaction
use constants
use endf_header, only: Constant1D
use error, only: fatal_error
use geometry_header
use global
@ -15,8 +15,6 @@ module tally
use search, only: binary_search
use string, only: to_str
use tally_header, only: TallyResult, TallyMapItem, TallyMapElement
use fission, only: nu_total, nu_delayed, yield_delayed
use interpolation, only: interpolate_tab1
#ifdef MPI
use message_passing
@ -245,9 +243,9 @@ contains
! Only analog estimators are available.
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP
! For scattering production, we need to use the pre-collision
! weight times the multiplicity as the estimate for the number of
! neutrons exiting a reaction with neutrons in the exit channel
! For scattering production, we need to use the pre-collision weight
! times the yield as the estimate for the number of neutrons exiting a
! reaction with neutrons in the exit channel
if (p % event_MT == ELASTIC .or. p % event_MT == N_LEVEL .or. &
(p % event_MT >= N_N1 .and. p % event_MT <= N_NC)) then
! Don't waste time on very common reactions we know have multiplicities
@ -257,16 +255,17 @@ contains
m = nuclides(p%event_nuclide)%reaction_index% &
get_key(p % event_MT)
! Get multiplicity and apply to score
! Get yield and apply to score
associate (rxn => nuclides(p%event_nuclide)%reactions(m))
if (rxn % multiplicity_with_E) then
! Then the multiplicity was already incorporated in to p % wgt
! per the scattering routine,
select type (yield => rxn % products(1) % yield)
type is (Constant1D)
! Grab the yield from the reaction
score = p % last_wgt * yield % y
class default
! the yield was already incorporated in to p % wgt per the
! scattering routine
score = p % wgt
else
! Grab the multiplicity from the rxn
score = p % last_wgt * rxn % multiplicity
end if
end select
end associate
end if
@ -279,7 +278,7 @@ contains
cycle SCORE_LOOP
end if
! For scattering production, we need to use the pre-collision
! weight times the multiplicity as the estimate for the number of
! weight times the yield as the estimate for the number of
! neutrons exiting a reaction with neutrons in the exit channel
if (p % event_MT == ELASTIC .or. p % event_MT == N_LEVEL .or. &
(p % event_MT >= N_N1 .and. p % event_MT <= N_NC)) then
@ -290,16 +289,17 @@ contains
m = nuclides(p%event_nuclide)%reaction_index% &
get_key(p % event_MT)
! Get multiplicity and apply to score
! Get yield and apply to score
associate (rxn => nuclides(p%event_nuclide)%reactions(m))
if (rxn % multiplicity_with_E) then
! Then the multiplicity was already incorporated in to p % wgt
! per the scattering routine,
select type (yield => rxn % products(1) % yield)
type is (Constant1D)
! Grab the yield from the reaction
score = p % last_wgt * yield % y
class default
! the yield was already incorporated in to p % wgt per the
! scattering routine
score = p % wgt
else
! Grab the multiplicity from the rxn
score = p % last_wgt * rxn % multiplicity
end if
end select
end associate
end if
@ -312,7 +312,7 @@ contains
cycle SCORE_LOOP
end if
! For scattering production, we need to use the pre-collision
! weight times the multiplicity as the estimate for the number of
! weight times the yield as the estimate for the number of
! neutrons exiting a reaction with neutrons in the exit channel
if (p % event_MT == ELASTIC .or. p % event_MT == N_LEVEL .or. &
(p % event_MT >= N_N1 .and. p % event_MT <= N_NC)) then
@ -323,16 +323,17 @@ contains
m = nuclides(p%event_nuclide)%reaction_index% &
get_key(p % event_MT)
! Get multiplicity and apply to score
! Get yield and apply to score
associate (rxn => nuclides(p%event_nuclide)%reactions(m))
if (rxn % multiplicity_with_E) then
! Then the multiplicity was already incorporated in to p % wgt
! per the scattering routine,
select type (yield => rxn % products(1) % yield)
type is (Constant1D)
! Grab the yield from the reaction
score = p % last_wgt * yield % y
class default
! the yield was already incorporated in to p % wgt per the
! scattering routine
score = p % wgt
else
! Grab the multiplicity from the rxn
score = p % last_wgt * rxn % multiplicity
end if
end select
end associate
end if
@ -499,12 +500,11 @@ contains
d = t % filters(dg_filter) % int_bins(d_bin)
! Compute the yield for this delayed group
yield = yield_delayed(nuclides(p % event_nuclide), E, d)
yield = nuclides(p % event_nuclide) % nu(E, EMISSION_DELAYED, d)
! Compute the score and tally to bin
score = p % absorb_wgt * yield * micro_xs(p % event_nuclide) &
% fission * nu_delayed(nuclides(p % event_nuclide), E) / &
micro_xs(p % event_nuclide) % absorption
% fission / micro_xs(p % event_nuclide) % absorption
call score_fission_delayed_dg(t, d_bin, score, score_index)
end do
cycle SCORE_LOOP
@ -512,9 +512,9 @@ contains
! If the delayed group filter is not present, compute the score
! by multiplying the absorbed weight by the fraction of the
! delayed-nu-fission xs to the absorption xs
score = p % absorb_wgt * micro_xs(p % event_nuclide) &
% fission * nu_delayed(nuclides(p % event_nuclide), E) / &
micro_xs(p % event_nuclide) % absorption
score = p % absorb_wgt * micro_xs(p % event_nuclide) % fission &
* nuclides(p % event_nuclide) % nu(E, EMISSION_DELAYED) &
/ micro_xs(p % event_nuclide) % absorption
end if
end if
else
@ -564,11 +564,11 @@ contains
d = t % filters(dg_filter) % int_bins(d_bin)
! Compute the yield for this delayed group
yield = yield_delayed(nuclides(i_nuclide), E, d)
yield = nuclides(i_nuclide) % nu(E, EMISSION_DELAYED, d)
! Compute the score and tally to bin
score = micro_xs(i_nuclide) % fission * yield &
* nu_delayed(nuclides(i_nuclide), E) * atom_density * flux
score = micro_xs(i_nuclide) % fission * yield * &
atom_density * flux
call score_fission_delayed_dg(t, d_bin, score, score_index)
end do
cycle SCORE_LOOP
@ -576,8 +576,8 @@ contains
! If the delayed group filter is not present, compute the score
! by multiplying the delayed-nu-fission macro xs by the flux
score = micro_xs(i_nuclide) % fission * &
nu_delayed(nuclides(i_nuclide), E) * atom_density * flux
score = micro_xs(i_nuclide) % fission * nuclides(i_nuclide) % &
nu(E, EMISSION_DELAYED) * atom_density * flux
end if
! Tally is on total nuclides
@ -602,11 +602,10 @@ contains
d = t % filters(dg_filter) % int_bins(d_bin)
! Get the yield for the desired nuclide and delayed group
yield = yield_delayed(nuclides(i_nuc), E, d)
yield = nuclides(i_nuc) % nu(E, EMISSION_DELAYED, d)
! Compute the score and tally to bin
score = micro_xs(i_nuc) % fission * yield &
* nu_delayed(nuclides(i_nuc), E) * atom_density_ * flux
score = micro_xs(i_nuc) % fission * yield * atom_density_ * flux
call score_fission_delayed_dg(t, d_bin, score, score_index)
end do
end do
@ -625,8 +624,8 @@ contains
i_nuc = materials(p % material) % nuclide(l)
! Accumulate the contribution from each nuclide
score = score + micro_xs(i_nuc) % fission &
* nu_delayed(nuclides(i_nuc), E) * atom_density_ * flux
score = score + micro_xs(i_nuc) % fission * nuclides(i_nuc) % &
nu(E, EMISSION_DELAYED) * atom_density_ * flux
end do
end if
end if

20
src/urr_header.F90 Normal file
View file

@ -0,0 +1,20 @@
module urr_header
implicit none
!===============================================================================
! URRDATA contains probability tables for the unresolved resonance range.
!===============================================================================
type UrrData
integer :: n_energy ! # of incident neutron energies
integer :: n_prob ! # of probabilities
integer :: interp ! inteprolation (2=lin-lin, 5=log-log)
integer :: inelastic_flag ! inelastic competition flag
integer :: absorption_flag ! other absorption flag
logical :: multiply_smooth ! multiply by smooth cross section?
real(8), allocatable :: energy(:) ! incident energies
real(8), allocatable :: prob(:,:,:) ! actual probabibility tables
end type UrrData
end module urr_header

62
tests/test_mgxs_library_condense/results_true.dat
View file

@ -2,48 +2,48 @@
0 1 1 total 0.412084 0.02359 material group in nuclide mean std. dev.
0 1 1 total 0.076425 0.003691 material group in group out nuclide mean std. dev.
0 1 1 1 total 0.345643 0.021487 material group out nuclide mean std. dev.
0 1 1 total 1 0.055333 material group in nuclide mean std. dev.
0 1 1 total 1.0 0.055333 material group in nuclide mean std. dev.
0 2 1 total 0.241262 0.00841 material group in nuclide mean std. dev.
0 2 1 total 0 0 material group in group out nuclide mean std. dev.
0 2 1 total 0.0 0.0 material group in group out nuclide mean std. dev.
0 2 1 1 total 0.241262 0.00841 material group out nuclide mean std. dev.
0 2 1 total 0 0 material group in nuclide mean std. dev.
0 2 1 total 0.0 0.0 material group in nuclide mean std. dev.
0 3 1 total 0.400028 0.034667 material group in nuclide mean std. dev.
0 3 1 total 0 0 material group in group out nuclide mean std. dev.
0 3 1 total 0.0 0.0 material group in group out nuclide mean std. dev.
0 3 1 1 total 0.393462 0.033646 material group out nuclide mean std. dev.
0 3 1 total 0 0 material group in nuclide mean std. dev.
0 3 1 total 0.0 0.0 material group in nuclide mean std. dev.
0 4 1 total 0.377402 0.072937 material group in nuclide mean std. dev.
0 4 1 total 0 0 material group in group out nuclide mean std. dev.
0 4 1 total 0.0 0.0 material group in group out nuclide mean std. dev.
0 4 1 1 total 0.371473 0.071226 material group out nuclide mean std. dev.
0 4 1 total 0 0 material group in nuclide mean std. dev.
0 5 1 total 0 0 material group in nuclide mean std. dev.
0 5 1 total 0 0 material group in group out nuclide mean std. dev.
0 5 1 1 total 0 0 material group out nuclide mean std. dev.
0 5 1 total 0 0 material group in nuclide mean std. dev.
0 6 1 total 0 0 material group in nuclide mean std. dev.
0 6 1 total 0 0 material group in group out nuclide mean std. dev.
0 6 1 1 total 0 0 material group out nuclide mean std. dev.
0 6 1 total 0 0 material group in nuclide mean std. dev.
0 7 1 total 0 0 material group in nuclide mean std. dev.
0 7 1 total 0 0 material group in group out nuclide mean std. dev.
0 7 1 1 total 0 0 material group out nuclide mean std. dev.
0 7 1 total 0 0 material group in nuclide mean std. dev.
0 8 1 total 0 0 material group in nuclide mean std. dev.
0 8 1 total 0 0 material group in group out nuclide mean std. dev.
0 8 1 1 total 0 0 material group out nuclide mean std. dev.
0 8 1 total 0 0 material group in nuclide mean std. dev.
0 4 1 total 0.0 0.0 material group in nuclide mean std. dev.
0 5 1 total 0.0 0.0 material group in nuclide mean std. dev.
0 5 1 total 0.0 0.0 material group in group out nuclide mean std. dev.
0 5 1 1 total 0.0 0.0 material group out nuclide mean std. dev.
0 5 1 total 0.0 0.0 material group in nuclide mean std. dev.
0 6 1 total 0.0 0.0 material group in nuclide mean std. dev.
0 6 1 total 0.0 0.0 material group in group out nuclide mean std. dev.
0 6 1 1 total 0.0 0.0 material group out nuclide mean std. dev.
0 6 1 total 0.0 0.0 material group in nuclide mean std. dev.
0 7 1 total 0.0 0.0 material group in nuclide mean std. dev.
0 7 1 total 0.0 0.0 material group in group out nuclide mean std. dev.
0 7 1 1 total 0.0 0.0 material group out nuclide mean std. dev.
0 7 1 total 0.0 0.0 material group in nuclide mean std. dev.
0 8 1 total 0.0 0.0 material group in nuclide mean std. dev.
0 8 1 total 0.0 0.0 material group in group out nuclide mean std. dev.
0 8 1 1 total 0.0 0.0 material group out nuclide mean std. dev.
0 8 1 total 0.0 0.0 material group in nuclide mean std. dev.
0 9 1 total 0.600536 0.748875 material group in nuclide mean std. dev.
0 9 1 total 0 0 material group in group out nuclide mean std. dev.
0 9 1 total 0.0 0.0 material group in group out nuclide mean std. dev.
0 9 1 1 total 0.600536 0.748875 material group out nuclide mean std. dev.
0 9 1 total 0 0 material group in nuclide mean std. dev.
0 9 1 total 0.0 0.0 material group in nuclide mean std. dev.
0 10 1 total 0.235515 0.613974 material group in nuclide mean std. dev.
0 10 1 total 0 0 material group in group out nuclide mean std. dev.
0 10 1 total 0.0 0.0 material group in group out nuclide mean std. dev.
0 10 1 1 total 0.235515 0.613974 material group out nuclide mean std. dev.
0 10 1 total 0 0 material group in nuclide mean std. dev.
0 10 1 total 0.0 0.0 material group in nuclide mean std. dev.
0 11 1 total 0.510145 0.741941 material group in nuclide mean std. dev.
0 11 1 total 0 0 material group in group out nuclide mean std. dev.
0 11 1 total 0.0 0.0 material group in group out nuclide mean std. dev.
0 11 1 1 total 0.491857 0.715554 material group out nuclide mean std. dev.
0 11 1 total 0 0 material group in nuclide mean std. dev.
0 11 1 total 0.0 0.0 material group in nuclide mean std. dev.
0 12 1 total 0.73836 0.825631 material group in nuclide mean std. dev.
0 12 1 total 0 0 material group in group out nuclide mean std. dev.
0 12 1 total 0.0 0.0 material group in group out nuclide mean std. dev.
0 12 1 1 total 0.723265 0.808231 material group out nuclide mean std. dev.
0 12 1 total 0 0
0 12 1 total 0.0 0.0

4
tests/test_mgxs_library_distribcell/results_true.dat
View file

@ -1,5 +1,5 @@
avg(distribcell) group in nuclide mean std. dev.
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.718919 0.520644 avg(distribcell) group in nuclide mean std. dev.
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0 0 avg(distribcell) group in group out nuclide mean std. dev.
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.0 0.0 avg(distribcell) group in group out nuclide mean std. dev.
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total 0.695166 0.510606 avg(distribcell) group out nuclide mean std. dev.
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0 0
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.0 0.0

140
tests/test_mgxs_library_no_nuclides/results_true.dat
View file

@ -7,115 +7,115 @@
2 1 1 2 total 0.001559 0.000510
1 1 2 1 total 0.000000 0.000000
0 1 2 2 total 0.422051 0.021617 material group out nuclide mean std. dev.
1 1 1 total 1 0.055333
0 1 2 total 0 0.000000 material group in nuclide mean std. dev.
1 1 1 total 1.0 0.055333
0 1 2 total 0.0 0.000000 material group in nuclide mean std. dev.
1 2 1 total 0.237254 0.008184
0 2 2 total 0.285930 0.048796 material group in nuclide mean std. dev.
1 2 1 total 0 0
0 2 2 total 0 0 material group in group out nuclide mean std. dev.
1 2 1 total 0.0 0.0
0 2 2 total 0.0 0.0 material group in group out nuclide mean std. dev.
3 2 1 1 total 0.237254 0.008184
2 2 1 2 total 0.000000 0.000000
1 2 2 1 total 0.000000 0.000000
0 2 2 2 total 0.285930 0.048796 material group out nuclide mean std. dev.
1 2 1 total 0 0
0 2 2 total 0 0 material group in nuclide mean std. dev.
1 2 1 total 0.0 0.0
0 2 2 total 0.0 0.0 material group in nuclide mean std. dev.
1 3 1 total 0.286906 0.027401
0 3 2 total 1.418151 0.265308 material group in nuclide mean std. dev.
1 3 1 total 0 0
0 3 2 total 0 0 material group in group out nuclide mean std. dev.
1 3 1 total 0.0 0.0
0 3 2 total 0.0 0.0 material group in group out nuclide mean std. dev.
3 3 1 1 total 0.259937 0.026115
2 3 1 2 total 0.026187 0.001665
1 3 2 1 total 0.000000 0.000000
0 3 2 2 total 1.359521 0.258505 material group out nuclide mean std. dev.
1 3 1 total 0 0
0 3 2 total 0 0 material group in nuclide mean std. dev.
1 3 1 total 0.0 0.0
0 3 2 total 0.0 0.0 material group in nuclide mean std. dev.
1 4 1 total 0.242447 0.061031
0 4 2 total 1.253959 0.388363 material group in nuclide mean std. dev.
1 4 1 total 0 0
0 4 2 total 0 0 material group in group out nuclide mean std. dev.
1 4 1 total 0.0 0.0
0 4 2 total 0.0 0.0 material group in group out nuclide mean std. dev.
3 4 1 1 total 0.217930 0.058565
2 4 1 2 total 0.023662 0.003083
1 4 2 1 total 0.000000 0.000000
0 4 2 2 total 1.215074 0.381025 material group out nuclide mean std. dev.
1 4 1 total 0 0
0 4 2 total 0 0 material group in nuclide mean std. dev.
1 5 1 total 0 0
0 5 2 total 0 0 material group in nuclide mean std. dev.
1 5 1 total 0 0
0 5 2 total 0 0 material group in group out nuclide mean std. dev.
3 5 1 1 total 0 0
2 5 1 2 total 0 0
1 5 2 1 total 0 0
0 5 2 2 total 0 0 material group out nuclide mean std. dev.
1 5 1 total 0 0
0 5 2 total 0 0 material group in nuclide mean std. dev.
1 6 1 total 0 0
0 6 2 total 0 0 material group in nuclide mean std. dev.
1 6 1 total 0 0
0 6 2 total 0 0 material group in group out nuclide mean std. dev.
3 6 1 1 total 0 0
2 6 1 2 total 0 0
1 6 2 1 total 0 0
0 6 2 2 total 0 0 material group out nuclide mean std. dev.
1 6 1 total 0 0
0 6 2 total 0 0 material group in nuclide mean std. dev.
1 7 1 total 0 0
0 7 2 total 0 0 material group in nuclide mean std. dev.
1 7 1 total 0 0
0 7 2 total 0 0 material group in group out nuclide mean std. dev.
3 7 1 1 total 0 0
2 7 1 2 total 0 0
1 7 2 1 total 0 0
0 7 2 2 total 0 0 material group out nuclide mean std. dev.
1 7 1 total 0 0
0 7 2 total 0 0 material group in nuclide mean std. dev.
1 8 1 total 0 0
0 8 2 total 0 0 material group in nuclide mean std. dev.
1 8 1 total 0 0
0 8 2 total 0 0 material group in group out nuclide mean std. dev.
3 8 1 1 total 0 0
2 8 1 2 total 0 0
1 8 2 1 total 0 0
0 8 2 2 total 0 0 material group out nuclide mean std. dev.
1 8 1 total 0 0
0 8 2 total 0 0 material group in nuclide mean std. dev.
1 4 1 total 0.0 0.0
0 4 2 total 0.0 0.0 material group in nuclide mean std. dev.
1 5 1 total 0.0 0.0
0 5 2 total 0.0 0.0 material group in nuclide mean std. dev.
1 5 1 total 0.0 0.0
0 5 2 total 0.0 0.0 material group in group out nuclide mean std. dev.
3 5 1 1 total 0.0 0.0
2 5 1 2 total 0.0 0.0
1 5 2 1 total 0.0 0.0
0 5 2 2 total 0.0 0.0 material group out nuclide mean std. dev.
1 5 1 total 0.0 0.0
0 5 2 total 0.0 0.0 material group in nuclide mean std. dev.
1 6 1 total 0.0 0.0
0 6 2 total 0.0 0.0 material group in nuclide mean std. dev.
1 6 1 total 0.0 0.0
0 6 2 total 0.0 0.0 material group in group out nuclide mean std. dev.
3 6 1 1 total 0.0 0.0
2 6 1 2 total 0.0 0.0
1 6 2 1 total 0.0 0.0
0 6 2 2 total 0.0 0.0 material group out nuclide mean std. dev.
1 6 1 total 0.0 0.0
0 6 2 total 0.0 0.0 material group in nuclide mean std. dev.
1 7 1 total 0.0 0.0
0 7 2 total 0.0 0.0 material group in nuclide mean std. dev.
1 7 1 total 0.0 0.0
0 7 2 total 0.0 0.0 material group in group out nuclide mean std. dev.
3 7 1 1 total 0.0 0.0
2 7 1 2 total 0.0 0.0
1 7 2 1 total 0.0 0.0
0 7 2 2 total 0.0 0.0 material group out nuclide mean std. dev.
1 7 1 total 0.0 0.0
0 7 2 total 0.0 0.0 material group in nuclide mean std. dev.
1 8 1 total 0.0 0.0
0 8 2 total 0.0 0.0 material group in nuclide mean std. dev.
1 8 1 total 0.0 0.0
0 8 2 total 0.0 0.0 material group in group out nuclide mean std. dev.
3 8 1 1 total 0.0 0.0
2 8 1 2 total 0.0 0.0
1 8 2 1 total 0.0 0.0
0 8 2 2 total 0.0 0.0 material group out nuclide mean std. dev.
1 8 1 total 0.0 0.0
0 8 2 total 0.0 0.0 material group in nuclide mean std. dev.
1 9 1 total 0.600536 0.748875
0 9 2 total 0.000000 0.000000 material group in nuclide mean std. dev.
1 9 1 total 0 0
0 9 2 total 0 0 material group in group out nuclide mean std. dev.
1 9 1 total 0.0 0.0
0 9 2 total 0.0 0.0 material group in group out nuclide mean std. dev.
3 9 1 1 total 0.600536 0.748875
2 9 1 2 total 0.000000 0.000000
1 9 2 1 total 0.000000 0.000000
0 9 2 2 total 0.000000 0.000000 material group out nuclide mean std. dev.
1 9 1 total 0 0
0 9 2 total 0 0 material group in nuclide mean std. dev.
1 9 1 total 0.0 0.0
0 9 2 total 0.0 0.0 material group in nuclide mean std. dev.
1 10 1 total 0.235515 0.613974
0 10 2 total 0.000000 0.000000 material group in nuclide mean std. dev.
1 10 1 total 0 0
0 10 2 total 0 0 material group in group out nuclide mean std. dev.
1 10 1 total 0.0 0.0
0 10 2 total 0.0 0.0 material group in group out nuclide mean std. dev.
3 10 1 1 total 0.235515 0.613974
2 10 1 2 total 0.000000 0.000000
1 10 2 1 total 0.000000 0.000000
0 10 2 2 total 0.000000 0.000000 material group out nuclide mean std. dev.
1 10 1 total 0 0
0 10 2 total 0 0 material group in nuclide mean std. dev.
1 10 1 total 0.0 0.0
0 10 2 total 0.0 0.0 material group in nuclide mean std. dev.
1 11 1 total 0.186324 0.632129
0 11 2 total 0.945986 1.591133 material group in nuclide mean std. dev.
1 11 1 total 0 0
0 11 2 total 0 0 material group in group out nuclide mean std. dev.
1 11 1 total 0.0 0.0
0 11 2 total 0.0 0.0 material group in group out nuclide mean std. dev.
3 11 1 1 total 0.154449 0.597686
2 11 1 2 total 0.031875 0.045078
1 11 2 1 total 0.000000 0.000000
0 11 2 2 total 0.903085 1.532144 material group out nuclide mean std. dev.
1 11 1 total 0 0
0 11 2 total 0 0 material group in nuclide mean std. dev.
1 11 1 total 0.0 0.0
0 11 2 total 0.0 0.0 material group in nuclide mean std. dev.
1 12 1 total 0.213292 0.271444
0 12 2 total 1.390975 2.137346 material group in nuclide mean std. dev.
1 12 1 total 0 0
0 12 2 total 0 0 material group in group out nuclide mean std. dev.
1 12 1 total 0.0 0.0
0 12 2 total 0.0 0.0 material group in group out nuclide mean std. dev.
3 12 1 1 total 0.186052 0.257633
2 12 1 2 total 0.027240 0.029555
1 12 2 1 total 0.000000 0.000000
0 12 2 2 total 1.357118 2.089846 material group out nuclide mean std. dev.
1 12 1 total 0 0
0 12 2 total 0 0
1 12 1 total 0.0 0.0
0 12 2 total 0.0 0.0

2520
tests/test_mgxs_library_nuclides/results_true.dat
View file

File diff suppressed because it is too large Load diff

2
tests/test_tallies/results_true.dat
View file

@ -1 +1 @@
80bb207ab79131ff264a205703fcc798e3353dbead81e39dadf262979d6d6ad786123588e330c8d0bccddbcb7b7ce9af8447c73a317174019977d2392edf31f6
f1b2b43197e1bbb305000d5a84c228361afb876d23ed866cdb073fe7410335c87fb16066c031d0e4397225321632566c00f48eac6187d59bdeab9a8c60986c3c