crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-26 21:25:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

added treatment of delayed neutrons to the multi-group solver

Browse source
This commit is contained in:
Sam Shaner 2016-10-13 14:38:06 -04:00
parent 9035f6951e
commit 89149ea7eb
15 changed files with 3019 additions and 805 deletions
Download patch file Download diff file Expand all files Collapse all files

12
docs/source/io_formats/mgxs_library.rst
View file

@ -15,8 +15,9 @@ to be used in the ``mgxs.h5`` file.
The multi-group library is provided in the HDF5_ format. This library must
provide some meta-data about the library itself (such as the number of
groups and the group structure, etc.) as well as the actual cross section
data itself for each of the necessary nuclides or materials.
energy groups, delayed groups, and the energy group structure, etc.) as
well as the actual cross section data itself for each of the necessary
nuclides or materials.
.. _HDF5: http://www.hdfgroup.org/HDF5/
@ -28,7 +29,8 @@ MGXS Library Specification
**/**
:Attributes: - **groups** (*int*) -- Number of energy groups
:Attributes: - **energy_groups** (*int*) -- Number of energy groups
- **delayed_groups** (*int*) -- Number of delayed groups (optional)
- **group structure** (*double[]*) -- Monotonically increasing
list of group boundaries, in units of MeV. The length of this
array should be the number of groups plus 1.
@ -110,7 +112,7 @@ Temperature-dependent data, provided for temperature <TTT>K.
spectra as well and thus will have one additional dimension
for the outgoing energy group. In this case, `nu-fission` has the
same dimensionality as `multiplicity matrix`.
- **inverse velocities** (*double[]*) -- Average inverse velocity
- **inverse-velocity** (*double[]*) -- Average inverse velocity
for each of the groups in the library. This dataset is optional.
**/<library name>/<TTT>K/scatter data/**
@ -159,4 +161,4 @@ Data specific to neutron scattering for the temperature <TTT>K
The former is used when `representation` is "isotropic", and the
latter when `representation` is "angle". Finally, the g_out
dimension has a dimensionality of `g_out bounds`[0] to
`g_out bounds`[1].
`g_out bounds`[1].

12
docs/source/usersguide/input.rst
View file

@ -1681,8 +1681,6 @@ The ``<tally>`` element accepts the following sub-elements:
<filter type="delayedgroup" bins="1 2 3 4 5 6" />
.. note:: This filter type is not used in the multi-group :ref:`energy_mode`.
:nuclides:
If specified, the scores listed will be for particular nuclides, not the
summation of reactions from all nuclides. The format for nuclides should be
@ -1863,12 +1861,10 @@ The ``<tally>`` element accepts the following sub-elements:
|Score | Description |
+======================+===================================================+
|delayed-nu-fission |Total production of delayed neutrons due to |
| |fission. This score type is not used in the |
| |multi-group :ref:`energy_mode`. |
| |fission. |
+----------------------+---------------------------------------------------+
|prompt-nu-fission |Total production of prompt neutrons due to |
| |fission. This score type is not used in the |
| |multi-group :ref:`energy_mode`. |
| |fission. |
+----------------------+---------------------------------------------------+
|nu-fission |Total production of neutrons due to fission. |
+----------------------+---------------------------------------------------+
@ -1897,8 +1893,6 @@ The ``<tally>`` element accepts the following sub-elements:
+----------------------+---------------------------------------------------+
|inverse-velocity |The flux-weighted inverse velocity where the |
| |velocity is in units of centimeters per second. |
| |This score type is not used in the |
| |multi-group :ref:`energy_mode`. |
+----------------------+---------------------------------------------------+
|kappa-fission |The recoverable energy production rate due to |
| |fission. The recoverable energy is defined as the |
@ -1934,8 +1928,6 @@ The ``<tally>`` element accepts the following sub-elements:
+----------------------+---------------------------------------------------+
|decay-rate |The delayed-nu-fission-weighted decay rate where |
| |the decay rate is in units of inverse seconds. |
| |This score type is not used in the |
| |multi-group :ref:`energy_mode`. |
+----------------------+---------------------------------------------------+
.. note::

6
examples/python/pincell_multigroup/build-xml.py
View file

@ -19,7 +19,7 @@ groups = openmc.mgxs.EnergyGroups(group_edges=[1E-11, 0.0635E-6, 10.0E-6,
1.0E-4, 1.0E-3, 0.5, 1.0, 20.0])
# Instantiate the 7-group (C5G7) cross section data
uo2_xsdata = openmc.XSdata('UO2', groups)
uo2_xsdata = openmc.XSdata('UO2', groups, delayed_groups=8)
uo2_xsdata.order = 0
uo2_xsdata.set_total(
[0.1779492, 0.3298048, 0.4803882, 0.5543674, 0.3118013, 0.3951678,
@ -42,6 +42,8 @@ uo2_xsdata.set_nu_fission([2.005998E-02, 2.027303E-03, 1.570599E-02,
5.257105E-01])
uo2_xsdata.set_chi([5.8791E-01, 4.1176E-01, 3.3906E-04, 1.1761E-07, 0.0000E+00,
0.0000E+00, 0.0000E+00])
uo2_xsdata.set_beta([2.13333e-4, 1.04514e-3, 6.03969e-4, 1.33963e-3,
2.29386e-3, 7.05174e-4, 6.00381e-4, 2.07736e-4])
h2o_xsdata = openmc.XSdata('LWTR', groups)
h2o_xsdata.order = 0
@ -59,7 +61,7 @@ h2o_xsdata.set_scatter_matrix(
[0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0022157, 0.6999130, 0.5373200],
[0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.1324400, 2.4807000]]])
mg_cross_sections_file = openmc.MGXSLibrary(groups)
mg_cross_sections_file = openmc.MGXSLibrary(groups, delayed_groups=8)
mg_cross_sections_file.add_xsdatas([uo2_xsdata, h2o_xsdata])
mg_cross_sections_file.export_to_hdf5()

79
openmc/mgxs/library.py
View file

@ -898,6 +898,7 @@ class Library(object):
if nuclide is not 'total':
name += '_' + nuclide
xsdata = openmc.XSdata(name, self.energy_groups)
xsdata.delayed_groups = self.num_delayed_groups
if order is None:
# Set the order to the Library's order (the defualt behavior)
@ -924,41 +925,79 @@ class Library(object):
mymgxs = self.get_mgxs(domain, 'nu-transport')
xsdata.set_total_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide],
subdomains=subdomain)
elif 'total' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'total')
xsdata.set_total_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide],
subdomain=subdomain)
if 'absorption' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'absorption')
xsdata.set_absorption_mgxs(mymgxs, xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
if 'fission' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'fission')
xsdata.set_fission_mgxs(mymgxs, xs_type=xs_type,
nuclide=[nuclide], subdomain=subdomain)
if 'kappa-fission' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'kappa-fission')
xsdata.set_kappa_fission_mgxs(mymgxs, xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
# For chi and nu-fission we can either have only a nu-fission matrix
# provided, or vectors of chi and nu-fission provided
if 'inverse-velocity' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'inverse-velocity')
xsdata.set_inverse_velocity_mgxs(mymgxs, xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
if 'nu-fission matrix' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'nu-fission matrix')
xsdata.set_nu_fission_mgxs(mymgxs, xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
else:
if 'chi' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'chi')
xsdata.set_chi_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide],
subdomain=subdomain)
if 'nu-fission' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'nu-fission')
xsdata.set_nu_fission_mgxs(mymgxs, xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
if 'chi' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'chi')
xsdata.set_chi_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide],
subdomain=subdomain)
if 'chi-prompt' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'chi-prompt')
xsdata.set_chi_prompt_mgxs(mymgxs, xs_type=xs_type,
nuclide=[nuclide], subdomain=subdomain)
if 'chi-delayed' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'chi-delayed')
xsdata.set_chi_delayed_mgxs(mymgxs, xs_type=xs_type,
nuclide=[nuclide], subdomain=subdomain)
if 'nu-fission' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'nu-fission')
xsdata.set_nu_fission_mgxs(mymgxs, xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
if 'prompt-nu-fission' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'prompt-nu-fission')
xsdata.set_prompt_nu_fission_mgxs(mymgxs, xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
if 'delayed-nu-fission' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'delayed-nu-fission')
xsdata.set_delayed_nu_fission_mgxs(mymgxs, xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
if 'beta' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'nu-fission')
xsdata.set_beta_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide],
subdomain=subdomain)
# If multiplicity matrix is available, prefer that
if 'multiplicity matrix' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'multiplicity matrix')
@ -966,6 +1005,7 @@ class Library(object):
nuclide=[nuclide],
subdomain=subdomain)
using_multiplicity = True
# multiplicity will fall back to using scatter and nu-scatter
elif ((('scatter matrix' in self.mgxs_types) and
('nu-scatter matrix' in self.mgxs_types))):
@ -976,6 +1016,7 @@ class Library(object):
nuclide=[nuclide],
subdomain=subdomain)
using_multiplicity = True
else:
using_multiplicity = False
@ -991,14 +1032,13 @@ class Library(object):
nuclide=[nuclide],
subdomain=subdomain)
# Since we are not using multiplicity, then
# scattering multiplication (nu-scatter) must be
# accounted for approximately by using an adjusted
# absorption cross section.
# Since we are not using multiplicity, then scattering
# multiplication (nu-scatter) must be accounted for
# approximately by using an adjusted absorption cross section.
if 'total' in self.mgxs_types:
for i in range(len(xsdata.temperatures)):
xsdata._absorption[i] = \
np.subtract(xsdata._total[i], np.sum(
xsdata._absorption[i] \
= np.subtract(xsdata._total[i], np.sum(
xsdata._scatter_matrix[i][0, :, :], axis=1))
return xsdata
@ -1051,7 +1091,8 @@ class Library(object):
xs_type = 'macro'
# Initialize file
mgxs_file = openmc.MGXSLibrary(self.energy_groups)
mgxs_file = openmc.MGXSLibrary(self.energy_groups,
delayed_groups=self.num_delayed_groups)
if self.domain_type == 'mesh':
# Create the xsdata objects and add to the mgxs_file

934
openmc/mgxs_library.py
View file

File diff suppressed because it is too large Load diff

2
src/cmfd_input.F90
View file

@ -107,7 +107,7 @@ contains
if (.not. run_CE) then
do i = 1, ng
found = .false.
do g = 1, energy_groups + 1
do g = 1, num_energy_groups + 1
if (cmfd % egrid(i) == energy_bins(g)) then
found = .true.
exit

5
src/global.F90
View file

@ -115,7 +115,10 @@ module global
type(MgxsContainer), target, allocatable :: macro_xs(:)
! Number of energy groups
integer :: energy_groups
integer :: num_energy_groups
! Number of delayed groups
integer :: num_delayed_groups
! Energy group structure
real(8), allocatable :: energy_bins(:)

21
src/hdf5_interface.F90
View file

@ -99,6 +99,7 @@ module hdf5_interface
public :: open_dataset
public :: close_dataset
public :: get_shape
public :: get_ndims
public :: write_attribute_string
public :: get_groups
public :: get_datasets
@ -2601,6 +2602,26 @@ contains
end if
end subroutine get_shape
subroutine get_ndims(obj_id, ndims)
integer(HID_T), intent(in) :: obj_id
integer(HID_T), intent(out) :: ndims
integer :: hdf5_err
integer :: type
integer(HID_T) :: space_id
call h5iget_type_f(obj_id, type, hdf5_err)
if (type == H5I_DATASET_F) then
call h5dget_space_f(obj_id, space_id, hdf5_err)
call h5sget_simple_extent_ndims_f(space_id, ndims, hdf5_err)
call h5sclose_f(space_id, hdf5_err)
elseif (type == H5I_ATTR_F) then
call h5aget_space_f(obj_id, space_id, hdf5_err)
call h5sget_simple_extent_ndims_f(space_id, ndims, hdf5_err)
call h5sclose_f(space_id, hdf5_err)
end if
end subroutine get_ndims
function using_mpio_device(obj_id) result(mpio)
integer(HID_T), intent(in) :: obj_id
logical :: mpio

58
src/input_xml.F90
View file

@ -3060,8 +3060,8 @@ contains
! index is simply the group (after flipping for the different
! ordering of the library and tallying systems).
if (.not. run_CE) then
if (n_words == energy_groups + 1) then
if (all(filt % bins == energy_bins(energy_groups + 1:1:-1))) &
if (n_words == num_energy_groups + 1) then
if (all(filt % bins == energy_bins(num_energy_groups + 1:1:-1))) &
then
filt % matches_transport_groups = .true.
end if
@ -3086,8 +3086,8 @@ contains
! index is simply the group (after flipping for the different
! ordering of the library and tallying systems).
if (.not. run_CE) then
if (n_words == energy_groups + 1) then
if (all(filt % bins == energy_bins(energy_groups + 1:1:-1))) &
if (n_words == num_energy_groups + 1) then
if (all(filt % bins == energy_bins(num_energy_groups + 1:1:-1))) &
then
filt % matches_transport_groups = .true.
end if
@ -3101,14 +3101,6 @@ contains
t % estimator = ESTIMATOR_ANALOG
case ('delayedgroup')
! Check to see if running in MG mode, because if so, the current
! system isnt set up yet to support delayed group data and thus
! these tallies
if (.not. run_CE) then
call fatal_error("delayedgroup filter on tally " &
// trim(to_str(t % id)) // " not yet supported&
& for multi-group mode.")
end if
! Allocate and declare the filter type
allocate(DelayedGroupFilter::t % filters(j) % obj)
@ -3599,7 +3591,13 @@ contains
end if
case ('decay-rate')
t % score_bins(j) = SCORE_DECAY_RATE
t % estimator = ESTIMATOR_ANALOG
! Set tally estimator to analog for CE mode
! (MG mode has all data available without a collision being
! necessary)
if (run_CE) then
t % estimator = ESTIMATOR_ANALOG
end if
case ('delayed-nu-fission')
t % score_bins(j) = SCORE_DELAYED_NU_FISSION
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
@ -3612,12 +3610,6 @@ contains
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
end if
! Disallow for MG mode since data not present
if (.not. run_CE) then
call fatal_error("Cannot tally delayed nu-fission rate in &
&multi-group mode")
end if
case ('kappa-fission')
t % score_bins(j) = SCORE_KAPPA_FISSION
case ('inverse-velocity')
@ -4654,15 +4646,25 @@ contains
! Open file for reading
file_id = file_open(path_cross_sections, 'r', parallel=.true.)
if (check_attribute(file_id, "groups")) then
! Get neutron group count
call read_attribute(energy_groups, file_id, "groups")
if (check_attribute(file_id, "energy_groups")) then
! Get neutron energy group count
call read_attribute(num_energy_groups, file_id, "energy_groups")
else
call fatal_error("'groups' attribute must exist!")
call fatal_error("'energy_groups' attribute must exist!")
end if
allocate(rev_energy_bins(energy_groups + 1))
allocate(energy_bins(energy_groups + 1))
if (check_attribute(file_id, "delayed_groups")) then
! Get neutron delayed group count
call read_attribute(num_delayed_groups, file_id, "delayed_groups")
else
num_delayed_groups = 0
call write_message("WARNING: delayed_groups element not provided so &
&number of delayed groups set to 0")
end if
allocate(rev_energy_bins(num_energy_groups + 1))
allocate(energy_bins(num_energy_groups + 1))
if (check_attribute(file_id, "group structure")) then
! Get neutron group structure
call read_attribute(energy_bins, file_id, "group structure")
@ -4671,10 +4673,10 @@ contains
end if
! First reverse the order of energy_groups
energy_bins = energy_bins(energy_groups + 1:1:-1)
energy_bins = energy_bins(num_energy_groups + 1:1:-1)
allocate(energy_bin_avg(energy_groups))
do i = 1, energy_groups
allocate(energy_bin_avg(num_energy_groups))
do i = 1, num_energy_groups
energy_bin_avg(i) = HALF * (energy_bins(i) + energy_bins(i + 1))
end do

88
src/mgxs_data.F90
View file

@ -21,26 +21,25 @@ contains
!===============================================================================
subroutine read_mgxs()
integer :: i ! index in materials array
integer :: j ! index over nuclides in material
integer :: i_xsdata ! index in xsdata_dict
integer :: i_nuclide ! index in nuclides array
character(20) :: name ! name of library to load
integer :: representation ! Data representation
integer :: i ! index in materials array
integer :: j ! index over nuclides in material
integer :: i_xsdata ! index in <xsdata> list
integer :: i_nuclide ! index in nuclides array
character(20) :: name ! name of library to load
integer :: representation ! Data representation
character(MAX_LINE_LEN) :: temp_str
type(Material), pointer :: mat
type(SetChar) :: already_read
integer(HID_T) :: file_id
integer(HID_T) :: xsdata_group
logical :: file_exists
logical :: get_kfiss, get_fiss
integer :: l
type(DictCharInt) :: xsdata_dict
type(SetChar) :: already_read
integer(HID_T) :: file_id
integer(HID_T) :: xsdata_group
logical :: file_exists
type(DictCharInt) :: xsdata_dict
type(VectorReal), allocatable :: temps(:)
! Check if MGXS Library exists
inquire(FILE=path_cross_sections, EXIST=file_exists)
if (.not. file_exists) then
! Could not find MGXS Library file
call fatal_error("Cross sections HDF5 file '" &
&// trim(path_cross_sections) // "' does not exist!")
@ -60,23 +59,6 @@ contains
allocate(micro_xs(n_nuclides_total))
!$omp end parallel
! Find out if we need fission & kappa fission
! (i.e., are there any SCORE_FISSION or SCORE_KAPPA_FISSION tallies?)
get_kfiss = .false.
get_fiss = .false.
do i = 1, n_tallies
do l = 1, tallies(i) % n_score_bins
if (tallies(i) % score_bins(l) == SCORE_KAPPA_FISSION) then
get_kfiss = .true.
end if
if (tallies(i) % score_bins(l) == SCORE_FISSION .or. &
tallies(i) % score_bins(l) == SCORE_NU_FISSION) then
get_fiss = .true.
end if
end do
if (get_kfiss .and. get_fiss) exit
end do
! ==========================================================================
! READ ALL MGXS CROSS SECTION TABLES
@ -103,7 +85,9 @@ contains
! First find out the data representation
if (check_attribute(xsdata_group, "representation")) then
call read_attribute(temp_str, xsdata_group, "representation")
if (trim(temp_str) == 'isotropic') then
representation = MGXS_ISOTROPIC
else if (trim(temp_str) == 'angle') then
@ -118,16 +102,19 @@ contains
! Now allocate accordingly
select case(representation)
case(MGXS_ISOTROPIC)
allocate(MgxsIso :: nuclides_MG(i_nuclide) % obj)
case(MGXS_ANGLE)
allocate(MgxsAngle :: nuclides_MG(i_nuclide) % obj)
end select
! Now read in the data specific to the type we just declared
call nuclides_MG(i_nuclide) % obj % from_hdf5(xsdata_group, &
energy_groups, temps(i_nuclide), temperature_method, &
temperature_tolerance, get_kfiss, get_fiss, max_order, &
num_energy_groups, num_delayed_groups, temps(i_nuclide), &
temperature_method, temperature_tolerance, max_order, &
legendre_to_tabular, legendre_to_tabular_points)
! Add name to dictionary
@ -147,6 +134,7 @@ contains
! Loop around nuclides in material
NUCLIDE_LOOP2: do j = 1, mat % n_nuclides
! Is this fissionable?
if (nuclides_MG(mat % nuclide(j)) % obj % fissionable) then
mat % fissionable = .true.
@ -161,12 +149,12 @@ contains
end subroutine read_mgxs
!===============================================================================
! CREATE_MACRO_XS generates the macroscopic x/s from the microscopic input data
! CREATE_MACRO_XS generates the macroscopic xs from the microscopic input data
!===============================================================================
subroutine create_macro_xs()
integer :: i_mat ! index in materials array
type(Material), pointer :: mat ! current material
integer :: i_mat ! index in materials array
type(Material), pointer :: mat ! current material
type(VectorReal), allocatable :: kTs(:)
allocate(macro_xs(n_materials))
@ -174,27 +162,32 @@ contains
! Get temperatures to read for each material
call get_mat_kTs(kTs)
! Force all nuclides in a material to be the same representation.
! Therefore type(nuclides(mat % nuclide(1)) % obj) dictates type(macroxs).
! At the same time, we will find the scattering type, as that will dictate
! how we allocate the scatter object within macroxs.allocate(macro_xs(n_materials))
do i_mat = 1, n_materials
! Get the material
mat => materials(i_mat)
! Check to see how our nuclides are represented
! Force all to be the same type
! Therefore type(nuclides(mat % nuclide(1)) % obj) dictates type(macroxs)
! Get the scattering type for the first nuclide
select type(nuc => nuclides_MG(mat % nuclide(1)) % obj)
type is (MgxsIso)
allocate(MgxsIso :: macro_xs(i_mat) % obj)
type is (MgxsAngle)
allocate(MgxsAngle :: macro_xs(i_mat) % obj)
end select
! Do not read materials which we do not actually use in the problem to
! save space
if (allocated(kTs(i_mat) % data)) then
call macro_xs(i_mat) % obj % combine(kTs(i_mat), mat, nuclides_MG, &
energy_groups, max_order, &
temperature_tolerance, &
temperature_method)
num_energy_groups, num_delayed_groups, max_order, &
temperature_tolerance, temperature_method)
end if
end do
end subroutine create_macro_xs
!===============================================================================
@ -203,16 +196,17 @@ contains
!===============================================================================
subroutine get_mat_kTs(kTs)
type(VectorReal), allocatable, intent(out) :: kTs(:)
integer :: i, j
type(VectorReal), allocatable, intent(out) :: kTs(:)
integer :: i, j ! Cell and material index
integer :: i_material ! Index in materials array
real(8) :: kT ! temperature in MeV
real(8) :: kT ! temperature in MeV
allocate(kTs(size(materials)))
do i = 1, size(cells)
do j = 1, size(cells(i) % material)
! Skip any non-material cells and void materials
if (cells(i) % material(j) == NONE .or. &
cells(i) % material(j) == MATERIAL_VOID) cycle
@ -243,17 +237,18 @@ contains
!===============================================================================
subroutine get_temperatures(temps)
type(VectorReal), allocatable, intent(out) :: temps(:)
type(VectorReal), allocatable, intent(out) :: temps(:)
integer :: i, j, k
integer :: i_nuclide ! index in nuclides array
integer :: i_material
integer :: i_material ! Index in materials array
real(8) :: temperature ! temperature in Kelvin
allocate(temps(n_nuclides_total))
do i = 1, size(cells)
do j = 1, size(cells(i) % material)
! Skip any non-material cells and void materials
if (cells(i) % material(j) == NONE .or. &
cells(i) % material(j) == MATERIAL_VOID) cycle
@ -268,6 +263,7 @@ contains
i_material = material_dict % get_key(cells(i) % material(j))
associate (mat => materials(i_material))
NUC_NAMES_LOOP: do k = 1, size(mat % names)
! Get index in temps array
i_nuclide = nuclide_dict % get_key(to_lower(mat % names(k)))

1975
src/mgxs_header.F90
View file

File diff suppressed because it is too large Load diff

26
src/physics_mg.F90
View file

@ -171,7 +171,10 @@ contains
type(Bank), intent(inout) :: bank_array(:)
integer(8), intent(inout) :: size_bank
integer :: nu_d(MAX_DELAYED_GROUPS) ! number of delayed neutrons born
integer :: i ! loop index
integer :: dg ! delayed group
integer :: gout ! group out
integer :: nu ! actual number of neutrons produced
integer :: ijk(3) ! indices in ufs mesh
real(8) :: nu_t ! total nu
@ -190,8 +193,10 @@ contains
! the expected number of fission sites produced
if (ufs) then
! Determine indices on ufs mesh for current location
call get_mesh_indices(ufs_mesh, p % coord(1) % xyz, ijk, in_mesh)
if (.not. in_mesh) then
call write_particle_restart(p)
call fatal_error("Source site outside UFS mesh!")
@ -209,6 +214,7 @@ contains
! Determine expected number of neutrons produced
nu_t = p % wgt / keff * weight * &
material_xs % nu_fission / material_xs % total
! Sample number of neutrons produced
if (prn() > nu_t - int(nu_t)) then
nu = int(nu_t)
@ -234,6 +240,10 @@ contains
! Bank source neutrons
if (nu == 0 .or. size_bank == size(bank_array)) return
! Initialize counter of delayed neutrons encountered for each delayed group
! to zero.
nu_d(:) = 0
p % fission = .true. ! Fission neutrons will be banked
do i = int(size_bank,4) + 1, int(min(size_bank + nu, int(size(bank_array),8)),4)
! Bank source neutrons by copying particle data
@ -254,8 +264,19 @@ contains
! Sample secondary energy distribution for fission reaction and set energy
! in fission bank
bank_array(i) % E = &
real(xs % sample_fission_energy(p % g, bank_array(i) % uvw), 8)
call xs % sample_fission_energy(p % g, bank_array(i) % uvw, dg, gout)
bank_array(i) % E = real(gout, 8)
bank_array(i) % delayed_group = dg
! Set delayed group on particle too
p % delayed_group = dg
! Increment the number of neutrons born delayed
if (p % delayed_group > 0) then
nu_d(p % delayed_group) = nu_d(p % delayed_group) + 1
end if
end do
! increment number of bank sites
@ -264,6 +285,7 @@ contains
! Store total weight banked for analog fission tallies
p % n_bank = nu
p % wgt_bank = nu/weight
p % n_delayed_bank(:) = nu_d(:)
end subroutine create_fission_sites

6
src/source.F90
View file

@ -199,10 +199,10 @@ contains
if (.not. run_CE) then
if (site % E <= energy_bins(1)) then
site % E = real(1, 8)
else if (site % E > energy_bins(energy_groups + 1)) then
site % E = real(energy_groups, 8)
else if (site % E > energy_bins(num_energy_groups + 1)) then
site % E = real(num_energy_groups, 8)
else
site % E = real(binary_search(energy_bins, energy_groups + 1, &
site % E = real(binary_search(energy_bins, num_energy_groups + 1, &
site % E), 8)
end if
end if

596
src/tally.F90
View file

@ -164,6 +164,7 @@ contains
case (SCORE_INVERSE_VELOCITY)
! make sure the correct energy is used
if (t % estimator == ESTIMATOR_TRACKLENGTH) then
E = p % E
@ -382,7 +383,7 @@ contains
! neutrons were emitted with different energies, multiple
! outgoing energy bins may have been scored to. The following
! logic treats this special case and results to multiple bins
call score_fission_eout_ce(p, t, score_index, score_bin)
call score_fission_eout(p, t, score_index, score_bin)
cycle SCORE_LOOP
end if
end if
@ -432,7 +433,7 @@ contains
! neutrons were emitted with different energies, multiple
! outgoing energy bins may have been scored to. The following
! logic treats this special case and results to multiple bins
call score_fission_eout_ce(p, t, score_index, score_bin)
call score_fission_eout(p, t, score_index, score_bin)
cycle SCORE_LOOP
end if
end if
@ -504,7 +505,7 @@ contains
! neutrons were emitted with different energies, multiple
! outgoing energy bins may have been scored to. The following
! logic treats this special case and results to multiple bins
call score_fission_eout_ce(p, t, score_index, score_bin)
call score_fission_eout(p, t, score_index, score_bin)
cycle SCORE_LOOP
end if
end if
@ -512,7 +513,8 @@ contains
! No fission events occur if survival biasing is on -- need to
! calculate fraction of absorptions that would have resulted in
! delayed-nu-fission
if (micro_xs(p % event_nuclide) % absorption > ZERO) then
if (micro_xs(p % event_nuclide) % absorption > ZERO .and. &
nuclides(p % event_nuclide) % fissionable) then
! Check if the delayed group filter is present
if (dg_filter > 0) then
@ -572,7 +574,9 @@ contains
d = filt % groups(d_bin)
! Compute the score and tally to bin
score = keff * p % wgt_bank / p % n_bank * p % n_delayed_bank(d)
score = keff * p % wgt_bank / p % n_bank * &
p % n_delayed_bank(d)
call score_fission_delayed_dg(t, d_bin, score, score_index)
end do
cycle SCORE_LOOP
@ -1106,6 +1110,11 @@ contains
integer :: q ! loop index for scoring bins
integer :: score_bin ! scoring bin, e.g. SCORE_FLUX
integer :: score_index ! scoring bin index
integer :: d ! delayed neutron index
integer :: g ! delayed neutron index
integer :: k ! loop index for bank sites
integer :: d_bin ! delayed group bin index
integer :: dg_filter ! index of delayed group filter
real(8) :: score ! analog tally score
real(8) :: p_uvw(3) ! Particle's current uvw
integer :: p_g ! Particle group to use for getting info
@ -1117,7 +1126,9 @@ contains
! this only depends on if we
if (t % estimator == ESTIMATOR_ANALOG .or. &
t % estimator == ESTIMATOR_COLLISION) then
if (survival_biasing) then
! Then we either are alive and had a scatter (and so g changed),
! or are dead and g did not change
if (p % alive) then
@ -1128,16 +1139,19 @@ contains
p_g = p % g
end if
else if (p % event == EVENT_SCATTER) then
! Then the energy group has been changed by the scattering routine
! meaning gin is now in p % last_g
p_uvw = p % last_uvw
p_g = p % last_g
else
! No scatter, no change in g.
p_uvw = p % coord(p % n_coord) % uvw
p_g = p % g
end if
else
! No actual collision so g has not changed.
p_uvw = p % coord(p % n_coord) % uvw
p_g = p % g
@ -1146,6 +1160,7 @@ contains
! To significantly reduce de-referencing, point matxs to the
! macroscopic Mgxs for the material of interest
matxs => macro_xs(p % material) % obj
! Do same for nucxs, point it to the microscopic nuclide data of interest
if (i_nuclide > 0) then
nucxs => nuclides_MG(i_nuclide) % obj
@ -1172,6 +1187,7 @@ contains
! All events score to a flux bin. We actually use a collision
! estimator in place of an analog one since there is no way to count
! 'events' exactly for the flux
if (survival_biasing) then
! We need to account for the fact that some weight was already
! absorbed
@ -1199,6 +1215,7 @@ contains
else
score = p % last_wgt
end if
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('total', p_g, UVW=p_uvw) / &
@ -1228,20 +1245,21 @@ contains
else
score = p % last_wgt
end if
if (i_nuclide > 0) then
score = score * nucxs % get_xs('inv_vel', p_g, UVW=p_uvw) / &
matxs % get_xs('total', p_g, UVW=p_uvw) * flux
score = score * nucxs % get_xs('inverse_velocity', p_g, UVW=p_uvw) &
/ nucxs % get_xs('total', p_g, UVW=p_uvw) * flux
else
score = matxs % get_xs('inv_vel', p_g, UVW=p_uvw) * flux
score = score * matxs % get_xs('inverse_velocity', p_g, UVW=p_uvw) &
/ matxs % get_xs('total', p_g, UVW=p_uvw) * flux
end if
else
! For inverse velocity, we need no cross section
if (i_nuclide > 0) then
score = score * nucxs % get_xs('inv_vel', p_g, UVW=p_uvw) * &
atom_density * flux
score = flux * nucxs % get_xs('inverse_velocity', p_g, UVW=p_uvw)
else
score = flux * matxs % get_xs('inv_vel', p_g, UVW=p_uvw)
score = flux * matxs % get_xs('inverse_velocity', p_g, UVW=p_uvw)
end if
end if
@ -1358,6 +1376,14 @@ contains
case (SCORE_FISSION)
! make sure the correct energy is used
if (t % estimator == ESTIMATOR_TRACKLENGTH) then
p_g = p % g
else
p_g = p % last_g
end if
if (t % estimator == ESTIMATOR_ANALOG) then
if (survival_biasing) then
! No fission events occur if survival biasing is on -- need to
@ -1386,13 +1412,20 @@ contains
score = nucxs % get_xs('fission', p_g, UVW=p_uvw) * &
atom_density * flux
else
score = flux * material_xs % fission
score = matxs % get_xs('fission', p_g, UVW=p_uvw) * flux
end if
end if
case (SCORE_NU_FISSION)
! make sure the correct energy is used
if (t % estimator == ESTIMATOR_TRACKLENGTH) then
p_g = p % g
else
p_g = p % last_g
end if
if (t % estimator == ESTIMATOR_ANALOG) then
if (survival_biasing .or. p % fission) then
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
@ -1401,8 +1434,7 @@ contains
! neutrons were emitted with different energies, multiple
! outgoing energy bins may have been scored to. The following
! logic treats this special case and results to multiple bins
call score_fission_eout_mg(p, t, score_index, i_nuclide, &
atom_density)
call score_fission_eout(p, t, score_index, score_bin)
cycle SCORE_LOOP
end if
end if
@ -1441,12 +1473,418 @@ contains
score = nucxs % get_xs('nu_fission', p_g, UVW=p_uvw) * &
atom_density * flux
else
score = material_xs % nu_fission * flux
score = matxs % get_xs('nu_fission', p_g, UVW=p_uvw) * flux
end if
end if
case (SCORE_PROMPT_NU_FISSION)
! make sure the correct energy is used
if (t % estimator == ESTIMATOR_TRACKLENGTH) then
p_g = p % g
else
p_g = p % last_g
end if
if (t % estimator == ESTIMATOR_ANALOG) then
if (survival_biasing .or. p % fission) then
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
! Normally, we only need to make contributions to one scoring
! bin. However, in the case of fission, since multiple fission
! neutrons were emitted with different energies, multiple
! outgoing energy bins may have been scored to. The following
! logic treats this special case and results to multiple bins
call score_fission_eout(p, t, score_index, score_bin)
cycle SCORE_LOOP
end if
end if
if (survival_biasing) then
! No fission events occur if survival biasing is on -- need to
! calculate fraction of absorptions that would have resulted in
! nu-fission
score = p % absorb_wgt * flux
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('prompt_nu_fission', p_g, UVW=p_uvw) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
else
score = score * &
matxs % get_xs('prompt_nu_fission', p_g, UVW=p_uvw) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
end if
else
! Skip any non-fission events
if (.not. p % fission) cycle SCORE_LOOP
! If there is no outgoing energy filter, than we only need to
! score to one bin. For the score to be 'analog', we need to
! score the number of particles that were banked in the fission
! bank. Since this was weighted by 1/keff, we multiply by keff
! to get the proper score.
score = keff * p % wgt_bank * (ONE - sum(p % n_delayed_bank) &
/ real(p % n_bank, 8))
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('fission', p_g, UVW=p_uvw) / &
matxs % get_xs('fission', p_g, UVW=p_uvw)
end if
end if
else
if (i_nuclide > 0) then
score = nucxs % get_xs('prompt_nu_fission', p_g, UVW=p_uvw) * &
atom_density * flux
else
score = matxs % get_xs('prompt_nu_fission', p_g, UVW=p_uvw) * flux
end if
end if
case (SCORE_DELAYED_NU_FISSION)
! make sure the correct energy is used
if (t % estimator == ESTIMATOR_TRACKLENGTH) then
p_g = p % g
else
p_g = p % last_g
end if
! Set the delayedgroup filter index and the number of delayed group bins
dg_filter = t % find_filter(FILTER_DELAYEDGROUP)
if (t % estimator == ESTIMATOR_ANALOG) then
if (survival_biasing .or. p % fission) then
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
! Normally, we only need to make contributions to one scoring
! bin. However, in the case of fission, since multiple fission
! neutrons were emitted with different energies, multiple
! outgoing energy bins may have been scored to. The following
! logic treats this special case and results to multiple bins
call score_fission_eout(p, t, score_index, score_bin)
cycle SCORE_LOOP
end if
end if
if (survival_biasing) then
! No fission events occur if survival biasing is on -- need to
! calculate fraction of absorptions that would have resulted in
! nu-fission
if (matxs % get_xs('absorption', p_g, UVW=p_uvw) > ZERO) then
if (dg_filter > 0) then
select type(filt => t % filters(dg_filter) % obj)
type is (DelayedGroupFilter)
! Loop over all delayed group bins and tally to them
! individually
do d_bin = 1, filt % n_bins
! Get the delayed group for this bin
d = filt % groups(d_bin)
score = p % absorb_wgt * flux
if (i_nuclide > 0) then
score = score * nucxs % get_xs('delayed_nu_fission', &
p_g, UVW=p_uvw, dg=d) / &
nucxs % get_xs('absorption', p_g, UVW=p_uvw)
else
score = score * matxs % get_xs('delayed_nu_fission', &
p_g, UVW=p_uvw, dg=d) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
end if
call score_fission_delayed_dg(t, d_bin, score, score_index)
end do
cycle SCORE_LOOP
end select
else
score = p % absorb_wgt * flux
if (i_nuclide > 0) then
score = score * nucxs % get_xs('delayed_nu_fission', p_g, &
UVW=p_uvw) / nucxs % get_xs('absorption', p_g, UVW=p_uvw)
else
score = score * matxs % get_xs('delayed_nu_fission', p_g, &
UVW=p_uvw) / matxs % get_xs('absorption', p_g, UVW=p_uvw)
end if
end if
end if
else
! Skip any non-fission events
if (.not. p % fission) cycle SCORE_LOOP
! If there is no outgoing energy filter, than we only need to
! score to one bin. For the score to be 'analog', we need to
! score the number of particles that were banked in the fission
! bank. Since this was weighted by 1/keff, we multiply by keff
! to get the proper score.
! Check if the delayed group filter is present
if (dg_filter > 0) then
select type(filt => t % filters(dg_filter) % obj)
type is (DelayedGroupFilter)
! Loop over all delayed group bins and tally to them
! individually
do d_bin = 1, filt % n_bins
! Get the delayed group for this bin
d = filt % groups(d_bin)
score = keff * p % wgt_bank / p % n_bank * &
p % n_delayed_bank(d)
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('fission', p_g, UVW=p_uvw) / &
matxs % get_xs('fission', p_g, UVW=p_uvw)
end if
call score_fission_delayed_dg(t, d_bin, score, score_index)
end do
cycle SCORE_LOOP
end select
else
score = keff * p % wgt_bank / p % n_bank * sum(p % n_delayed_bank)
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('fission', p_g, UVW=p_uvw) / &
matxs % get_xs('fission', p_g, UVW=p_uvw)
end if
end if
end if
else
! Check if the delayed group filter is present
if (dg_filter > 0) then
select type(filt => t % filters(dg_filter) % obj)
type is (DelayedGroupFilter)
! Loop over all delayed group bins and tally to them
! individually
do d_bin = 1, filt % n_bins
! Get the delayed group for this bin
d = filt % groups(d_bin)
if (i_nuclide > 0) then
score = nucxs % get_xs('delayed_nu_fission', p_g, &
UVW=p_uvw, dg=d) * atom_density * flux
else
score = matxs % get_xs('delayed_nu_fission', p_g, &
UVW=p_uvw, dg=d) * flux
end if
call score_fission_delayed_dg(t, d_bin, score, score_index)
end do
cycle SCORE_LOOP
end select
else
if (i_nuclide > 0) then
score = nucxs % get_xs('delayed_nu_fission', p_g, UVW=p_uvw) &
* atom_density * flux
else
score = matxs % get_xs('delayed_nu_fission', p_g, UVW=p_uvw) &
* flux
end if
end if
end if
case (SCORE_DECAY_RATE)
! make sure the correct energy is used
if (t % estimator == ESTIMATOR_TRACKLENGTH) then
p_g = p % g
else
p_g = p % last_g
end if
! Set the delayedgroup filter index and the number of delayed group bins
dg_filter = t % find_filter(FILTER_DELAYEDGROUP)
if (t % estimator == ESTIMATOR_ANALOG) then
if (survival_biasing) then
! No fission events occur if survival biasing is on -- need to
! calculate fraction of absorptions that would have resulted in
! nu-fission
if (matxs % get_xs('absorption', p_g, UVW=p_uvw) > ZERO) then
if (dg_filter > 0) then
select type(filt => t % filters(dg_filter) % obj)
type is (DelayedGroupFilter)
! Loop over all delayed group bins and tally to them
! individually
do d_bin = 1, filt % n_bins
! Get the delayed group for this bin
d = filt % groups(d_bin)
score = p % absorb_wgt * flux
if (i_nuclide > 0) then
score = score * nucxs % get_xs('decay_rate', p_g, &
UVW=p_uvw, dg=d) * &
nucxs % get_xs('delayed_nu_fission', p_g, &
UVW=p_uvw, dg=d) / nucxs % get_xs('absorption', &
p_g, UVW=p_uvw)
else
score = score * matxs % get_xs('decay_rate', p_g, &
UVW=p_uvw, dg=d) * &
matxs % get_xs('delayed_nu_fission', p_g, &
UVW=p_uvw, dg=d) / matxs % get_xs('absorption', &
p_g, UVW=p_uvw)
end if
call score_fission_delayed_dg(t, d_bin, score, score_index)
end do
cycle SCORE_LOOP
end select
else
score = ZERO
! If the delayed group filter is not present, compute the score
! by accumulating the absorbed weight times the decay rate times
! the fraction of the delayed-nu-fission xs to the absorption xs
! for all delayed groups.
do d = 1, num_delayed_groups
if (i_nuclide > 0) then
score = score + p % absorb_wgt * &
nucxs % get_xs('decay_rate', p_g, UVW=p_uvw, dg=d) * &
nucxs % get_xs('delayed_nu_fission', p_g, UVW=p_uvw, &
dg=d) / nucxs % get_xs('absorption', p_g, UVW=p_uvw)
else
score = score + p % absorb_wgt * &
matxs % get_xs('decay_rate', p_g, UVW=p_uvw, dg=d) * &
matxs % get_xs('delayed_nu_fission', p_g, UVW=p_uvw, &
dg=d) / matxs % get_xs('absorption', p_g, UVW=p_uvw)
end if
end do
end if
end if
else
! Skip any non-fission events
if (.not. p % fission) cycle SCORE_LOOP
! If there is no outgoing energy filter, than we only need to
! score to one bin. For the score to be 'analog', we need to
! score the number of particles that were banked in the fission
! bank. Since this was weighted by 1/keff, we multiply by keff
! to get the proper score.
score = ZERO
! loop over number of particles banked
do k = 1, p % n_bank
! get the delayed group
g = fission_bank(n_bank - p % n_bank + k) % delayed_group
! Case for tallying delayed emissions
if (g /= 0) then
! determine score based on bank site weight and keff.
if (i_nuclide > 0) then
score = score + keff * atom_density * &
fission_bank(n_bank - p % n_bank + k) % wgt * &
nucxs % get_xs('decay_rate', p_g, UVW=p_uvw, dg=g) * &
nucxs % get_xs('fission', p_g, UVW=p_uvw) / &
matxs % get_xs('fission', p_g, UVW=p_uvw)
else
score = score + keff * &
fission_bank(n_bank - p % n_bank + k) % wgt * &
matxs % get_xs('decay_rate', p_g, UVW=p_uvw, dg=g)
end if
! if the delayed group filter is present, tally to corresponding
! delayed group bin if it exists
if (dg_filter > 0) then
! declare the delayed group filter type
select type(filt => t % filters(dg_filter) % obj)
type is (DelayedGroupFilter)
! loop over delayed group bins until the corresponding bin
! is found
do d_bin = 1, filt % n_bins
d = filt % groups(d_bin)
! check whether the delayed group of the particle is equal
! to the delayed group of this bin
if (d == g) then
call score_fission_delayed_dg(t, d_bin, score, &
score_index)
end if
end do
end select
! Reset the score to zero
score = ZERO
end if
end if
end do
! If the delayed group filter is present, cycle because the
! score_fission_delayed_dg(...) has already tallied the score
if (dg_filter > 0) then
cycle SCORE_LOOP
end if
end if
else
! Check if the delayed group filter is present
if (dg_filter > 0) then
select type(filt => t % filters(dg_filter) % obj)
type is (DelayedGroupFilter)
! Loop over all delayed group bins and tally to them
! individually
do d_bin = 1, filt % n_bins
! Get the delayed group for this bin
d = filt % groups(d_bin)
if (i_nuclide > 0) then
score = nucxs % get_xs('decay_rate', p_g, UVW=p_uvw, dg=d) * &
nucxs % get_xs('delayed_nu_fission', p_g, UVW=p_uvw, &
dg=d) * atom_density * flux
else
score = matxs % get_xs('decay_rate', p_g, UVW=p_uvw, dg=d) * &
matxs % get_xs('delayed_nu_fission', p_g, UVW=p_uvw, &
dg=d) * flux
end if
call score_fission_delayed_dg(t, d_bin, score, score_index)
end do
cycle SCORE_LOOP
end select
else
score = ZERO
! If the delayed group filter is not present, compute the score
! by accumulating the absorbed weight times the decay rate times
! the fraction of the delayed-nu-fission xs to the absorption xs
! for all delayed groups.
do d = 1, num_delayed_groups
if (i_nuclide > 0) then
score = score + atom_density * flux * &
nucxs % get_xs('decay_rate', p_g, UVW=p_uvw, dg=d) * &
nucxs % get_xs('delayed_nu_fission', p_g, UVW=p_uvw, dg=d)
else
score = score + flux * &
matxs % get_xs('decay_rate', p_g, UVW=p_uvw, dg=d) * &
matxs % get_xs('delayed_nu_fission', p_g, UVW=p_uvw, dg=d)
end if
end do
end if
end if
case (SCORE_KAPPA_FISSION)
! make sure the correct energy is used
if (t % estimator == ESTIMATOR_TRACKLENGTH) then
p_g = p % g
else
p_g = p % last_g
end if
if (t % estimator == ESTIMATOR_ANALOG) then
if (survival_biasing) then
! No fission events occur if survival biasing is on -- need to
@ -1472,10 +1910,10 @@ contains
end if
else
if (i_nuclide > 0) then
score = flux * atom_density * &
nucxs % get_xs('kappa_fission', p_g, UVW=p_uvw)
score = nucxs % get_xs('kappa_fission', p_g, UVW=p_uvw) * &
atom_density * flux
else
score = flux * matxs % get_xs('kappa_fission', p_g, UVW=p_uvw)
score = matxs % get_xs('kappa_fission', p_g, UVW=p_uvw) * flux
end if
end if
@ -1931,7 +2369,7 @@ contains
! neutrons produced with different energies.
!===============================================================================
subroutine score_fission_eout_ce(p, t, i_score, score_bin)
subroutine score_fission_eout(p, t, i_score, score_bin)
type(Particle), intent(in) :: p
type(TallyObject), intent(inout) :: t
@ -1950,6 +2388,7 @@ contains
real(8) :: filter_weight ! combined weight of all filters
real(8) :: score ! actual score
real(8) :: E_out ! energy of fission bank site
integer :: g_out ! energy group of fission bank site
! save original outgoing energy bin and score index
i = t % find_filter(FILTER_ENERGYOUT)
@ -1976,14 +2415,31 @@ contains
! determine score based on bank site weight and keff
score = keff * fission_bank(n_bank - p % n_bank + k) % wgt
! determine outgoing energy from fission bank
E_out = fission_bank(n_bank - p % n_bank + k) % E
if (.not. run_CE .and. eo_filt % matches_transport_groups) then
! check if outgoing energy is within specified range on filter
if (E_out < eo_filt % bins(1) .or. E_out > eo_filt % bins(n)) cycle
! determine outgoing energy from fission bank
g_out = int(fission_bank(n_bank - p % n_bank + k) % E)
! change outgoing energy bin
Matching_bins(i) = binary_search(eo_filt % bins, n, E_out)
! change outgoing energy bin
matching_bins(i) = g_out
else
! determine outgoing energy from fission bank
if (run_CE) then
E_out = fission_bank(n_bank - p % n_bank + k) % E
else
E_out = energy_bin_avg(int(fission_bank(n_bank - p % n_bank + k) &
% E))
end if
! check if outgoing energy is within specified range on filter
if (E_out < eo_filt % bins(1) .or. E_out > eo_filt % bins(n)) cycle
! change outgoing energy bin
matching_bins(i) = binary_search(eo_filt % bins, n, E_out)
end if
! Case for tallying prompt neutrons
if (score_bin == SCORE_NU_FISSION .or. &
@ -2046,91 +2502,7 @@ contains
! reset outgoing energy bin and score index
matching_bins(i) = bin_energyout
end subroutine score_fission_eout_ce
subroutine score_fission_eout_mg(p, t, i_score, i_nuclide, atom_density)
type(Particle), intent(in) :: p
type(TallyObject), intent(inout) :: t
integer, intent(in) :: i_score ! index for score
integer, intent(in) :: i_nuclide ! index for nuclide
real(8), intent(in) :: atom_density
integer :: i ! index of outgoing energy filter
integer :: n ! number of energies on filter
integer :: k ! loop index for bank sites
integer :: bin_energyout ! original outgoing energy bin
integer :: i_filter ! index for matching filter bin combination
real(8) :: filter_weight ! combined weight of all filters
real(8) :: score ! actual score
integer :: gout ! energy group of fission bank site
integer :: gin ! energy group of incident particle
real(8) :: E_out
! save original outgoing energy bin and score index
i = t % find_filter(FILTER_ENERGYOUT)
bin_energyout = matching_bins(i)
! Declare the filter type
select type(filt => t % filters(i) % obj)
type is (EnergyoutFilter)
! Get number of energies on filter
n = size(filt % bins)
! Since the creation of fission sites is weighted such that it is
! expected to create n_particles sites, we need to multiply the
! score by keff to get the true nu-fission rate. Otherwise, the sum
! of all nu-fission rates would be ~1.0.
! loop over number of particles banked
do k = 1, p % n_bank
! determine score based on bank site weight and keff
score = keff * fission_bank(n_bank - p % n_bank + k) % wgt
if (i_nuclide > 0) then
if (survival_biasing) then
gin = p % g
else
gin = p % last_g
end if
score = score * atom_density * &
nuclides_MG(i_nuclide) % obj % get_xs('fission', gin, &
UVW=p % last_uvw) / &
macro_xs(p % material) % obj % get_xs('fission', gin, &
UVW=p % last_uvw)
end if
if (filt % matches_transport_groups) then
! determine outgoing energy from fission bank
gout = int(fission_bank(n_bank - p % n_bank + k) % E)
! change outgoing energy bin
matching_bins(i) = gout
else
! determine outgoing energy from fission bank
E_out = energy_bin_avg(int(fission_bank(n_bank - p % n_bank + k) % E))
! check if outgoing energy is within specified range on filter
if (E_out < filt % bins(1) .or. E_out > filt % bins(n)) cycle
! change outgoing energy bin
matching_bins(i) = binary_search(filt % bins, n, E_out)
end if
! determine scoring index and weight for this filter combination
i_filter = sum((matching_bins(1:size(t%filters)) - 1) * t % stride) + 1
filter_weight = product(filter_weights(:size(t % filters)))
! Add score to tally
!$omp atomic
t % results(i_score, i_filter) % value = &
t % results(i_score, i_filter) % value + score * filter_weight
end do
! reset outgoing energy bin and score index
matching_bins(i) = bin_energyout
end select
end subroutine score_fission_eout_mg
end subroutine score_fission_eout
!===============================================================================
! SCORE_FISSION_DELAYED_DG helper function used to increment the tally when a

4
src/tally_filter.F90
View file

@ -890,7 +890,7 @@ contains
! Tallies are ordered in increasing groups, group indices
! however are the opposite, so switch
next_bin = energy_groups - next_bin + 1
next_bin = num_energy_groups - next_bin + 1
else
! Make sure the correct energy is used.
@ -959,7 +959,7 @@ contains
! Tallies are ordered in increasing groups, group indices
! however are the opposite, so switch
next_bin = energy_groups - next_bin + 1
next_bin = num_energy_groups - next_bin + 1
else
! Check if energy of the particle is within energy bins.