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

Merge pull request #739 from samuelshaner/mg-mode-delayed

Browse source 
Treatment of delayed neutrons in multi-group mode
This commit is contained in:
Adam Nelson 2016-11-06 13:35:06 -05:00 • committed by GitHub
commit d9f1bdc761
32 changed files with 3167 additions and 1023 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 eV. The length of this
array should be the number of groups plus 1.
@ -119,8 +121,8 @@ Temperature-dependent data, provided for temperature <TTT>K.
the `nu-fission` data must represent the fission neutron energy
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
same dimensionality as `multiplicity matrix`.
- **inverse-velocity** (*double[]*) -- Average inverse velocity
for each of the groups in the library. This dataset is optional.
**/<library name>/<TTT>K/scatter_data/**

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

@ -1629,8 +1629,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
@ -1811,12 +1809,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. |
+----------------------+---------------------------------------------------+
@ -1845,8 +1841,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 |
@ -1882,8 +1876,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::

22
examples/xml/pincell_multigroup/geometry.xml
View file

@ -1,16 +1,10 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<surface coeffs="0. 0. 0.540" id="1" type="z-cylinder" />
<surface boundary="reflective" coeffs="-0.63" id="20" type="x-plane" />
<surface boundary="reflective" coeffs=" 0.63" id="21" type="x-plane" />
<surface boundary="reflective" coeffs="-0.63" id="22" type="y-plane" />
<surface boundary="reflective" coeffs=" 0.63" id="23" type="y-plane" />
<cell id="1" material="1" region=" -1" />
<cell id="2" material="7" region="1 20 -21 22 -23" />
<cell id="1" material="1" name="cell 1" region="-1" universe="0" />
<cell id="2" material="2" name="cell 2" region="((((1 4) -5) 6) -7)" universe="0" />
<surface coeffs="0 0 0.54" id="1" name="Fuel OR" type="z-cylinder" />
<surface boundary="reflective" coeffs="-0.63" id="4" name="left" type="x-plane" />
<surface boundary="reflective" coeffs="0.63" id="5" name="right" type="x-plane" />
<surface boundary="reflective" coeffs="-0.63" id="6" name="bottom" type="y-plane" />
<surface boundary="reflective" coeffs="0.63" id="7" name="top" type="y-plane" />
</geometry>

62
examples/xml/pincell_multigroup/materials.xml
View file

@ -1,54 +1,12 @@
<?xml version="1.0"?>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<cross_sections>./mgxs.h5</cross_sections>
<!-- UO2 -->
<material id="1">
<density units="macro" value="1.0" />
<macroscopic name="UO2"/>
</material>
<!-- 4.3% MOX -->
<material id="2">
<density units="macro" value="1.0" />
<macroscopic name="MOX1"/>
</material>
<!-- 7.0 MOX -->
<material id="3">
<density units="macro" value="1.0" />
<macroscopic name="MOX2"/>
</material>
<!-- 8.0% MOX -->
<material id="4">
<density units="macro" value="1.0" />
<macroscopic name="MOX3"/>
</material>
<!-- Fission Chamber -->
<material id="5">
<density units="macro" value="1.0" />
<macroscopic name="FC"/>
</material>
<!-- Guide Tube -->
<material id="6">
<density units="macro" value="1.0" />
<macroscopic name="GT"/>
</material>
<!-- Water -->
<material id="7">
<density units="macro" value="1.0" />
<macroscopic name="LWTR"/>
</material>
<!-- Control Rod -->
<material id="8">
<density units="macro" value="1.0" />
<macroscopic name="CR"/>
</material>
<cross_sections>./mgxs.h5</cross_sections>
<material id="1" name="UO2 fuel">
<density units="macro" value="1.0" />
<macroscopic name="UO2" />
</material>
<material id="2" name="Water">
<density units="macro" value="1.0" />
<macroscopic name="LWTR" />
</material>
</materials>

BIN
examples/xml/pincell_multigroup/mgxs.h5
View file

Binary file not shown.

47
examples/xml/pincell_multigroup/settings.xml
View file

@ -1,37 +1,14 @@
<?xml version="1.0"?>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<energy_mode>multi-group</energy_mode>
<!--
Define how many particles to run and for how many batches
in an eigenvalue calculation mode
-->
<eigenvalue>
<batches>100</batches>
<inactive>10</inactive>
<particles>1000</particles>
</eigenvalue>
<!--
Start with uniformally distributed neutron source
with the default energy spectrum of a Maxwellian
and isotropic distribution.
-->
<source>
<space type="box">
<parameters>
-0.63 -0.63 -1E50
0.63 0.63 1E50
</parameters>
</space>
</source>
<output>
<summary>true</summary>
<tallies>true</tallies>
</output>
<survival_biasing>false</survival_biasing>
<eigenvalue>
<particles>1000</particles>
<batches>100</batches>
<inactive>10</inactive>
</eigenvalue>
<source strength="1.0">
<space type="box">
<parameters>-0.63 -0.63 -1 0.63 0.63 1</parameters>
</space>
</source>
<energy_mode>multi-group</energy_mode>
</settings>

4
examples/xml/pincell_multigroup/tallies.xml
View file

@ -1,4 +1,4 @@
<?xml version="1.0"?>
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<mesh id="1" type="regular">
<dimension>100 100 1</dimension>
@ -6,7 +6,7 @@
<upper_right>0.63 0.63 1e+50</upper_right>
</mesh>
<tally id="1" name="tally 1">
<filter bins="1e-11 6.35e-08 1e-05 0.0001 0.001 0.5 1.0 20.0" type="energy" />
<filter bins="1e-05 0.0635 10.0 100.0 1000.0 500000.0 1000000.0 20000000.0" type="energy" />
<filter bins="1" type="mesh" />
<scores>flux fission nu-fission</scores>
</tally>

121
openmc/mgxs/library.py
View file

@ -63,8 +63,8 @@ class Library(object):
The highest legendre moment in the scattering matrices (default is 0)
energy_groups : openmc.mgxs.EnergyGroups
Energy group structure for energy condensation
delayed_groups : list of int
Delayed groups to filter out the xs
num_delayed_groups : int
Number of delayed groups
estimator : str or None
The tally estimator used to compute multi-group cross sections. If None,
the default for each MGXS type is used.
@ -99,7 +99,7 @@ class Library(object):
self._domain_type = None
self._domains = 'all'
self._energy_groups = None
self._delayed_groups = None
self._num_delayed_groups = 0
self._correction = 'P0'
self._legendre_order = 0
self._tally_trigger = None
@ -132,7 +132,7 @@ class Library(object):
clone._correction = self.correction
clone._legendre_order = self.legendre_order
clone._energy_groups = copy.deepcopy(self.energy_groups, memo)
clone._delayed_groups = copy.deepcopy(self.delayed_groups, memo)
clone._num_delayed_groups = self.num_delayed_groups
clone._tally_trigger = copy.deepcopy(self.tally_trigger, memo)
clone._all_mgxs = copy.deepcopy(self.all_mgxs)
clone._sp_filename = self._sp_filename
@ -202,8 +202,8 @@ class Library(object):
return self._energy_groups
@property
def delayed_groups(self):
return self._delayed_groups
def num_delayed_groups(self):
return self._num_delayed_groups
@property
def correction(self):
@ -225,13 +225,6 @@ class Library(object):
def num_groups(self):
return self.energy_groups.num_groups
@property
def num_delayed_groups(self):
if self.delayed_groups == None:
return 0
else:
return len(self.delayed_groups)
@property
def all_mgxs(self):
return self._all_mgxs
@ -331,22 +324,14 @@ class Library(object):
cv.check_type('energy groups', energy_groups, openmc.mgxs.EnergyGroups)
self._energy_groups = energy_groups
@delayed_groups.setter
def delayed_groups(self, delayed_groups):
@num_delayed_groups.setter
def num_delayed_groups(self, num_delayed_groups):
if delayed_groups != None:
cv.check_type('delayed groups', delayed_groups, list, int)
cv.check_greater_than('num delayed groups', len(delayed_groups), 0)
# Check that the groups are within [1, MAX_DELAYED_GROUPS]
for group in delayed_groups:
cv.check_greater_than('delayed group', group, 0)
cv.check_less_than('delayed group', group,
openmc.mgxs.MAX_DELAYED_GROUPS,
equality=True)
self._delayed_groups = delayed_groups
cv.check_less_than('num delayed groups', num_delayed_groups,
openmc.mgxs.MAX_DELAYED_GROUPS, equality=True)
cv.check_greater_than('num delayed groups', num_delayed_groups, 0,
equality=True)
self._num_delayed_groups = num_delayed_groups
@correction.setter
def correction(self, correction):
@ -431,7 +416,12 @@ class Library(object):
mgxs.estimator = self.estimator
if mgxs_type in openmc.mgxs.MDGXS_TYPES:
mgxs.delayed_groups = self.delayed_groups
if self.num_delayed_groups == 0:
mgxs.delayed_groups = None
else:
delayed_groups \
= list(range(1,self.num_delayed_groups+1))
mgxs.delayed_groups = delayed_groups
# If a tally trigger was specified, add it to the MGXS
if self.tally_trigger is not None:
@ -896,6 +886,7 @@ class Library(object):
if nuclide != 'total':
name += '_' + nuclide
xsdata = openmc.XSdata(name, self.energy_groups)
xsdata.num_delayed_groups = self.num_delayed_groups
if order is None:
# Set the order to the Library's order (the defualt behavior)
@ -922,41 +913,79 @@ class Library(object):
mymgxs = self.get_mgxs(domain, 'nu-transport')
xsdata.set_total_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide],
subdomain=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')
@ -964,6 +993,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))):
@ -974,6 +1004,7 @@ class Library(object):
nuclide=[nuclide],
subdomain=subdomain)
using_multiplicity = True
else:
using_multiplicity = False
@ -995,8 +1026,8 @@ class Library(object):
# absorption cross section.
if 'total' in self.mgxs_types or 'transport' 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 +1082,9 @@ class Library(object):
xs_type = 'macro'
# Initialize file
mgxs_file = openmc.MGXSLibrary(self.energy_groups)
mgxs_file = openmc.MGXSLibrary(self.energy_groups,
num_delayed_groups=\
self.num_delayed_groups)
if self.domain_type == 'mesh':
# Create the xsdata objects and add to the mgxs_file

7
openmc/mgxs/mdgxs.py
View file

@ -184,7 +184,7 @@ class MDGXS(MGXS):
cv.check_less_than('delayed group', group, MAX_DELAYED_GROUPS,
equality=True)
self._delayed_groups = delayed_groups
self._delayed_groups = delayed_groups
@property
def filters(self):
@ -1464,6 +1464,11 @@ class Beta(MDGXS):
\beta_{d,g} &= \frac{\langle \nu^d \sigma_f \phi \rangle}
{\langle \nu \sigma_f \phi \rangle}
NOTE: The Beta MGXS is the delayed neutron fraction computed directly from
the nuclear data. Often the delayed neutron fraction is
"importance-weighted" by the adjoint flux and called "beta-effective". It
is important to make clear that this Beta is not importance-weighted.
Parameters
----------
domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh

952
openmc/mgxs_library.py
View file

File diff suppressed because it is too large Load diff

10
scripts/openmc-update-mgxs
View file

@ -85,13 +85,13 @@ if __name__ == '__main__':
temp = np.array(temp.split())
group_structure = temp.astype(np.float)
energy_groups = openmc.mgxs.EnergyGroups(group_structure)
temp = tree.find('inverse_velocities')
temp = tree.find('inverse-velocity')
if temp is not None:
temp = temp.text.strip()
temp = np.array(temp.split())
inverse_velocities = temp.astype(np.float)
inverse_velocity = temp.astype(np.float)
else:
inverse_velocities = None
inverse_velocity = None
xsd = []
names = []
@ -166,8 +166,8 @@ if __name__ == '__main__':
total.shape = xsd[i].vector_shape
xsd[i].set_total(total, temperature)
if inverse_velocities is not None:
xsd[i].set_inverse_velocities(inverse_velocities, temperature)
if inverse_velocity is not None:
xsd[i].set_inverse_velocity(inverse_velocity, temperature)
temp = get_data(xsdata_elem, 'absorption')
temp = np.array(temp.split())

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

@ -112,7 +112,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

@ -98,6 +98,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
@ -2586,6 +2587,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

82
src/input_xml.F90
View file

@ -1062,19 +1062,6 @@ contains
end select
end if
! Check whether create fission sites
if (run_mode == MODE_FIXEDSOURCE) then
if (check_for_node(doc, "create_fission_neutrons")) then
call get_node_value(doc, "create_fission_neutrons", temp_str)
temp_str = to_lower(temp_str)
if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') then
create_fission_neutrons = .true.
else if (trim(temp_str) == 'false' .or. trim(temp_str) == '0') then
create_fission_neutrons = .false.
end if
end if
end if
! Check for tabular_legendre options
if (check_for_node(doc, "tabular_legendre")) then
@ -1100,6 +1087,19 @@ contains
end if
end if
! Check whether create fission sites
if (run_mode == MODE_FIXEDSOURCE) then
if (check_for_node(doc, "create_fission_neutrons")) then
call get_node_value(doc, "create_fission_neutrons", temp_str)
temp_str = to_lower(temp_str)
if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') then
create_fission_neutrons = .true.
else if (trim(temp_str) == 'false' .or. trim(temp_str) == '0') then
create_fission_neutrons = .false.
end if
end if
end if
! Close settings XML file
call close_xmldoc(doc)
@ -3017,8 +3017,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
@ -3043,8 +3043,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
@ -3058,14 +3058,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)
@ -3556,7 +3548,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
@ -3569,12 +3567,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')
@ -4611,15 +4603,23 @@ contains
! Open file for reading
file_id = file_open(path_cross_sections, 'r', parallel=.true.)
if (attribute_exists(file_id, "groups")) then
! Get neutron group count
call read_attribute(energy_groups, file_id, "groups")
if (attribute_exists(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 (attribute_exists(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
end if
allocate(rev_energy_bins(num_energy_groups + 1))
allocate(energy_bins(num_energy_groups + 1))
if (attribute_exists(file_id, "group structure")) then
! Get neutron group structure
call read_attribute(energy_bins, file_id, "group structure")
@ -4628,10 +4628,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

84
src/mgxs_data.F90
View file

@ -22,26 +22,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!")
@ -62,23 +61,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
@ -105,7 +87,9 @@ contains
! First find out the data representation
if (attribute_exists(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
@ -120,20 +104,24 @@ 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
call already_read % add(name)
end if
end do NUCLIDE_LOOP
end do MATERIAL_LOOP
@ -149,6 +137,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.
@ -163,12 +152,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))
@ -176,27 +165,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
!===============================================================================
@ -205,16 +199,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 eV
real(8) :: kT ! temperature in eV
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
@ -238,5 +233,4 @@ contains
end subroutine get_mat_kTs
end module mgxs_data

2070
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

606
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
@ -1104,6 +1108,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
@ -1115,7 +1124,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
@ -1126,16 +1137,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
@ -1144,6 +1158,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) &
/ matxs % get_xs('absorption', 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('absorption', 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
@ -1413,11 +1445,11 @@ contains
score = p % absorb_wgt * flux
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('nu_fission', p_g, UVW=p_uvw) / &
nucxs % get_xs('nu-fission', p_g, UVW=p_uvw) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
else
score = score * &
matxs % get_xs('nu_fission', p_g, UVW=p_uvw) / &
matxs % get_xs('nu-fission', p_g, UVW=p_uvw) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
end if
else
@ -1438,15 +1470,421 @@ contains
else
if (i_nuclide > 0) then
score = nucxs % get_xs('nu_fission', p_g, UVW=p_uvw) * &
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) / &
matxs % 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) / matxs % 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) / matxs % 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) / matxs % 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
@ -1463,19 +1901,19 @@ contains
end if
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('kappa_fission', p_g, UVW=p_uvw) / &
nucxs % get_xs('kappa-fission', p_g, UVW=p_uvw) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
else
score = score * &
matxs % get_xs('kappa_fission', p_g, UVW=p_uvw) / &
matxs % get_xs('kappa-fission', p_g, UVW=p_uvw) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
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
@ -1937,7 +2375,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
@ -1956,6 +2394,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)
@ -1982,14 +2421,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. &
@ -2052,91 +2508,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.

BIN
tests/1d_mgxs.h5
View file

Binary file not shown.

2
tests/test_mg_basic/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
1.160653E+00 6.227102E-02
1.086852E+00 2.677280E-02

2
tests/test_mg_max_order/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
1.098323E+00 1.979833E-02
1.077247E+00 2.412834E-02

2
tests/test_mg_nuclide/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
1.160653E+00 6.227102E-02
1.086852E+00 2.677280E-02

2
tests/test_mg_tallies/results_true.dat
View file

@ -1 +1 @@
6ab4a178b642f886eea9fddbdf67a8d6d3e0d1968f06ee52d7fd5ecdc1b4a796398e9b3a8c15ac72fb1a011d13a5fcc59a868781d38fb11787e060bb2e602a7b
3d2ce1b8bdd558fe9f8560e8bb91455e1fedd80a47349344399e22f5b0fac1391f1721e7f12a42bb0c458fd6c87ebf9fa38fe9539fcf69292b21cebf8e5f8989

2
tests/test_mgxs_library_ce_to_mg/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
1.140724E+00 2.952379E-02
1.133362E+00 1.746521E-02

5
tests/test_mgxs_library_condense/test_mgxs_library_condense.py
View file

@ -22,9 +22,6 @@ class MGXSTestHarness(PyAPITestHarness):
# Initialize a two-group structure
energy_groups = openmc.mgxs.EnergyGroups(group_edges=[0, 0.625, 20.e6])
# Initialize a six-delayed-group structure
delayed_groups = list(range(1,7))
# Initialize MGXS Library for a few cross section types
self.mgxs_lib = openmc.mgxs.Library(self._input_set.geometry)
self.mgxs_lib.by_nuclide = False
@ -33,7 +30,7 @@ class MGXSTestHarness(PyAPITestHarness):
self.mgxs_lib.mgxs_types = openmc.mgxs.MGXS_TYPES + \
openmc.mgxs.MDGXS_TYPES
self.mgxs_lib.energy_groups = energy_groups
self.mgxs_lib.delayed_groups = delayed_groups
self.mgxs_lib.num_delayed_groups = 6
self.mgxs_lib.legendre_order = 3
self.mgxs_lib.domain_type = 'material'
self.mgxs_lib.build_library()

5
tests/test_mgxs_library_distribcell/test_mgxs_library_distribcell.py
View file

@ -22,9 +22,6 @@ class MGXSTestHarness(PyAPITestHarness):
# Initialize a one-group structure
energy_groups = openmc.mgxs.EnergyGroups(group_edges=[0, 20.e6])
# Initialize a six-delayed-group structure
delayed_groups = list(range(1,7))
# Initialize MGXS Library for a few cross section types
# for one material-filled cell in the geometry
self.mgxs_lib = openmc.mgxs.Library(self._input_set.geometry)
@ -34,7 +31,7 @@ class MGXSTestHarness(PyAPITestHarness):
self.mgxs_lib.mgxs_types = openmc.mgxs.MGXS_TYPES + \
openmc.mgxs.MDGXS_TYPES
self.mgxs_lib.energy_groups = energy_groups
self.mgxs_lib.delayed_groups = delayed_groups
self.mgxs_lib.num_delayed_groups = 6
self.mgxs_lib.legendre_order = 3
self.mgxs_lib.domain_type = 'distribcell'
cells = self.mgxs_lib.openmc_geometry.get_all_material_cells()

5
tests/test_mgxs_library_hdf5/test_mgxs_library_hdf5.py
View file

@ -29,9 +29,6 @@ class MGXSTestHarness(PyAPITestHarness):
# Initialize a two-group structure
energy_groups = openmc.mgxs.EnergyGroups(group_edges=[0, 0.625, 20.e6])
# Initialize a six-delayed-group structure
delayed_groups = list(range(1,7))
# Initialize MGXS Library for a few cross section types
self.mgxs_lib = openmc.mgxs.Library(self._input_set.geometry)
self.mgxs_lib.by_nuclide = False
@ -40,7 +37,7 @@ class MGXSTestHarness(PyAPITestHarness):
self.mgxs_lib.mgxs_types = openmc.mgxs.MGXS_TYPES + \
openmc.mgxs.MDGXS_TYPES
self.mgxs_lib.energy_groups = energy_groups
self.mgxs_lib.delayed_groups = delayed_groups
self.mgxs_lib.num_delayed_groups = 6
self.mgxs_lib.legendre_order = 3
self.mgxs_lib.domain_type = 'material'
self.mgxs_lib.build_library()

5
tests/test_mgxs_library_mesh/test_mgxs_library_mesh.py
View file

@ -18,9 +18,6 @@ class MGXSTestHarness(PyAPITestHarness):
# Initialize a one-group structure
energy_groups = openmc.mgxs.EnergyGroups(group_edges=[0, 20.e6])
# Initialize a six-delayed-group structure
delayed_groups = list(range(1,7))
# Initialize MGXS Library for a few cross section types
# for one material-filled cell in the geometry
self.mgxs_lib = openmc.mgxs.Library(self._input_set.geometry)
@ -30,7 +27,7 @@ class MGXSTestHarness(PyAPITestHarness):
self.mgxs_lib.mgxs_types = openmc.mgxs.MGXS_TYPES + \
openmc.mgxs.MDGXS_TYPES
self.mgxs_lib.energy_groups = energy_groups
self.mgxs_lib.delayed_groups = delayed_groups
self.mgxs_lib.num_delayed_groups = 6
self.mgxs_lib.legendre_order = 3
self.mgxs_lib.domain_type = 'mesh'

5
tests/test_mgxs_library_no_nuclides/test_mgxs_library_no_nuclides.py
View file

@ -22,9 +22,6 @@ class MGXSTestHarness(PyAPITestHarness):
# Initialize a two-group structure
energy_groups = openmc.mgxs.EnergyGroups(group_edges=[0, 0.625, 20.e6])
# Initialize a six-delayed-group structure
delayed_groups = list(range(1,7))
# Initialize MGXS Library for a few cross section types
self.mgxs_lib = openmc.mgxs.Library(self._input_set.geometry)
self.mgxs_lib.by_nuclide = False
@ -33,7 +30,7 @@ class MGXSTestHarness(PyAPITestHarness):
self.mgxs_lib.mgxs_types = openmc.mgxs.MGXS_TYPES + \
openmc.mgxs.MDGXS_TYPES
self.mgxs_lib.energy_groups = energy_groups
self.mgxs_lib.delayed_groups = delayed_groups
self.mgxs_lib.num_delayed_groups = 6
self.mgxs_lib.legendre_order = 3
self.mgxs_lib.domain_type = 'material'
self.mgxs_lib.build_library()