mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 22:26:08 -04:00
fixed conflict, merged develop
This commit is contained in:
walshjon
commit
2f1952922d
10 changed files with 318 additions and 111 deletions
|
|
@ -790,6 +790,7 @@ outgoing angle is
|
|||
|
||||
\mu = \frac{1}{A} \ln \left ( \xi_4 e^A + (1 - \xi_4) e^{-A} \right ).
|
||||
|
||||
.. _ace-law-61:
|
||||
|
||||
ACE Law 61 - Correlated Energy and Angle Distribution
|
||||
+++++++++++++++++++++++++++++++++++++++++++++++++++++
|
||||
|
|
@ -952,7 +953,7 @@ as
|
|||
|
||||
v_n \bar{\sigma} (v_n, T) = \int d\mathbf{v}_T v_r \sigma(v_r)
|
||||
M (\mathbf{v}_T)
|
||||
|
||||
|
||||
where :math:`v_n` is the magnitude of the velocity of the neutron,
|
||||
:math:`\bar{\sigma}` is an effective cross section, :math:`T` is the temperature
|
||||
of the target material, :math:`\mathbf{v}_T` is the velocity of the target
|
||||
|
|
@ -1321,7 +1322,7 @@ given analytically by
|
|||
|
||||
\mu = 1 - \frac{E_i}{E}
|
||||
|
||||
where :math:`E_i` is the energy of the Bragg edge that scattered the neutron.
|
||||
where :math:`E_i` is the energy of the Bragg edge that scattered the neutron.
|
||||
|
||||
Outgoing Angle for Incoherent Elastic Scattering
|
||||
------------------------------------------------
|
||||
|
|
@ -1348,18 +1349,24 @@ where the interpolation factor is defined as
|
|||
Outgoing Energy and Angle for Inelastic Scattering
|
||||
--------------------------------------------------
|
||||
|
||||
On each |sab| table, there is a correlated angle-energy secondary distribution
|
||||
for neutron thermal inelastic scattering. While the documentation for the ACE
|
||||
format implies that there are a series of equiprobable outgoing energies, the
|
||||
outgoing energies may have non-uniform probability distribution. In particular,
|
||||
if the thermal data were processed with :math:`iwt = 0` in NJOY, then the first
|
||||
and last outgoing energies have a relative probability of 1, the second and
|
||||
second to last energies have a relative probability of 4, and all other energies
|
||||
have a relative probability of 10. The procedure to determine the outgoing
|
||||
energy and angle is as such. First, the interpolation factor is determined from
|
||||
equation :eq:`sab-interpolation-factor`. Then, an outgoing energy bin is sampled
|
||||
either from a uniform distribution or from the aforementioned skewed
|
||||
distribution. The outgoing energy is then interpolated between values
|
||||
Each |sab| table provides a correlated angle-energy secondary distribution for
|
||||
neutron thermal inelastic scattering. There are three representations used
|
||||
in the ACE thermal scattering data: equiprobable discrete outgoing
|
||||
energies, non-uniform yet still discrete outgoing energies, and continuous
|
||||
outgoing energies with corresponding probability and cumulative distribution
|
||||
functions provided in tabular format. These three representations all
|
||||
represent the angular distribution in a common format, using a series of
|
||||
discrete equiprobable outgoing cosines.
|
||||
|
||||
Equi-Probable Outgoing Energies
|
||||
+++++++++++++++++++++++++++++++
|
||||
|
||||
If the thermal data was processed with :math:`iwt = 1` in NJOY, then the
|
||||
outgoing energy spectra is represented in the ACE data as a set of discrete and
|
||||
equiprobable outgoing energies. The procedure to determine the outgoing energy
|
||||
and angle is as such. First, the interpolation factor is determined from
|
||||
equation :eq:`sab-interpolation-factor`. Then, an outgoing energy bin is
|
||||
sampled from a uniform distribution and then interpolated between values
|
||||
corresponding to neighboring incoming energies:
|
||||
|
||||
.. math::
|
||||
|
|
@ -1380,6 +1387,37 @@ uniformly and then the final cosine is interpolated on the incoming energy grid:
|
|||
where :math:`\mu_{i,j,k}` is the k-th outgoing cosine corresponding to the j-th
|
||||
outgoing energy and the i-th incoming energy.
|
||||
|
||||
Skewed Equi-Probable Outgoing Energies
|
||||
++++++++++++++++++++++++++++++++++++++
|
||||
|
||||
If the thermal data was processed with :math:`iwt=0` in NJOY, then the
|
||||
outgoing energy spectra is represented in the ACE data according to the
|
||||
following: the first and last outgoing energies have a relative probability of
|
||||
1, the second and second-to-last energies have a relative probability of 4, and
|
||||
all other energies have a relative probability of 10. The procedure to
|
||||
determine the outgoing energy and angle is similar to the method discussed
|
||||
above, except that the sampled probability distribution is now skewed
|
||||
accordingly.
|
||||
|
||||
Continuous Outgoing Energies
|
||||
++++++++++++++++++++++++++++
|
||||
|
||||
If the thermal data was processed with :math:`iwt=2` in NJOY, then the
|
||||
outgoing energy spectra is represented by a continuous outgoing energy spectra
|
||||
in tabular form with linear-linear interpolation. The sampling of the outgoing
|
||||
energy portion of this format is very similar to :ref:`ACE Law 61<ace-law-61>`,
|
||||
but the sampling of the correlated angle is performed as it was in the other
|
||||
two representations discussed in this sub-section. In the Law 61 algorithm,
|
||||
we found an interpolation factor :math:`f`, statistically sampled an incoming
|
||||
energy bin :math:`\ell`, and sampled an outgoing energy bin :math:`j` based on
|
||||
the tabulated cumulative distribution function. Once the outgoing energy has
|
||||
been determined with equation :eq:`ace-law-4-energy`, we then need to decide
|
||||
which angular distribution data to use. Like the linear-linear interpolation
|
||||
case in Law 61, the angular distribution closest to the sampled value of the
|
||||
cumulative distribution function for the outgoing energy is utilized. The
|
||||
actual algorithm utilized to sample the outgoing angle is shown in equation
|
||||
:eq:`inelastic-angle`.
|
||||
|
||||
.. _probability_tables:
|
||||
|
||||
----------------------------------------------
|
||||
|
|
|
|||
|
|
@ -78,7 +78,7 @@ Prerequisites
|
|||
./configure --prefix=/opt/hdf5/1.8.11-gnu --enable-fortran \
|
||||
--enable-fortran2003 --enable-parallel
|
||||
|
||||
You may omit '--enable-parallel' if you want to compile HDF5_ in serial.
|
||||
You may omit ``--enable-parallel`` if you want to compile HDF5_ in serial.
|
||||
|
||||
* PETSc_ for CMFD acceleration
|
||||
|
||||
|
|
@ -93,7 +93,7 @@ Prerequisites
|
|||
--with-fortran-datatypes
|
||||
|
||||
The BLAS/LAPACK library is not required to be downloaded and can be linked
|
||||
explicitly (e.g., Intel MLK library).
|
||||
explicitly (e.g., Intel MKL library).
|
||||
|
||||
* git_ version control software for obtaining source code
|
||||
|
||||
|
|
|
|||
83
src/ace.F90
83
src/ace.F90
|
|
@ -1250,15 +1250,10 @@ contains
|
|||
integer :: NE_out ! number of outgoing energies
|
||||
integer :: NMU ! number of outgoing angles
|
||||
integer :: JXS4 ! location of elastic energy table
|
||||
integer(8), allocatable :: LOCC(:) ! Location of inelastic data
|
||||
|
||||
! read secondary energy mode for inelastic scattering and check
|
||||
! read secondary energy mode for inelastic scattering
|
||||
table % secondary_mode = NXS(7)
|
||||
if (table % secondary_mode /= SAB_SECONDARY_EQUAL .and. &
|
||||
table % secondary_mode /= SAB_SECONDARY_SKEWED) then
|
||||
message = "Unsupported secondary mode on S(a,b) table " // &
|
||||
trim(adjustl(table % name)) // ": " // to_str(table % secondary_mode)
|
||||
call fatal_error()
|
||||
end if
|
||||
|
||||
! read number of inelastic energies and allocate arrays
|
||||
NE_in = int(XSS(JXS(1)))
|
||||
|
|
@ -1276,29 +1271,67 @@ contains
|
|||
|
||||
! allocate space for outgoing energy/angle for inelastic
|
||||
! scattering
|
||||
NE_out = NXS(4)
|
||||
NMU = NXS(3) + 1
|
||||
table % n_inelastic_e_out = NE_out
|
||||
table % n_inelastic_mu = NMU
|
||||
allocate(table % inelastic_e_out(NE_out, NE_in))
|
||||
allocate(table % inelastic_mu(NMU, NE_out, NE_in))
|
||||
if (table % secondary_mode == SAB_SECONDARY_EQUAL .or. &
|
||||
table % secondary_mode == SAB_SECONDARY_SKEWED) then
|
||||
NMU = NXS(3) + 1
|
||||
table % n_inelastic_mu = NMU
|
||||
NE_out = NXS(4)
|
||||
table % n_inelastic_e_out = NE_out
|
||||
allocate(table % inelastic_e_out(NE_out, NE_in))
|
||||
allocate(table % inelastic_mu(NMU, NE_out, NE_in))
|
||||
else if (table % secondary_mode == SAB_SECONDARY_CONT) then
|
||||
NMU = NXS(3) - 1
|
||||
table % n_inelastic_mu = NMU
|
||||
allocate(table % inelastic_data(NE_in))
|
||||
allocate(LOCC(NE_in))
|
||||
! NE_out will be determined later
|
||||
end if
|
||||
|
||||
! read outgoing energy/angle distribution for inelastic scattering
|
||||
lc = JXS(3) - 1
|
||||
do i = 1, NE_in
|
||||
do j = 1, NE_out
|
||||
! read outgoing energy
|
||||
table % inelastic_e_out(j,i) = XSS(lc + 1)
|
||||
if (table % secondary_mode == SAB_SECONDARY_EQUAL .or. &
|
||||
table % secondary_mode == SAB_SECONDARY_SKEWED) then
|
||||
lc = JXS(3) - 1
|
||||
do i = 1, NE_in
|
||||
do j = 1, NE_out
|
||||
! read outgoing energy
|
||||
table % inelastic_e_out(j,i) = XSS(lc + 1)
|
||||
|
||||
! read outgoing angles for this outgoing energy
|
||||
do k = 1, NMU
|
||||
table % inelastic_mu(k,j,i) = XSS(lc + 1 + k)
|
||||
! read outgoing angles for this outgoing energy
|
||||
do k = 1, NMU
|
||||
table % inelastic_mu(k,j,i) = XSS(lc + 1 + k)
|
||||
end do
|
||||
|
||||
! advance pointer
|
||||
lc = lc + 1 + NMU
|
||||
end do
|
||||
|
||||
! advance pointer
|
||||
lc = lc + 1 + NMU
|
||||
end do
|
||||
end do
|
||||
else if (table % secondary_mode == SAB_SECONDARY_CONT) then
|
||||
! Get the location pointers to each Ein's DistEnergySAB data
|
||||
LOCC = get_int(NE_in)
|
||||
! Get the number of outgoing energies and allocate space accordingly
|
||||
do i = 1, NE_in
|
||||
NE_out = int(XSS(XSS_index + i - 1))
|
||||
table % inelastic_data(i) % n_e_out = NE_out
|
||||
allocate(table % inelastic_data(i) % e_out (NE_out))
|
||||
allocate(table % inelastic_data(i) % e_out_pdf (NE_out))
|
||||
allocate(table % inelastic_data(i) % e_out_cdf (NE_out))
|
||||
allocate(table % inelastic_data(i) % mu (NMU, NE_out))
|
||||
end do
|
||||
|
||||
! Now we can fill the inelastic_data(i) attributes
|
||||
do i = 1, NE_in
|
||||
XSS_index = LOCC(i)
|
||||
NE_out = table % inelastic_data(i) % n_e_out
|
||||
do j = 1, NE_out
|
||||
table % inelastic_data(i) % e_out(j) = XSS(XSS_index + 1)
|
||||
table % inelastic_data(i) % e_out_pdf(j) = XSS(XSS_index + 2)
|
||||
table % inelastic_data(i) % e_out_cdf(j) = XSS(XSS_index + 3)
|
||||
table % inelastic_data(i) % mu(:, j) = &
|
||||
XSS(XSS_index + 4: XSS_index + 4 + NMU - 1)
|
||||
XSS_index = XSS_index + 4 + NMU - 1
|
||||
end do
|
||||
end do
|
||||
end if
|
||||
|
||||
! read number of elastic energies and allocate arrays
|
||||
JXS4 = JXS(4)
|
||||
|
|
|
|||
|
|
@ -32,7 +32,7 @@ module ace_header
|
|||
type(Tab1) :: p_valid ! probability of law validity
|
||||
real(8), allocatable :: data(:) ! energy distribution data
|
||||
|
||||
! For reactions that may have multiple energy distributions such as (n.2n),
|
||||
! For reactions that may have multiple energy distributions such as (n,2n),
|
||||
! this pointer allows multiple laws to be stored
|
||||
type(DistEnergy), pointer :: next => null()
|
||||
|
||||
|
|
@ -171,6 +171,20 @@ module ace_header
|
|||
|
||||
end type Nuclide0K
|
||||
|
||||
!===============================================================================
|
||||
! DISTENERGYSAB contains the secondary energy/angle distributions for inelastic
|
||||
! thermal scattering collisions which utilize a continuous secondary energy
|
||||
! representation.
|
||||
!===============================================================================
|
||||
|
||||
type DistEnergySab
|
||||
integer :: n_e_out
|
||||
real(8), allocatable :: e_out(:)
|
||||
real(8), allocatable :: e_out_pdf(:)
|
||||
real(8), allocatable :: e_out_cdf(:)
|
||||
real(8), allocatable :: mu(:,:)
|
||||
end type DistEnergySab
|
||||
|
||||
!===============================================================================
|
||||
! SALPHABETA contains S(a,b) data for thermal neutron scattering, typically off
|
||||
! of light isotopes such as water, graphite, Be, etc
|
||||
|
|
@ -191,11 +205,17 @@ module ace_header
|
|||
integer :: n_inelastic_e_in ! # of incoming E for inelastic
|
||||
integer :: n_inelastic_e_out ! # of outgoing E for inelastic
|
||||
integer :: n_inelastic_mu ! # of outgoing angles for inelastic
|
||||
integer :: secondary_mode ! secondary mode (equal/skewed)
|
||||
integer :: secondary_mode ! secondary mode (equal/skewed/continuous)
|
||||
real(8), allocatable :: inelastic_e_in(:)
|
||||
real(8), allocatable :: inelastic_sigma(:)
|
||||
! The following are used only if secondary_mode is 0 or 1
|
||||
real(8), allocatable :: inelastic_e_out(:,:)
|
||||
real(8), allocatable :: inelastic_mu(:,:,:)
|
||||
! The following is used only if secondary_mode is 3
|
||||
! The different implementation is necessary because the continuous
|
||||
! representation has a variable number of outgoing energy points for each
|
||||
! incoming energy
|
||||
type(DistEnergySab), allocatable :: inelastic_data(:) ! One for each Ein
|
||||
|
||||
! Elastic scattering data
|
||||
integer :: elastic_mode ! elastic mode (discrete/exact)
|
||||
|
|
|
|||
|
|
@ -155,7 +155,8 @@ module constants
|
|||
! Secondary energy mode for S(a,b) inelastic scattering
|
||||
integer, parameter :: &
|
||||
SAB_SECONDARY_EQUAL = 0, & ! Equally-likely outgoing energy bins
|
||||
SAB_SECONDARY_SKEWED = 1 ! Skewed outgoing energy bins
|
||||
SAB_SECONDARY_SKEWED = 1, & ! Skewed outgoing energy bins
|
||||
SAB_SECONDARY_CONT = 2 ! Continuous, linear-linear interpolation
|
||||
|
||||
! Elastic mode for S(a,b) elastic scattering
|
||||
integer, parameter :: &
|
||||
|
|
|
|||
|
|
@ -477,7 +477,8 @@ contains
|
|||
allocate(temp_int_array(n_tracks))
|
||||
call get_node_array(doc, "track", temp_int_array)
|
||||
|
||||
! Reshape into track_identifiers -- note automatic array allocation
|
||||
! Reshape into track_identifiers
|
||||
allocate(track_identifiers(3, n_tracks/3))
|
||||
track_identifiers = reshape(temp_int_array, [3, n_tracks/3])
|
||||
end if
|
||||
|
||||
|
|
@ -3011,7 +3012,19 @@ contains
|
|||
end if
|
||||
|
||||
! set filetype, record length, and number of entries
|
||||
listing % filetype = filetype
|
||||
if (check_for_node(node_ace, "filetype")) then
|
||||
temp_str = ''
|
||||
call get_node_value(node_ace, "filetype", temp_str)
|
||||
if (temp_str == 'ascii') then
|
||||
listing % filetype = ASCII
|
||||
else if (temp_str == 'binary') then
|
||||
listing % filetype = BINARY
|
||||
end if
|
||||
else
|
||||
listing % filetype = filetype
|
||||
end if
|
||||
|
||||
! Set record length and entries for binary files
|
||||
if (filetype == BINARY) then
|
||||
listing % recl = recl
|
||||
listing % entries = entries
|
||||
|
|
|
|||
225
src/physics.F90
225
src/physics.F90
|
|
@ -49,7 +49,7 @@ contains
|
|||
if (verbosity >= 10 .or. trace) then
|
||||
message = " " // trim(reaction_name(p % event_MT)) // " with " // &
|
||||
trim(adjustl(nuclides(p % event_nuclide) % name)) // &
|
||||
". Energy = " // trim(to_str(p % E * 1e6_8)) // " eV."
|
||||
". Energy = " // trim(to_str(p % E * 1e6_8)) // " eV."
|
||||
call write_message()
|
||||
end if
|
||||
|
||||
|
|
@ -160,7 +160,7 @@ contains
|
|||
call fatal_error()
|
||||
end if
|
||||
|
||||
! Find atom density
|
||||
! Find atom density
|
||||
i_nuclide = mat % nuclide(i)
|
||||
atom_density = mat % atom_density(i)
|
||||
|
||||
|
|
@ -211,7 +211,7 @@ contains
|
|||
i_reaction = nuc % index_fission(1)
|
||||
return
|
||||
end if
|
||||
|
||||
|
||||
! Get grid index and interpolatoin factor and sample fission cdf
|
||||
i_grid = micro_xs(i_nuclide) % index_grid
|
||||
f = micro_xs(i_nuclide) % interp_factor
|
||||
|
|
@ -228,7 +228,7 @@ contains
|
|||
if (i_grid < rxn % threshold) cycle
|
||||
|
||||
! add to cumulative probability
|
||||
prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) &
|
||||
prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) &
|
||||
+ f*(rxn%sigma(i_grid - rxn%threshold + 2)))
|
||||
|
||||
! Create fission bank sites if fission occus
|
||||
|
|
@ -384,7 +384,7 @@ contains
|
|||
if (i_grid < rxn % threshold) cycle
|
||||
|
||||
! add to cumulative probability
|
||||
prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) &
|
||||
prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) &
|
||||
+ f*(rxn%sigma(i_grid - rxn%threshold + 2)))
|
||||
end do
|
||||
|
||||
|
|
@ -494,13 +494,24 @@ contains
|
|||
integer :: k ! outgoing cosine bin
|
||||
integer :: n_energy_out ! number of outgoing energy bins
|
||||
real(8) :: f ! interpolation factor
|
||||
real(8) :: r ! used for skewed sampling
|
||||
real(8) :: r ! used for skewed sampling & continuous
|
||||
real(8) :: E_ij ! outgoing energy j for E_in(i)
|
||||
real(8) :: E_i1j ! outgoing energy j for E_in(i+1)
|
||||
real(8) :: mu_ijk ! outgoing cosine k for E_in(i) and E_out(j)
|
||||
real(8) :: mu_i1jk ! outgoing cosine k for E_in(i+1) and E_out(j)
|
||||
real(8) :: prob ! probability for sampling Bragg edge
|
||||
type(SAlphaBeta), pointer, save :: sab => null()
|
||||
! Following are needed only for SAB_SECONDARY_CONT scattering
|
||||
integer :: l ! sampled incoming E bin (is i or i + 1)
|
||||
real(8) :: E_i_1, E_i_J ! endpoints on outgoing grid i
|
||||
real(8) :: E_i1_1, E_i1_J ! endpoints on outgoing grid i+1
|
||||
real(8) :: E_1, E_J ! endpoints interpolated between i and i+1
|
||||
real(8) :: E_l_j, E_l_j1 ! adjacent E on outgoing grid l
|
||||
real(8) :: p_l_j, p_l_j1 ! adjacent p on outgoing grid l
|
||||
real(8) :: c_j, c_j1 ! cumulative probability
|
||||
real(8) :: frac ! interpolation factor on outgoing energy
|
||||
real(8) :: r1 ! RNG for outgoing energy
|
||||
|
||||
!$omp threadprivate(sab)
|
||||
|
||||
! Get pointer to S(a,b) table
|
||||
|
|
@ -517,7 +528,7 @@ contains
|
|||
f = ZERO
|
||||
else
|
||||
i = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, E)
|
||||
f = (E - sab%elastic_e_in(i)) / &
|
||||
f = (E - sab%elastic_e_in(i)) / &
|
||||
(sab%elastic_e_in(i+1) - sab%elastic_e_in(i))
|
||||
end if
|
||||
|
||||
|
|
@ -558,8 +569,7 @@ contains
|
|||
! Outgoing energy is same as incoming energy -- no need to do anything
|
||||
|
||||
else
|
||||
! Determine number of outgoing energy and angle bins
|
||||
n_energy_out = sab % n_inelastic_e_out
|
||||
! Perform inelastic calculations
|
||||
|
||||
! Get index and interpolation factor for inelastic grid
|
||||
if (E < sab % inelastic_e_in(1)) then
|
||||
|
|
@ -567,61 +577,150 @@ contains
|
|||
f = ZERO
|
||||
else
|
||||
i = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, E)
|
||||
f = (E - sab%inelastic_e_in(i)) / &
|
||||
f = (E - sab%inelastic_e_in(i)) / &
|
||||
(sab%inelastic_e_in(i+1) - sab%inelastic_e_in(i))
|
||||
end if
|
||||
|
||||
! Now that we have an incoming energy bin, we need to determine the
|
||||
! outgoing energy bin. This will depend on the "secondary energy
|
||||
! mode". If the mode is 0, then the outgoing energy bin is chosen from a
|
||||
! set of equally-likely bins. However, if the mode is 1, then the first
|
||||
! set of equally-likely bins. If the mode is 1, then the first
|
||||
! two and last two bins are skewed to have lower probabilities than the
|
||||
! other bins (0.1 for the first and last bins and 0.4 for the second and
|
||||
! second to last bins, relative to a normal bin probability of 1)
|
||||
! second to last bins, relative to a normal bin probability of 1).
|
||||
! Finally, if the mode is 2, then a continuous distribution (with
|
||||
! accompanying PDF and CDF is utilized)
|
||||
|
||||
if (sab % secondary_mode == SAB_SECONDARY_EQUAL) then
|
||||
! All bins equally likely
|
||||
j = 1 + int(prn() * n_energy_out)
|
||||
elseif (sab % secondary_mode == SAB_SECONDARY_SKEWED) then
|
||||
r = prn() * (n_energy_out - 3)
|
||||
if (r > ONE) then
|
||||
! equally likely N-4 middle bins
|
||||
j = int(r) + 2
|
||||
elseif (r > 0.6) then
|
||||
! second to last bin has relative probability of 0.4
|
||||
j = n_energy_out - 1
|
||||
elseif (r > 0.5) then
|
||||
! last bin has relative probability of 0.1
|
||||
j = n_energy_out
|
||||
elseif (r > 0.1) then
|
||||
! second bin has relative probability of 0.4
|
||||
j = 2
|
||||
else
|
||||
! first bin has relative probability of 0.1
|
||||
j = 1
|
||||
if ((sab % secondary_mode == SAB_SECONDARY_EQUAL) .or. &
|
||||
(sab % secondary_mode == SAB_SECONDARY_SKEWED)) then
|
||||
if (sab % secondary_mode == SAB_SECONDARY_EQUAL) then
|
||||
! All bins equally likely
|
||||
|
||||
j = 1 + int(prn() * sab % n_inelastic_e_out)
|
||||
elseif (sab % secondary_mode == SAB_SECONDARY_SKEWED) then
|
||||
! Distribution skewed away from edge points
|
||||
|
||||
! Determine number of outgoing energy and angle bins
|
||||
n_energy_out = sab % n_inelastic_e_out
|
||||
|
||||
r = prn() * (n_energy_out - 3)
|
||||
if (r > ONE) then
|
||||
! equally likely N-4 middle bins
|
||||
j = int(r) + 2
|
||||
elseif (r > 0.6) then
|
||||
! second to last bin has relative probability of 0.4
|
||||
j = n_energy_out - 1
|
||||
elseif (r > 0.5) then
|
||||
! last bin has relative probability of 0.1
|
||||
j = n_energy_out
|
||||
elseif (r > 0.1) then
|
||||
! second bin has relative probability of 0.4
|
||||
j = 2
|
||||
else
|
||||
! first bin has relative probability of 0.1
|
||||
j = 1
|
||||
end if
|
||||
end if
|
||||
|
||||
! Determine outgoing energy corresponding to E_in(i) and E_in(i+1)
|
||||
E_ij = sab % inelastic_e_out(j,i)
|
||||
E_i1j = sab % inelastic_e_out(j,i+1)
|
||||
|
||||
! Outgoing energy
|
||||
E = (1 - f)*E_ij + f*E_i1j
|
||||
|
||||
! Sample outgoing cosine bin
|
||||
k = 1 + int(prn() * sab % n_inelastic_mu)
|
||||
|
||||
! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1)
|
||||
mu_ijk = sab % inelastic_mu(k,j,i)
|
||||
mu_i1jk = sab % inelastic_mu(k,j,i+1)
|
||||
|
||||
! Cosine of angle between incoming and outgoing neutron
|
||||
mu = (1 - f)*mu_ijk + f*mu_i1jk
|
||||
|
||||
else if (sab % secondary_mode == SAB_SECONDARY_CONT) then
|
||||
! Continuous secondary energy - this is to be similar to
|
||||
! Law 61 interpolation on outgoing energy
|
||||
|
||||
! Sample between ith and (i+1)th bin
|
||||
r = prn()
|
||||
if (f > r) then
|
||||
l = i + 1
|
||||
else
|
||||
l = i
|
||||
end if
|
||||
|
||||
! Determine endpoints on grid i
|
||||
n_energy_out = sab % inelastic_data(i) % n_e_out
|
||||
E_i_1 = sab % inelastic_data(i) % e_out(1)
|
||||
E_i_J = sab % inelastic_data(i) % e_out(n_energy_out)
|
||||
|
||||
! Determine endpoints on grid i + 1
|
||||
n_energy_out = sab % inelastic_data(i + 1) % n_e_out
|
||||
E_i1_1 = sab % inelastic_data(i + 1) % e_out(1)
|
||||
E_i1_J = sab % inelastic_data(i + 1) % e_out(n_energy_out)
|
||||
|
||||
E_1 = E_i_1 + f * (E_i1_1 - E_i_1)
|
||||
E_J = E_i_J + f * (E_i1_J - E_i_J)
|
||||
|
||||
! Determine outgoing energy bin
|
||||
! (First reset n_energy_out to the right value)
|
||||
n_energy_out = sab % inelastic_data(l) % n_e_out
|
||||
r1 = prn()
|
||||
c_j = sab % inelastic_data(l) % e_out_cdf(1)
|
||||
do j = 1, n_energy_out - 1
|
||||
c_j1 = sab % inelastic_data(l) % e_out_cdf(j + 1)
|
||||
if (r1 < c_j1) exit
|
||||
c_j = c_j1
|
||||
end do
|
||||
|
||||
! check to make sure k is <= n_energy_out - 1
|
||||
j = min(j, n_energy_out - 1)
|
||||
|
||||
! Get the data to interpolate between
|
||||
E_l_j = sab % inelastic_data(l) % e_out(j)
|
||||
p_l_j = sab % inelastic_data(l) % e_out_pdf(j)
|
||||
|
||||
! Next part assumes linear-linear interpolation in standard
|
||||
E_l_j1 = sab % inelastic_data(l) % e_out(j + 1)
|
||||
p_l_j1 = sab % inelastic_data(l) % e_out_pdf(j + 1)
|
||||
|
||||
! Find secondary energy (variable E)
|
||||
frac = (p_l_j1 - p_l_j) / (E_l_j1 - E_l_j)
|
||||
if (frac == ZERO) then
|
||||
E = E_l_j + (r1 - c_j) / p_l_j
|
||||
else
|
||||
E = E_l_j + (sqrt(max(ZERO, p_l_j * p_l_j + &
|
||||
TWO * frac * (r1 - c_j))) - p_l_j) / frac
|
||||
end if
|
||||
|
||||
! Now interpolate between incident energy bins i and i + 1
|
||||
if (l == i) then
|
||||
E = E_1 + (E - E_i_1) * (E_J - E_1) / (E_i_J - E_i_1)
|
||||
else
|
||||
E = E_1 + (E - E_i1_1) * (E_J - E_1) / (E_i1_J - E_i1_1)
|
||||
end if
|
||||
|
||||
! Find angular distribution for closest outgoing energy bin
|
||||
if (r1 - c_j < c_j1 - r1) then
|
||||
j = j
|
||||
else
|
||||
j = j + 1
|
||||
end if
|
||||
|
||||
! Sample outgoing cosine bin
|
||||
k = 1 + int(prn() * sab % n_inelastic_mu)
|
||||
|
||||
! Will use mu from the randomly chosen incoming and closest outgoing
|
||||
! energy bins
|
||||
mu = sab % inelastic_data(l) % mu(k, j)
|
||||
|
||||
else
|
||||
message = "Invalid secondary energy mode on S(a,b) table " // &
|
||||
trim(sab % name)
|
||||
end if
|
||||
|
||||
! Determine outgoing energy corresponding to E_in(i) and E_in(i+1)
|
||||
E_ij = sab % inelastic_e_out(j,i)
|
||||
E_i1j = sab % inelastic_e_out(j,i+1)
|
||||
|
||||
! Outgoing energy
|
||||
E = (1 - f)*E_ij + f*E_i1j
|
||||
|
||||
! Sample outgoing cosine bin
|
||||
k = 1 + int(prn() * sab % n_inelastic_mu)
|
||||
|
||||
! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1)
|
||||
mu_ijk = sab % inelastic_mu(k,j,i)
|
||||
mu_i1jk = sab % inelastic_mu(k,j,i+1)
|
||||
|
||||
! Cosine of angle between incoming and outgoing neutron
|
||||
mu = (1 - f)*mu_ijk + f*mu_i1jk
|
||||
end if
|
||||
end if ! (inelastic secondary energy treatment)
|
||||
end if ! (elastic or inelastic)
|
||||
|
||||
! change direction of particle
|
||||
uvw = rotate_angle(uvw, mu)
|
||||
|
|
@ -1111,7 +1210,7 @@ contains
|
|||
E_cm = E
|
||||
|
||||
! determine outgoing energy in lab
|
||||
E = E_cm + (E_in + TWO * mu * (A+ONE) * sqrt(E_in * E_cm)) &
|
||||
E = E_cm + (E_in + TWO * mu * (A+ONE) * sqrt(E_in * E_cm)) &
|
||||
/ ((A+ONE)*(A+ONE))
|
||||
|
||||
! determine outgoing angle in lab
|
||||
|
|
@ -1174,7 +1273,7 @@ contains
|
|||
r = ONE
|
||||
else
|
||||
i = binary_search(rxn % adist % energy, n, E)
|
||||
r = (E - rxn % adist % energy(i)) / &
|
||||
r = (E - rxn % adist % energy(i)) / &
|
||||
(rxn % adist % energy(i+1) - rxn % adist % energy(i))
|
||||
end if
|
||||
|
||||
|
|
@ -1303,7 +1402,7 @@ contains
|
|||
end if
|
||||
|
||||
end function rotate_angle
|
||||
|
||||
|
||||
!===============================================================================
|
||||
! SAMPLE_ENERGY samples an outgoing energy distribution, either for a secondary
|
||||
! neutron from a collision or for a prompt/delayed fission neutron
|
||||
|
|
@ -1459,7 +1558,7 @@ contains
|
|||
! =======================================================================
|
||||
! CONTINUOUS TABULAR DISTRIBUTION
|
||||
|
||||
! read number of interpolation regions and incoming energies
|
||||
! read number of interpolation regions and incoming energies
|
||||
NR = int(edist % data(1))
|
||||
NE = int(edist % data(2 + 2*NR))
|
||||
if (NR == 1) then
|
||||
|
|
@ -1485,7 +1584,7 @@ contains
|
|||
r = ONE
|
||||
else
|
||||
i = binary_search(edist % data(lc+1:lc+NE), NE, E_in)
|
||||
r = (E_in - edist%data(lc+i)) / &
|
||||
r = (E_in - edist%data(lc+i)) / &
|
||||
(edist%data(lc+i+1) - edist%data(lc+i))
|
||||
end if
|
||||
|
||||
|
|
@ -1589,7 +1688,7 @@ contains
|
|||
! =======================================================================
|
||||
! MAXWELL FISSION SPECTRUM
|
||||
|
||||
! read number of interpolation regions and incoming energies
|
||||
! read number of interpolation regions and incoming energies
|
||||
NR = int(edist % data(1))
|
||||
NE = int(edist % data(2 + 2*NR))
|
||||
|
||||
|
|
@ -1621,7 +1720,7 @@ contains
|
|||
! =======================================================================
|
||||
! EVAPORATION SPECTRUM
|
||||
|
||||
! read number of interpolation regions and incoming energies
|
||||
! read number of interpolation regions and incoming energies
|
||||
NR = int(edist % data(1))
|
||||
NE = int(edist % data(2 + 2*NR))
|
||||
|
||||
|
|
@ -1702,7 +1801,7 @@ contains
|
|||
call fatal_error()
|
||||
end if
|
||||
|
||||
! read number of interpolation regions and incoming energies
|
||||
! read number of interpolation regions and incoming energies
|
||||
NR = int(edist % data(1))
|
||||
NE = int(edist % data(2 + 2*NR))
|
||||
if (NR > 0) then
|
||||
|
|
@ -1724,7 +1823,7 @@ contains
|
|||
r = ONE
|
||||
else
|
||||
i = binary_search(edist % data(lc+1:lc+NE), NE, E_in)
|
||||
r = (E_in - edist%data(lc+i)) / &
|
||||
r = (E_in - edist%data(lc+i)) / &
|
||||
(edist%data(lc+i+1) - edist%data(lc+i))
|
||||
end if
|
||||
|
||||
|
|
@ -1851,11 +1950,11 @@ contains
|
|||
|
||||
if (.not. present(mu_out)) then
|
||||
! call write_particle_restart()
|
||||
message = "Law 44 called without giving mu_out as argument."
|
||||
message = "Law 61 called without giving mu_out as argument."
|
||||
call fatal_error()
|
||||
end if
|
||||
|
||||
! read number of interpolation regions and incoming energies
|
||||
! read number of interpolation regions and incoming energies
|
||||
NR = int(edist % data(1))
|
||||
NE = int(edist % data(2 + 2*NR))
|
||||
if (NR > 0) then
|
||||
|
|
@ -1877,7 +1976,7 @@ contains
|
|||
r = ONE
|
||||
else
|
||||
i = binary_search(edist % data(lc+1:lc+NE), NE, E_in)
|
||||
r = (E_in - edist%data(lc+i)) / &
|
||||
r = (E_in - edist%data(lc+i)) / &
|
||||
(edist%data(lc+i+1) - edist%data(lc+i))
|
||||
end if
|
||||
|
||||
|
|
|
|||
|
|
@ -10,7 +10,9 @@ element cross_sections {
|
|||
(element temperature { xsd:double } | attribute temperature { xsd:double }) &
|
||||
(element path { xsd:string { maxLength = "255" } } |
|
||||
attribute path { xsd:string { maxLength = "255" } }) &
|
||||
(element location { xsd:int } | attribute location { xsd:int })?
|
||||
(element location { xsd:int } | attribute location { xsd:int })? &
|
||||
(element filetype { ( "ascii" | "binary" ) } |
|
||||
attribute filetype { ( "ascii" | "binary" ) })?
|
||||
}* &
|
||||
|
||||
element directory { xsd:string { maxLength = "255" } }? &
|
||||
|
|
|
|||
|
|
@ -97,6 +97,7 @@ contains
|
|||
|
||||
! Fill p with needed data
|
||||
p % coord0 % xyz = site % xyz
|
||||
p % coord0 % uvw = [ ONE, ZERO, ZERO ]
|
||||
|
||||
! Now search to see if location exists in geometry
|
||||
call find_cell(p, found)
|
||||
|
|
|
|||
|
|
@ -71,7 +71,7 @@ class Xsdir(object):
|
|||
# Handle continuation lines
|
||||
while words[-1] == '+':
|
||||
extraWords = self.f.readline().split()
|
||||
words = words + extraWords
|
||||
words = words[:-1] + extraWords
|
||||
assert len(words) >= 7
|
||||
|
||||
# Create XsdirTable object and add to line
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue