Merge branch 'develop' into tally-mean

This commit is contained in:
wbinventor@gmail.com 2016-01-23 20:08:21 -05:00
commit 235fdf4874
207 changed files with 2583 additions and 7947 deletions

3
.gitignore vendored
View file

@ -43,6 +43,8 @@ results_test.dat
# Test build files
tests/build/
tests/coverage/
tests/memcheck/
tests/ctestscript.run
# HDF5 files
@ -60,6 +62,7 @@ data/nndc
#Images
*.ppm
*.voxel
# PyCharm project configuration files
.idea

View file

@ -122,7 +122,7 @@ if(CMAKE_Fortran_COMPILER_ID STREQUAL GNU)
list(APPEND ldflags -pg)
endif()
if(optimize)
list(APPEND f90flags -O3)
list(APPEND f90flags -O3 -flto -fuse-linker-plugin)
endif()
if(openmp)
list(APPEND f90flags -fopenmp)
@ -263,7 +263,7 @@ endif()
# set compile flags. Note that this sets the COMPILE_OPTIONS property (also
# available only in 2.8.12+) rather than the COMPILE_FLAGS property, which is
# deprecated. The former can handle lists whereas the latter cannot.
if(CMAKE_VERSION VERSION_LESS 4.8.12)
if(CMAKE_VERSION VERSION_LESS 2.8.12)
string(REPLACE ";" " " f90flags "${f90flags}")
set_property(TARGET ${program} PROPERTY COMPILE_FLAGS "${f90flags}")
else()
@ -313,113 +313,36 @@ include(CTest)
# Get a list of all the tests to run
file(GLOB_RECURSE TESTS ${CMAKE_CURRENT_SOURCE_DIR}/tests/test_*.py)
# Check for MEM_CHECK and COVERAGE variables
if (DEFINED ENV{MEM_CHECK})
set(MEM_CHECK $ENV{MEM_CHECK})
else(DEFINED ENV{MEM_CHECK})
set(MEM_CHECK FALSE)
endif(DEFINED ENV{MEM_CHECK})
if (DEFINED ENV{COVERAGE})
set(COVERAGE $ENV{COVERAGE})
else(DEFINED ENV{COVERAGE})
set(COVERAGE FALSE)
endif(DEFINED ENV{COVERAGE})
# Loop through all the tests
foreach(test ${TESTS})
# Get test information
get_filename_component(TEST_NAME ${test} NAME)
get_filename_component(TEST_PATH ${test} PATH)
# Check for running standard tests (no valgrind, no gcov)
if(NOT ${MEM_CHECK} AND NOT ${COVERAGE})
if (DEFINED ENV{MEM_CHECK})
# Generate input files if needed
if (NOT EXISTS "${TEST_PATH}/geometry.xml")
execute_process(COMMAND ${PYTHON_EXECUTABLE} ${TEST_NAME} --build-inputs
WORKING_DIRECTORY ${TEST_PATH})
endif()
# Add serial test
add_test(NAME ${TEST_NAME}
WORKING_DIRECTORY ${TEST_PATH}
COMMAND $<TARGET_FILE:openmc>)
else()
# Check serial/parallel
if (${MPI_ENABLED})
# Preform a parallel test
add_test(NAME ${TEST_NAME}
WORKING_DIRECTORY ${TEST_PATH}
COMMAND ${PYTHON_EXECUTABLE} ${TEST_NAME} --exe $<TARGET_FILE:openmc>
--mpi_exec $ENV{MPI_DIR}/bin/mpiexec)
else(${MPI_ENABLED})
else()
# Perform a serial test
add_test(NAME ${TEST_NAME}
WORKING_DIRECTORY ${TEST_PATH}
COMMAND ${PYTHON_EXECUTABLE} ${TEST_NAME} --exe $<TARGET_FILE:openmc>)
endif(${MPI_ENABLED})
# Handle special case for valgrind and gcov (run openmc directly, no python)
else(NOT ${MEM_CHECK} AND NOT ${COVERAGE})
# If a plot test is encountered, run with "-p"
if (${test} MATCHES "test_plot")
# Perform serial valgrind and coverage test with plot flag
add_test(NAME ${TEST_NAME}
WORKING_DIRECTORY ${TEST_PATH}
COMMAND $<TARGET_FILE:openmc> -p ${TEST_PATH})
elseif(${test} MATCHES "test_filter_distribcell")
# Add each case for distribcell tests
add_test(NAME ${TEST_NAME}_case-1
WORKING_DIRECTORY ${TEST_PATH}/case-1
COMMAND $<TARGET_FILE:openmc> ${TEST_PATH}/case-1)
add_test(NAME ${TEST_NAME}_case-2
WORKING_DIRECTORY ${TEST_PATH}/case-2
COMMAND $<TARGET_FILE:openmc> ${TEST_PATH}/case-2)
add_test(NAME ${TEST_NAME}_case-3
WORKING_DIRECTORY ${TEST_PATH}/case-3
COMMAND $<TARGET_FILE:openmc> ${TEST_PATH}/case-3)
add_test(NAME ${TEST_NAME}_case-4
WORKING_DIRECTORY ${TEST_PATH}/case-4
COMMAND $<TARGET_FILE:openmc> ${TEST_PATH}/case-4)
# If a restart test is encounted, need to run with -r and restart file(s)
elseif(${test} MATCHES "restart")
# Handle restart tests separately
if(${test} MATCHES "test_statepoint_restart")
set(RESTART_FILE statepoint.07.h5)
elseif(${test} MATCHES "test_sourcepoint_restart")
set(RESTART_FILE statepoint.07.h5 source.07.h5)
elseif(${test} MATCHES "test_particle_restart_eigval")
set(RESTART_FILE particle_9_555.h5)
elseif(${test} MATCHES "test_particle_restart_fixed")
set(RESTART_FILE particle_7_928.h5)
else(${test} MATCHES "test_statepoint_restart")
message(FATAL_ERROR "Restart test ${test} not recognized")
endif(${test} MATCHES "test_statepoint_restart")
# Perform serial valgrind and coverage test
add_test(NAME ${TEST_NAME}
WORKING_DIRECTORY ${TEST_PATH}
COMMAND $<TARGET_FILE:openmc> ${TEST_PATH})
# Perform serial valgrind and coverage restart test
add_test(NAME ${TEST_NAME}_restart
WORKING_DIRECTORY ${TEST_PATH}
COMMAND $<TARGET_FILE:openmc> -r ${RESTART_FILE} ${TEST_PATH})
# Set test dependency
set_tests_properties(${TEST_NAME}_restart PROPERTIES DEPENDS ${TEST_NAME})
# Handle standard tests for valgrind and gcov
else(${test} MATCHES "test_plot")
# Perform serial valgrind and coverage test
add_test(NAME ${TEST_NAME}
WORKING_DIRECTORY ${TEST_PATH}
COMMAND $<TARGET_FILE:openmc> ${TEST_PATH})
endif(${test} MATCHES "test_plot")
endif(NOT ${MEM_CHECK} AND NOT ${COVERAGE})
endif()
endif()
endforeach(test)

View file

@ -28,7 +28,7 @@ Benchmarking
- Khurrum S. Chaudri and Sikander M. Mirza, "Burnup dependent Monte Carlo
neutron physics calculations of IAEA MTR benchmark," *Prog. Nucl. Energy*,
**81**, 43-52 (2015). `<http://dx.doi.org/j.pnucene.2014.12.018>`_
**81**, 43-52 (2015). `<http://dx.doi.org/10.1016/j.pnucene.2014.12.018>`_
- Daniel J. Kelly, Brian N. Aviles, Paul K. Romano, Bryan R. Herman,
Nicholas E. Horelik, and Benoit Forget, "Analysis of select BEAVRS PWR

File diff suppressed because one or more lines are too long

View file

@ -1452,7 +1452,6 @@
],
"source": [
"# Generate tracks for OpenMOC\n",
"openmoc_geometry.initializeFlatSourceRegions()\n",
"track_generator = openmoc.TrackGenerator(openmoc_geometry, num_azim=32, spacing=0.1)\n",
"track_generator.generateTracks()\n",
"\n",
@ -1638,7 +1637,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.10"
"version": "2.7.6"
}
},
"nbformat": 4,

View file

@ -12,24 +12,55 @@ OpenMC, see :ref:`usersguide_install` in the User's Manual.
Installing on Ubuntu through PPA
--------------------------------
For users with Ubuntu 11.10 or later, a binary package for OpenMC is available
through a `Personal Package Archive`_ (PPA) and can be installed through the `APT
package manager`_. Simply enter the following commands into the terminal:
For users with Ubuntu 15.04 or later, a binary package for OpenMC is available
through a `Personal Package Archive`_ (PPA) and can be installed through the
`APT package manager`_. First, add the following PPA to the repository sources:
.. code-block:: sh
sudo apt-add-repository ppa:paulromano/staging
Next, resynchronize the package index files:
.. code-block:: sh
sudo apt-get update
Now OpenMC should be recognized within the repository and can be installed:
.. code-block:: sh
sudo apt-get install openmc
Currently, the binary package does not allow for parallel simulations or use of
HDF5_. Users who need such capabilities should build OpenMC from source as is
described in :ref:`usersguide_install`.
Binary packages from this PPA may exist for earlier versions of Ubuntu, but they
are no longer supported.
.. _Personal Package Archive: https://launchpad.net/~paulromano/+archive/staging
.. _APT package manager: https://help.ubuntu.com/community/AptGet/Howto
.. _HDF5: http://www.hdfgroup.org/HDF5/
---------------------------------------
Installing from Source on Ubuntu 15.04+
---------------------------------------
To build OpenMC from source, several :ref:`prerequisites <prerequisites>` are
needed. If you are using Ubuntu 15.04 or higher, all prerequisites can be
installed directly from the package manager.
.. code-block:: sh
sudo apt-get install gfortran
sudo apt-get install cmake
sudo apt-get install libhdf5-dev
After the packages have been installed, follow the instructions below for
building and installing OpenMC from source.
.. note:: Before Ubuntu 15.04, the HDF5 package included in the Ubuntu Package
archive was not built with support for the Fortran 2003 HDF5
interface, which is needed by OpenMC. If you are using Ubuntu 14.10 or
before you will need to build HDF5 from source.
-------------------------------------------
Installing from Source on Linux or Mac OS X
-------------------------------------------
@ -42,7 +73,6 @@ entering the following commands in a terminal:
git clone https://github.com/mit-crpg/openmc.git
cd openmc
git checkout -b master origin/master
mkdir build && cd build
cmake ..
make

View file

@ -9,8 +9,8 @@ Installing on Ubuntu with PPA
-----------------------------
For users with Ubuntu 15.04 or later, a binary package for OpenMC is available
through a Personal Package Archive (PPA) and can be installed through the APT
package manager. First, add the following PPA to the repository sources:
through a `Personal Package Archive`_ (PPA) and can be installed through the
`APT package manager`_. First, add the following PPA to the repository sources:
.. code-block:: sh
@ -31,10 +31,15 @@ Now OpenMC should be recognized within the repository and can be installed:
Binary packages from this PPA may exist for earlier versions of Ubuntu, but they
are no longer supported.
.. _Personal Package Archive: https://launchpad.net/~paulromano/+archive/staging
.. _APT package manager: https://help.ubuntu.com/community/AptGet/Howto
--------------------
Building from Source
--------------------
.. _prerequisites:
Prerequisites
-------------

View file

@ -237,7 +237,7 @@ class AggregateFilter(object):
self.aggregate_op = aggregate_op
def __hash__(self):
return hash((self.type, self.bins, self.aggregate_op))
return hash(repr(self))
def __eq__(self, other):
return str(other) == str(self)
@ -405,5 +405,5 @@ class AggregateFilter(object):
aggregate_bin_array = np.repeat(aggregate_bin_array, datasize)
# Construct Pandas DataFrame for the AggregateFilter
df = pd.DataFrame({self.aggregate_filter.type: aggregate_bin_array})
return df
df = pd.DataFrame({self.type: aggregate_bin_array})
return df

View file

@ -843,50 +843,20 @@ class MGXS(object):
elif self.domain_type == 'distribcell':
subdomains = np.arange(self.num_subdomains)
else:
subdomains = [0]
subdomains = None
# Clone this MGXS to initialize the subdomain-averaged version
avg_xs = copy.deepcopy(self)
avg_xs._rxn_rate_tally = None
avg_xs._xs_tally = None
avg_xs._sparse = False
# If domain is distribcell, make the new domain 'cell'
if self.domain_type == 'distribcell':
avg_xs.domain_type = 'cell'
# Average each of the tallies across subdomains
for tally_type, tally in avg_xs.tallies.items():
tally_avg = tally.summation(filter_type=self.domain_type,
filter_bins=subdomains)
avg_xs.tallies[tally_type] = tally_avg
# Make condensed tally derived and null out sum, sum_sq
tally._derived = True
tally._sum = None
tally._sum_sq = None
# Get tally data arrays reshaped with one dimension per filter
mean = tally.get_reshaped_data(value='mean')
std_dev = tally.get_reshaped_data(value='std_dev')
# Get the mean, std. dev. across requested subdomains
mean = np.sum(mean[subdomains, ...], axis=0)
std_dev = np.sum(std_dev[subdomains, ...]**2, axis=0)
std_dev = np.sqrt(std_dev)
# If domain is distribcell, make subdomain-averaged a 'cell' domain
domain_filter = tally.find_filter(self._domain_type)
if domain_filter.type == 'distribcell':
domain_filter.type = 'cell'
domain_filter.num_bins = 1
# Reshape averaged data arrays with one dimension for all filters
mean = np.reshape(mean, tally.shape)
std_dev = np.reshape(std_dev, tally.shape)
# Override tally's data with the new condensed data
tally._mean = mean
tally._std_dev = std_dev
# Compute the subdomain-averaged multi-group cross section
avg_xs._domain_type = 'sum({0})'.format(self.domain_type)
avg_xs.sparse = self.sparse
return avg_xs
@ -1249,7 +1219,7 @@ class MGXS(object):
df = self.xs_tally.get_pandas_dataframe(summary=summary)
# Remove the score column since it is homogeneous and redundant
if summary and self.domain_type == 'distribcell':
if summary and 'distribcell' in self.domain_type:
df = df.drop('score', level=0, axis=1)
else:
df = df.drop('score', axis=1)

View file

@ -1003,16 +1003,17 @@ class Tally(object):
A list of filter type strings
(e.g., ['mesh', 'energy']; default is [])
filter_bins : list of Iterables
A list of the filter bins corresponding to the filter_types
parameter (e.g., [(1,), (0., 0.625e-6)]; default is []). Each bin
in the list is the integer ID for 'material', 'surface', 'cell',
'cellborn', and 'universe' Filters. Each bin is an integer for the
cell instance ID for 'distribcell' Filters. Each bin is a 2-tuple of
floats for 'energy' and 'energyout' filters corresponding to the
energy boundaries of the bin of interest. The bin is a (x,y,z)
3-tuple for 'mesh' filters corresponding to the mesh cell of
interest. The order of the bins in the list must correspond to the
filter_types parameter.
A list of tuples of filter bins corresponding to the filter_types
parameter (e.g., [(1,), ((0., 0.625e-6),)]; default is []). Each
tuple contains bins for the corresponding filter type in the filters
parameter. Each bins is the integer ID for 'material', 'surface',
'cell', 'cellborn', and 'universe' Filters. Each bin is an integer
for the cell instance ID for 'distribcell' Filters. Each bin is a
2-tuple of floats for 'energy' and 'energyout' filters corresponding
to the energy boundaries of the bin of interest. The bin is an
(x,y,z) 3-tuple for 'mesh' filters corresponding to the mesh cell
of interest. The order of the bins in the list must correspond to
the filter_types parameter.
Returns
-------
@ -1174,16 +1175,17 @@ class Tally(object):
A list of filter type strings
(e.g., ['mesh', 'energy']; default is [])
filter_bins : list of Iterables
A list of the filter bins corresponding to the filter_types
parameter (e.g., [(1,), (0., 0.625e-6)]; default is []). Each bin
in the list is the integer ID for 'material', 'surface', 'cell',
'cellborn', and 'universe' Filters. Each bin is an integer for the
cell instance ID for 'distribcell' Filters. Each bin is a 2-tuple of
floats for 'energy' and 'energyout' filters corresponding to the
energy boundaries of the bin of interest. The bin is a (x,y,z)
3-tuple for 'mesh' filters corresponding to the mesh cell of
interest. The order of the bins in the list must correspond to the
filter_types parameter.
A list of tuples of filter bins corresponding to the filter_types
parameter (e.g., [(1,), ((0., 0.625e-6),)]; default is []). Each
tuple contains bins for the corresponding filter type in the filters
parameter. Each bins is the integer ID for 'material', 'surface',
'cell', 'cellborn', and 'universe' Filters. Each bin is an integer
for the cell instance ID for 'distribcell' Filters. Each bin is a
2-tuple of floats for 'energy' and 'energyout' filters corresponding
to the energy boundaries of the bin of interest. The bin is an
(x,y,z) 3-tuple for 'mesh' filters corresponding to the mesh cell
of interest. The order of the bins in the list must correspond to
the filter_types parameter.
nuclides : list of str
A list of nuclide name strings
(e.g., ['U-235', 'U-238']; default is [])
@ -2641,16 +2643,17 @@ class Tally(object):
A list of filter type strings
(e.g., ['mesh', 'energy']; default is [])
filter_bins : list of Iterables
A list of the filter bins corresponding to the filter_types
parameter (e.g., [(1,), (0., 0.625e-6)]; default is []). Each bin
in the list is the integer ID for 'material', 'surface', 'cell',
'cellborn', and 'universe' Filters. Each bin is an integer for the
cell instance ID for 'distribcell' Filters. Each bin is a 2-tuple of
floats for 'energy' and 'energyout' filters corresponding to the
energy boundaries of the bin of interest. The bin is a (x,y,z)
3-tuple for 'mesh' filters corresponding to the mesh cell of
interest. The order of the bins in the list must correspond to the
filter_types parameter.
A list of tuples of filter bins corresponding to the filter_types
parameter (e.g., [(1,), ((0., 0.625e-6),)]; default is []). Each
tuple contains bins to slice for the corresponding filter type in
the filters parameter. Each bins is the integer ID for 'material',
'surface', 'cell', 'cellborn', and 'universe' Filters. Each bin is
an integer for the cell instance ID for 'distribcell' Filters. Each
bin is a 2-tuple of floats for 'energy' and 'energyout' filters
corresponding to the energy boundaries of the bin of interest. The
bin is an (x,y,z) 3-tuple for 'mesh' filters corresponding to the
mesh cell of interest. The order of the bins in the list must
correspond to the filter_types parameter.
nuclides : list of str
A list of nuclide name strings
(e.g., ['U-235', 'U-238']; default is [])
@ -2736,6 +2739,7 @@ class Tally(object):
for filter_bin in filter_bins[i]:
bin_index = find_filter.get_bin_index(filter_bin)
if filter_type in ['energy', 'energyout']:
bin_indices.extend([bin_index])
bin_indices.extend([bin_index, bin_index+1])
num_bins += 1
elif filter_type == 'distribcell':
@ -2745,7 +2749,7 @@ class Tally(object):
bin_indices.append(bin_index)
num_bins += 1
find_filter.bins = find_filter.bins[bin_indices]
find_filter.bins = set(find_filter.bins[bin_indices])
find_filter.num_bins = num_bins
# Update the new tally's filter strides

View file

@ -397,28 +397,28 @@ class Cell(object):
def create_xml_subelement(self, xml_element):
element = ET.Element("cell")
element.set("id", str(self._id))
element.set("id", str(self.id))
if len(self._name) > 0:
element.set("name", str(self._name))
element.set("name", str(self.name))
if isinstance(self._fill, openmc.Material):
element.set("material", str(self._fill._id))
if isinstance(self.fill, basestring):
element.set("material", "void")
elif isinstance(self._fill, Iterable):
elif isinstance(self.fill, openmc.Material):
element.set("material", str(self.fill.id))
elif isinstance(self.fill, Iterable):
element.set("material", ' '.join([m if m == 'void' else str(m.id)
for m in self.fill]))
elif isinstance(self._fill, (Universe, Lattice)):
element.set("fill", str(self._fill._id))
self._fill.create_xml_subelement(xml_element)
elif self._fill.strip().lower() == "void":
element.set("material", "void")
elif isinstance(self.fill, (Universe, Lattice)):
element.set("fill", str(self.fill.id))
self.fill.create_xml_subelement(xml_element)
else:
element.set("fill", str(self._fill))
self._fill.create_xml_subelement(xml_element)
element.set("fill", str(self.fill))
self.fill.create_xml_subelement(xml_element)
if self.region is not None:
# Set the region attribute with the region specification
@ -445,11 +445,11 @@ class Cell(object):
# Call the recursive function from the top node
create_surface_elements(self.region, xml_element)
if self._translation is not None:
element.set("translation", ' '.join(map(str, self._translation)))
if self.translation is not None:
element.set("translation", ' '.join(map(str, self.translation)))
if self._rotation is not None:
element.set("rotation", ' '.join(map(str, self._rotation)))
if self.rotation is not None:
element.set("rotation", ' '.join(map(str, self.rotation)))
return element

View file

@ -1,17 +1,23 @@
module ace
use ace_header, only: Nuclide, Reaction, SAlphaBeta, XsListing, &
DistEnergy
use ace_header, only: Nuclide, Reaction, SAlphaBeta, XsListing
use constants
use endf, only: reaction_name, is_fission, is_disappearance
use error, only: fatal_error, warning
use fission, only: nu_total
use distribution_univariate, only: Uniform, Equiprobable, Tabular
use endf, only: reaction_name, is_fission, is_disappearance
use energy_distribution, only: TabularEquiprobable, LevelInelastic, &
ContinuousTabular, MaxwellEnergy, Evaporation, WattEnergy, NBodyPhaseSpace
use error, only: fatal_error, warning
use fission, only: nu_total
use global
use list_header, only: ListInt
use material_header, only: Material
use output, only: write_message
use set_header, only: SetChar
use string, only: to_str, to_lower
use list_header, only: ListInt
use material_header, only: Material
use output, only: write_message
use set_header, only: SetChar
use secondary_header, only: AngleEnergy
use secondary_correlated, only: CorrelatedAngleEnergy
use secondary_kalbach, only: KalbachMann
use secondary_uncorrelated, only: UncorrelatedAngleEnergy
use string, only: to_str, to_lower
implicit none
@ -372,8 +378,8 @@ contains
else
call read_nu_data(nuc)
call read_reactions(nuc)
call read_angular_dist(nuc)
call read_energy_dist(nuc)
call read_angular_dist(nuc)
call read_unr_res(nuc)
end if
@ -516,9 +522,10 @@ contains
integer :: LED ! location of energy distribution locators
integer :: LDIS ! location of all energy distributions
integer :: LOCC ! location of energy distributions for given MT
integer :: LAW
integer :: IDAT
integer :: lc ! locator
integer :: length ! length of data to allocate
type(DistEnergy), pointer :: edist
JXS2 = JXS(2)
JXS24 = JXS(24)
@ -661,11 +668,15 @@ contains
! Loop over all delayed neutron precursor groups
do i = 1, NPCR
! find location of energy distribution data
LOCC = int(XSS(LED + i - 1))
LOCC = nint(XSS(LED + i - 1))
! Determine law and location of data
LAW = nint(XSS(LDIS + LOCC))
IDAT = nint(XSS(LDIS + LOCC + 1))
! read energy distribution data
edist => nuc % nu_d_edist(i)
call get_energy_dist(edist, LOCC, .true.)
call get_energy_dist(nuc%nu_d_edist(i)%obj, LAW, LDIS, IDAT, &
ZERO, ZERO)
end do
! =======================================================================
@ -738,8 +749,8 @@ contains
rxn%multiplicity = 1
rxn%threshold = 1
rxn%scatter_in_cm = .true.
rxn%has_angle_dist = .false.
rxn%has_energy_dist = .false.
allocate(rxn%secondary%distribution(1))
allocate(UncorrelatedAngleEnergy :: rxn%secondary%distribution(1)%obj)
end associate
! Add contribution of elastic scattering to total cross section
@ -754,10 +765,6 @@ contains
do i = 1, NMT
associate (rxn => nuc % reactions(i+1))
! set defaults
rxn % has_angle_dist = .false.
rxn % has_energy_dist = .false.
! read MT number, Q-value, and neutrons produced
rxn % MT = int(XSS(LMT + i - 1))
rxn % Q_value = XSS(JXS4 + i - 1)
@ -887,86 +894,96 @@ contains
subroutine read_angular_dist(nuc)
type(Nuclide), intent(inout) :: nuc
integer :: JXS8 ! location of angular distribution locators
integer :: JXS9 ! location of angular distributions
integer :: LOCB ! location of angular distribution for given MT
integer :: NE ! number of incoming energies
integer :: NP ! number of points for cosine distribution
integer :: LC ! locator
integer :: i ! index in reactions array
integer :: j ! index over incoming energies
integer :: length ! length of data array to allocate
JXS8 = JXS(8)
JXS9 = JXS(9)
integer :: k ! index over energy distributions
integer :: interp
integer, allocatable :: LC(:) ! locator
! loop over all reactions with secondary neutrons -- NXS(5) does not include
! elastic scattering
do i = 1, NXS(5) + 1
associate (rxn => nuc%reactions(i))
! find location of angular distribution
LOCB = int(XSS(JXS8 + i - 1))
if (LOCB == -1) then
! Angular distribution data are specified through LAWi = 44 in the DLW
! block
cycle
elseif (LOCB == 0) then
! No angular distribution data are given for this reaction, isotropic
! scattering is assumed (in CM if TY < 0 and in LAB if TY > 0)
cycle
end if
rxn % has_angle_dist = .true.
LOCB = int(XSS(JXS(8) + i - 1))
! allocate space for incoming energies and locations
NE = int(XSS(JXS9 + LOCB - 1))
rxn % adist % n_energy = NE
allocate(rxn % adist % energy(NE))
allocate(rxn % adist % type(NE))
allocate(rxn % adist % location(NE))
! Angular distribution given as part of a correlated angle-energy distribution
if (LOCB == -1) cycle
! read incoming energy grid and location of nucs
XSS_index = JXS9 + LOCB
rxn % adist % energy = get_real(NE)
rxn % adist % location = get_int(NE)
! No angular distribution data are given for this reaction, isotropic
! scattering is assumed (in CM if TY < 0 and in LAB if TY > 0)
if (LOCB == 0) cycle
! determine dize of data block
length = 0
do j = 1, NE
LC = rxn % adist % location(j)
if (LC == 0) then
! isotropic
rxn % adist % type(j) = ANGLE_ISOTROPIC
elseif (LC > 0) then
! 32 equiprobable bins
rxn % adist % type(j) = ANGLE_32_EQUI
length = length + 33
elseif (LC < 0) then
! tabular distribution
rxn % adist % type(j) = ANGLE_TABULAR
NP = int(XSS(JXS9 + abs(LC)))
length = length + 2 + 3*NP
end if
end do
! Loop over each separate energy distribution. Even though there is only
! "one" angular distribution, it is repeated as many times as there are
! energy distributions for this reaction since the
! UncorrelatedAngleEnergy type holds one angle and energy distribution.
do k = 1, size(rxn%secondary%distribution)
select type (aedist => rxn%secondary%distribution(k)%obj)
type is (UncorrelatedAngleEnergy)
! allocate space for incoming energies and locations
NE = int(XSS(JXS(9) + LOCB - 1))
allocate(aedist%angle%energy(NE))
allocate(aedist%angle%distribution(NE))
allocate(LC(NE))
! allocate angular distribution data and read
allocate(rxn % adist % data(length))
! read incoming energy grid and location of nucs
XSS_index = JXS(9) + LOCB
aedist%angle%energy(:) = get_real(NE)
LC(:) = get_int(NE)
! read angular distribution -- currently this does not actually parse the
! angular distribution tables for each incoming energy, that must be done
! on-the-fly
XSS_index = JXS9 + LOCB + 2 * NE
rxn % adist % data = get_real(length)
! determine dize of data block
do j = 1, NE
if (LC(j) == 0) then
! isotropic
allocate(Uniform :: aedist%angle%distribution(j)%obj)
select type (adist => aedist%angle%distribution(j)%obj)
type is (Uniform)
adist%a = -ONE
adist%b = ONE
end select
! change location pointers since they are currently relative to JXS(9)
LC = LOCB + 2 * NE + 1
do j = 1, NE
! For consistency, leave location as 0 if type is isotropic.
! This is not necessary for current correctness, but can avoid
! future issues
if (rxn % adist % location(j) /= 0) then
rxn % adist % location(j) = abs(rxn % adist % location(j)) - LC
end if
elseif (LC(j) > 0) then
! 32 equiprobable bins
allocate(Equiprobable :: aedist%angle%distribution(j)%obj)
select type (adist => aedist%angle%distribution(j)%obj)
type is (Equiprobable)
allocate(adist%x(33))
end select
elseif (LC(j) < 0) then
! tabular distribution
allocate(Tabular :: aedist%angle%distribution(j)%obj)
end if
end do
! read angular distribution -- currently this does not actually parse the
! angular distribution tables for each incoming energy, that must be done
! on-the-fly
do j = 1, NE
XSS_index = JXS(9) + abs(LC(j)) - 1
select type(adist => aedist%angle%distribution(j)%obj)
type is (Equiprobable)
adist%x(:) = get_real(33)
type is (Tabular)
! determine interpolation and number of points
interp = nint(XSS(XSS_index))
NP = nint(XSS(XSS_index + 1))
! Get probability density data
XSS_index = XSS_index + 2
allocate(adist%x(NP), adist%p(NP), adist%c(NP))
adist%x(:) = get_real(NP)
adist%p(:) = get_real(NP)
adist%c(:) = get_real(NP)
end select
end do
deallocate(LC)
end select
end do
end associate
end do
@ -981,25 +998,62 @@ contains
subroutine read_energy_dist(nuc)
type(Nuclide), intent(inout) :: nuc
integer :: LED ! location of energy distribution locators
integer :: LOCC ! location of energy distributions for given MT
integer :: i ! loop index
LED = JXS(10)
integer :: n
integer :: IDAT ! locator for distribution data
integer :: LNW ! location of next energy law
integer :: LAW ! Type of energy law
! Loop over all reactions
do i = 1, NXS(5)
associate (rxn => nuc % reactions(i+1)) ! skip over elastic scattering
rxn % has_energy_dist = .true.
! Determine how many energy distributions are present for this reaction
LNW = nint(XSS(JXS(10) + i - 1))
n = 0
do while (LNW > 0)
n = n + 1
LNW = nint(XSS(JXS(11) + LNW - 1))
end do
! find location of energy distribution data
LOCC = int(XSS(LED + i - 1))
! Allocate space for distributions and probability of validity
associate (secondary => nuc%reactions(i + 1)%secondary)
allocate(secondary%applicability(n))
allocate(secondary%distribution(n))
! allocate energy distribution
allocate(rxn % edist)
LNW = nint(XSS(JXS(10) + i - 1))
n = 0
do while (LNW > 0)
n = n + 1
! read data for energy distribution
call get_energy_dist(rxn % edist, LOCC)
! Determine energy law and location of data
LAW = nint(XSS(JXS(11) + LNW))
IDAT = nint(XSS(JXS(11) + LNW + 1))
! Read probability of law validity
call secondary%applicability(n)%from_ace(XSS, JXS(11) + LNW + 2)
! Read energy law data
call get_energy_dist(secondary%distribution(n)%obj, LAW, &
JXS(11), IDAT, nuc%awr, nuc%reactions(i + 1)%Q_value)
! <<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<
! Before the secondary distribution refactor, when the angle/energy
! distribution was uncorrelated, no angle was actually sampled. With
! the refactor, an angle is always sampled for an uncorrelated
! distribution even when no angle distribution exists in the ACE file
! (isotropic is assumed). To preserve the RNG stream, we explicitly
! mark fission reactions so that we avoid the angle sampling.
if (any(nuc%reactions(i + 1)%MT == &
[N_FISSION, N_F, N_NF, N_2NF, N_3NF])) then
select type (aedist => secondary%distribution(n)%obj)
type is (UncorrelatedAngleEnergy)
aedist%fission = .true.
end select
end if
! <<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<
! Get locator for next distribution
LNW = nint(XSS(JXS(11) + LNW - 1))
end do
end associate
end do
@ -1011,281 +1065,319 @@ contains
! single reaction
!===============================================================================
recursive subroutine get_energy_dist(edist, loc_law, delayed_n)
type(DistEnergy), intent(inout) :: edist ! energy distribution
integer, intent(in) :: loc_law ! locator for data
logical, intent(in), optional :: delayed_n ! is this for delayed neutrons?
recursive subroutine get_energy_dist(aedist, law, LDIS, IDAT, awr, Q_value)
class(AngleEnergy), allocatable, intent(inout) :: aedist
integer, intent(in) :: law
integer, intent(in) :: LDIS
integer, intent(in) :: IDAT
real(8), intent(in) :: awr
real(8), intent(in) :: Q_value
integer :: LDIS ! location of all energy distributions
integer :: LNW ! location of next energy distribution if multiple
integer :: LAW ! secondary energy distribution law
integer :: i, j
integer :: NR ! number of interpolation regions
integer :: NE ! number of incoming energies
integer :: IDAT ! location of first energy distribution for given MT
integer :: lc ! locator
integer :: length ! length of data to allocate
integer :: length_interp_data ! length of interpolation data
integer :: NP ! number of outgoing energies/angles
integer :: interp
integer, allocatable :: L(:) ! locations of distributions for each Ein
integer, allocatable :: LC(:) ! locations of distributions for each Ein
! determine location of energy distribution
if (present(delayed_n)) then
LDIS = JXS(27)
XSS_index = LDIS + IDAT - 1
if (law == 44) then
allocate(KalbachMann :: aedist)
elseif (law == 61) then
allocate(CorrelatedAngleEnergy :: aedist)
else
LDIS = JXS(11)
allocate(UncorrelatedAngleEnergy :: aedist)
end if
! locator for next law and information on this law
LNW = int(XSS(LDIS + loc_law - 1))
LAW = int(XSS(LDIS + loc_law))
IDAT = int(XSS(LDIS + loc_law + 1))
NR = int(XSS(LDIS + loc_law + 2))
edist % law = LAW
edist % p_valid % n_regions = NR
select type (aedist)
type is (UncorrelatedAngleEnergy)
! ========================================================================
! UNCORRELATED ENERGY DISTRIBUTIONS
! allocate space for ENDF interpolation parameters
if (NR > 0) then
allocate(edist % p_valid % nbt(NR))
allocate(edist % p_valid % int(NR))
end if
! read ENDF interpolation parameters
XSS_index = LDIS + loc_law + 3
if (NR > 0) then
edist % p_valid % nbt = int(get_real(NR))
edist % p_valid % int = int(get_real(NR))
end if
! allocate space for law validity data
NE = int(XSS(LDIS + loc_law + 3 + 2*NR))
edist % p_valid % n_pairs = NE
allocate(edist % p_valid % x(NE))
allocate(edist % p_valid % y(NE))
length_interp_data = 5 + 2*(NR + NE)
! read law validity data
XSS_index = LDIS + loc_law + 4 + 2*NR
edist % p_valid % x = get_real(NE)
edist % p_valid % y = get_real(NE)
! Set index to beginning of IDAT array
lc = LDIS + IDAT - 2
! determine length of energy distribution
length = length_energy_dist(lc, LAW, loc_law, length_interp_data)
! allocate secondary energy distribution array
allocate(edist % data(length))
! read secondary energy distribution
XSS_index = lc + 1
edist % data = get_real(length)
! read next energy distribution if present
if (LNW > 0) then
allocate(edist % next)
call get_energy_dist(edist % next, LNW)
end if
end subroutine get_energy_dist
!===============================================================================
! LENGTH_ENERGY_DIST determines how many values are contained in an LDAT energy
! distribution array based on the secondary energy law and location in XSS
!===============================================================================
function length_energy_dist(lc, law, LOCC, lid) result(length)
integer, intent(in) :: lc ! location in XSS array
integer, intent(in) :: law ! energy distribution law
integer, intent(in) :: LOCC ! location of energy distribution
integer, intent(in) :: lid ! length of interpolation data
integer :: length ! length of energy distribution (LDAT)
integer :: i ! loop index for incoming energies
integer :: j ! loop index for outgoing energies
integer :: k ! dummy index in XSS
integer :: NR ! number of interpolation regions
integer :: NE ! number of incoming energies
integer :: NP ! number of points in outgoing energy distribution
integer :: NMU ! number of points in outgoing cosine distribution
integer :: NRa ! number of interpolation regions for Watt 'a'
integer :: NEa ! number of energies for Watt 'a'
integer :: NRb ! number of interpolation regions for Watt 'b'
integer :: NEb ! number of energies for Watt 'b'
real(8), allocatable :: L(:) ! locations of distributions for each Ein
! initialize length
length = 0
select case (law)
case (1)
! Tabular equiprobable energy bins
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
NP = int(XSS(lc + 3 + 2*NR + NE))
length = 3 + 2*NR + NE + 3*NP*NE
case (2)
! Discrete photon energy
length = 2
case (3)
! Level scattering
length = 2
case (4)
! Continuous tabular distribution
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
allocate(L(NE))
L(:) = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1))
! Continue with finding data length
length = length + 2 + 2*NR + 2*NE
do i = 1,NE
! Some older data sets use the same LDAT for multiple Ein tables.
! If this is the case, we should skip incrementing length when it is
! not needed.
if (i < NE) then
if (any(L(i) == L(i + 1: NE))) then
! adjust location for this block
j = lc + 2 + 2*NR + NE + i
XSS(j) = XSS(j) - LOCC - lid
cycle
select case (law)
case (1)
allocate(TabularEquiprobable :: aedist%energy)
select type (edist => aedist%energy)
type is (TabularEquiprobable)
NR = nint(XSS(XSS_index))
NE = nint(XSS(XSS_index + 1 + 2*NR))
if (NR > 0) then
call fatal_error("Multiple interpolation regions not yet supported &
&for tabular equiprobable energy distributions.")
end if
end if
! determine length
NP = int(XSS(lc + length + 2))
length = length + 2 + 3*NP
edist%n_region = NR
! adjust location for this block
j = lc + 2 + 2*NR + NE + i
XSS(j) = XSS(j) - LOCC - lid
! Read incoming energies for which outgoing energies are tabulated
allocate(edist%energy_in(NE))
XSS_index = XSS_index + 2 + 2*NR
edist%energy_in(:) = get_real(NE)
! Read outgoing energy tables
NP = nint(XSS(XSS_index))
allocate(edist%energy_out(NP, NE))
XSS_index = XSS_index + 1
do i = 1, NE
edist%energy_out(:, i) = get_real(NP)
end do
end select
case (3)
allocate(LevelInelastic :: aedist%energy)
select type (edist => aedist%energy)
type is (LevelInelastic)
edist%threshold = XSS(XSS_index)
edist%mass_ratio = XSS(XSS_index + 1)
end select
case (4)
allocate(ContinuousTabular :: aedist%energy)
select type (edist => aedist%energy)
type is (ContinuousTabular)
NR = nint(XSS(XSS_index))
XSS_index = XSS_index + 1
if (NR > 1) then
call fatal_error("Multiple interpolation regions not yet supported &
&for continuous tabular energy distributions.")
end if
edist%n_region = NR
! Read breakpoints and interpolation parameters
if (NR > 0) then
allocate(edist%breakpoints(NR))
allocate(edist%interpolation(NR))
edist%breakpoints(:) = get_int(NR)
edist%interpolation(:) = get_int(NR)
end if
! Read incoming energies for which outgoing energies are tabulated and
! locators
NE = nint(XSS(XSS_index))
XSS_index = XSS_index + 1
allocate(edist%energy_in(NE))
allocate(L(NE))
edist%energy_in(:) = get_real(NE)
L(:) = get_int(NE)
! Read outgoing energy tables
allocate(edist%energy_out(NE))
do i = 1, NE
! Determine interpolation and number of discrete points
XSS_index = LDIS + L(i) - 1
interp = nint(XSS(XSS_index))
edist%energy_out(i)%interpolation = mod(interp, 10)
edist%energy_out(i)%n_discrete = (interp - &
edist%energy_out(i)%interpolation)/10
! check for discrete lines present
if (edist%energy_out(i)%n_discrete > 0) then
call fatal_error("Discrete lines in continuous tabular &
&distribution not yet supported")
end if
! Determine number of points and allocate space
NP = nint(XSS(XSS_index + 1))
allocate(edist%energy_out(i)%e_out(NP))
allocate(edist%energy_out(i)%p(NP))
allocate(edist%energy_out(i)%c(NP))
! Read tabular PDF for outgoing energy
XSS_index = XSS_index + 2
edist%energy_out(i)%e_out(:) = get_real(NP)
edist%energy_out(i)%p(:) = get_real(NP)
edist%energy_out(i)%c(:) = get_real(NP)
end do
deallocate(L)
end select
case (7)
allocate(MaxwellEnergy :: aedist%energy)
select type (edist => aedist%energy)
type is (MaxwellEnergy)
call edist%theta%from_ace(XSS, XSS_index)
edist%u = XSS(XSS_index + 2 + 2*edist%theta%n_regions + &
2*edist%theta%n_pairs)
end select
case (9)
allocate(Evaporation :: aedist%energy)
select type(edist => aedist%energy)
type is (Evaporation)
call edist%theta%from_ace(XSS, XSS_index)
edist%u = XSS(XSS_index + 2 + 2*edist%theta%n_regions + &
2*edist%theta%n_pairs)
end select
case (11)
allocate(WattEnergy :: aedist%energy)
select type(edist => aedist%energy)
type is (WattEnergy)
call edist%a%from_ace(XSS, XSS_index)
XSS_index = XSS_index + 2 + 2*edist%a%n_regions + 2*edist%a%n_pairs
call edist%b%from_ace(XSS, XSS_index)
XSS_index = XSS_index + 2 + 2*edist%b%n_regions + 2*edist%b%n_pairs
edist%u = XSS(XSS_index)
end select
case (66)
allocate(NBodyPhaseSpace :: aedist%energy)
select type(edist => aedist%energy)
type is (NBodyPhaseSpace)
edist%n_bodies = int(XSS(XSS_index))
edist%mass_ratio = XSS(XSS_index + 1)
edist%A = awr
edist%Q = Q_value
end select
end select
type is (KalbachMann)
! ========================================================================
! CORRELATED KALBACH-MANN DISTRIBUTION
NR = int(XSS(XSS_index))
NE = int(XSS(XSS_index + 1 + 2*NR))
if (NR > 0) then
call fatal_error("Multiple interpolation regions not yet supported &
&for Kalbach-Mann energy distributions.")
end if
aedist%n_region = NR
! Read incoming energies for which outgoing energies are tabulated and locators
allocate(aedist%energy_in(NE))
allocate(L(NE))
XSS_index = XSS_index + 2 + 2*NR
aedist%energy_in(:) = get_real(NE)
L(:) = get_int(NE)
! Read outgoing energy tables
allocate(aedist%table(NE))
do i = 1, NE
! Determine interpolation and number of discrete points
XSS_index = LDIS + L(i) - 1
interp = nint(XSS(XSS_index))
aedist%table(i)%interpolation = mod(interp, 10)
aedist%table(i)%n_discrete = (interp - aedist%table(i)%interpolation)/10
! check for discrete lines present
if (aedist%table(i)%n_discrete > 0) then
call fatal_error("Discrete lines in Kalbach-Mann distribution not &
&yet supported")
end if
! Determine number of points and allocate space
NP = nint(XSS(XSS_index + 1))
allocate(aedist%table(i)%e_out(NP))
allocate(aedist%table(i)%p(NP))
allocate(aedist%table(i)%c(NP))
allocate(aedist%table(i)%r(NP))
allocate(aedist%table(i)%a(NP))
! Read tabular PDF for outgoing energy
XSS_index = XSS_index + 2
aedist%table(i)%e_out(:) = get_real(NP)
aedist%table(i)%p(:) = get_real(NP)
aedist%table(i)%c(:) = get_real(NP)
aedist%table(i)%r(:) = get_real(NP)
aedist%table(i)%a(:) = get_real(NP)
end do
deallocate(L)
case (5)
! General evaporation spectrum
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
NP = int(XSS(lc + 3 + 2*NR + 2*NE))
length = 3 + 2*NR + 2*NE + NP
type is (CorrelatedAngleEnergy)
! ========================================================================
! CORRELATED ANGLE-ENERGY DISTRIBUTION
case (7)
! Maxwell fission spectrum
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
length = 3 + 2*NR + 2*NE
NR = int(XSS(XSS_index))
NE = int(XSS(XSS_index + 1 + 2*NR))
if (NR > 0) then
call fatal_error("Multiple interpolation regions not yet supported &
&for correlated angle-energy distributions.")
end if
aedist%n_region = NR
case (9)
! Evaporation spectrum
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
length = 3 + 2*NR + 2*NE
case (11)
! Watt spectrum
NRa = int(XSS(lc + 1))
NEa = int(XSS(lc + 2 + 2*NRa))
NRb = int(XSS(lc + 3 + 2*(NRa+NEa)))
NEb = int(XSS(lc + 4 + 2*(NRa+NEa+NRb)))
length = 5 + 2*(NRa + NEa + NRb + NEb)
case (44)
! Kalbach-Mann correlated scattering
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
! Read incoming energies for which outgoing energies are tabulated and
! locators
allocate(aedist%energy_in(NE))
allocate(L(NE))
L(:) = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1))
XSS_index = XSS_index + 2 + 2*NR
aedist%energy_in(:) = get_real(NE)
L(:) = get_int(NE)
! Continue with finding data length
length = length + 2 + 2*NR + 2*NE
do i = 1,NE
! Some older data sets use the same LDAT for multiple Ein tables.
! If this is the case, we should skip incrementing length when it is
! not needed.
if (i < NE) then
if (any(L(i) == L(i + 1: NE))) then
! adjust location for this block
j = lc + 2 + 2*NR + NE + i
XSS(j) = XSS(j) - LOCC - lid
cycle
end if
end if
NP = int(XSS(lc + length + 2))
length = length + 2 + 5*NP
! Read outgoing energy tables
allocate(aedist%table(NE))
do i = 1, NE
! Determine interpolation and number of discrete points
XSS_index = LDIS + L(i) - 1
interp = nint(XSS(XSS_index))
aedist%table(i)%interpolation = mod(interp, 10)
aedist%table(i)%n_discrete = (interp - aedist%table(i)%interpolation)/10
! adjust location for this block
j = lc + 2 + 2*NR + NE + i
XSS(j) = XSS(j) - LOCC - lid
end do
deallocate(L)
case (61)
! Correlated energy and angle distribution
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
allocate(L(NE))
L(:) = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1))
! Continue with finding data length
length = length + 2 + 2*NR + 2*NE
do i = 1,NE
! Some older data sets use the same LDAT for multiple Ein tables.
! If this is the case, we should skip incrementing length when it is
! not needed.
if (i < NE) then
if (any(L(i) == L(i + 1: NE))) then
! adjust locators for energy distribution
j = lc + 2 + 2*NR + NE + i
XSS(j) = XSS(j) - LOCC - lid
cycle
end if
! check for discrete lines present
if (aedist%table(i)%n_discrete > 0) then
call fatal_error("Discrete lines in correlated angle-energy &
&distribution not yet supported")
end if
! outgoing energy distribution
NP = int(XSS(lc + length + 2))
! Determine number of points and allocate space
NP = nint(XSS(XSS_index + 1))
allocate(aedist%table(i)%e_out(NP))
allocate(aedist%table(i)%p(NP))
allocate(aedist%table(i)%c(NP))
allocate(LC(NP))
! adjust locators for angular distribution
! Read tabular PDF for outgoing energy
XSS_index = XSS_index + 2
aedist%table(i)%e_out(:) = get_real(NP)
aedist%table(i)%p(:) = get_real(NP)
aedist%table(i)%c(:) = get_real(NP)
LC(:) = get_int(NP)
! allocate angular distributions for each incoming/outgoing energy
allocate(aedist%table(i)%angle(NP))
do j = 1, NP
k = lc + length + 2 + 3*NP + j
if (XSS(k) /= 0) XSS(k) = XSS(k) - LOCC - lid
if (LC(j) == 0) then
! isotropic
allocate(Uniform :: aedist%table(i)%angle(j)%obj)
select type (adist => aedist%table(i)%angle(j)%obj)
type is (Uniform)
adist%a = -ONE
adist%b = ONE
end select
elseif (LC(j) > 0) then
! tabular distribution
allocate(Tabular :: aedist%table(i)%angle(j)%obj)
end if
end do
length = length + 2 + 4*NP
! read angular distributions
do j = 1, NP
! outgoing angle distribution -- NMU here is actually
! referred to as NP in the MCNP documentation
NMU = int(XSS(lc + length + 2))
length = length + 2 + 3*NMU
XSS_index = LDIS + abs(LC(j)) - 1
select type(adist => aedist%table(i)%angle(j)%obj)
type is (Tabular)
! determine interpolation and number of points
interp = nint(XSS(XSS_index))
NP = nint(XSS(XSS_index + 1))
! Get probability density data
XSS_index = XSS_index + 2
allocate(adist%x(NP), adist%p(NP), adist%c(NP))
adist%x(:) = get_real(NP)
adist%p(:) = get_real(NP)
adist%c(:) = get_real(NP)
end select
end do
deallocate(LC)
! adjust locators for energy distribution
j = lc + 2 + 2*NR + NE + i
XSS(j) = XSS(j) - LOCC - lid
end do
deallocate(L)
case (66)
! N-body phase space distribution
length = 2
case (67)
! Laboratory energy-angle law
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
! Before progressing, check to see if data set uses L(I) values
! in a way inconsistent with the current form of the ACE Format Guide
! (MCNP5 Manual, Vol 3)
allocate(L(NE))
L(:) = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1))
! Don't currently do anything with L
deallocate(L)
! Continue with finding data length
NMU = int(XSS(lc + 4 + 2*NR + 2*NE))
length = 4 + 2*(NR + NE + NMU)
end select
end function length_energy_dist
end subroutine get_energy_dist
!===============================================================================
! READ_UNR_RES reads in unresolved resonance probability tables if present.

View file

@ -3,42 +3,11 @@ module ace_header
use constants, only: MAX_FILE_LEN, ZERO
use dict_header, only: DictIntInt
use endf_header, only: Tab1
use secondary_header, only: SecondaryDistribution, AngleEnergyContainer
use stl_vector, only: VectorInt
implicit none
!===============================================================================
! DISTANGLE contains data for a tabular secondary angle distribution whether it
! be tabular or 32 equiprobable cosine bins
!===============================================================================
type DistAngle
integer :: n_energy ! # of incoming energies
real(8), allocatable :: energy(:) ! incoming energy grid
integer, allocatable :: type(:) ! type of distribution
integer, allocatable :: location(:) ! location of each table
real(8), allocatable :: data(:) ! angular distribution data
end type DistAngle
!===============================================================================
! DISTENERGY contains data for a secondary energy distribution for all
! scattering laws
!===============================================================================
type DistEnergy
integer :: law ! secondary distribution law
type(Tab1) :: p_valid ! probability of law validity
real(8), allocatable :: data(:) ! energy distribution data
! For reactions that may have multiple energy distributions such as (n,2n),
! this pointer allows multiple laws to be stored
type(DistEnergy), pointer :: next => null()
! Type-Bound procedures
contains
procedure :: clear => distenergy_clear ! Deallocates DistEnergy
end type DistEnergy
!===============================================================================
! REACTION contains the cross-section and secondary energy and angle
! distributions for a single reaction in a continuous-energy ACE-format table
@ -53,10 +22,7 @@ module ace_header
logical :: scatter_in_cm ! scattering system in center-of-mass?
logical :: multiplicity_with_E = .false. ! Flag to indicate E-dependent multiplicity
real(8), allocatable :: sigma(:) ! Cross section values
logical :: has_angle_dist ! Angle distribution present?
logical :: has_energy_dist ! Energy distribution present?
type(DistAngle) :: adist ! Secondary angular distribution
type(DistEnergy), pointer :: edist => null() ! Secondary energy distribution
type(SecondaryDistribution) :: secondary
! Type-Bound procedures
contains
@ -137,7 +103,7 @@ module ace_header
integer :: n_precursor ! # of delayed neutron precursors
real(8), allocatable :: nu_d_data(:)
real(8), allocatable :: nu_d_precursor_data(:)
type(DistEnergy), pointer :: nu_d_edist(:) => null()
type(AngleEnergyContainer), allocatable :: nu_d_edist(:)
! Unresolved resonance data
logical :: urr_present
@ -284,37 +250,14 @@ module ace_header
contains
!===============================================================================
! DISTENERGY_CLEAR resets and deallocates data in DistEnergy.
!===============================================================================
recursive subroutine distenergy_clear(this)
class(DistEnergy), intent(inout) :: this ! The DistEnergy object to clear
if (associated(this % next)) then
! recursively clear this item
call this % next % clear()
deallocate(this % next)
end if
end subroutine distenergy_clear
!===============================================================================
! REACTION_CLEAR resets and deallocates data in Reaction.
!===============================================================================
subroutine reaction_clear(this)
class(Reaction), intent(inout) :: this ! The Reaction object to clear
if (associated(this % multiplicity_E)) deallocate(this % multiplicity_E)
if (associated(this % edist)) then
call this % edist % clear()
deallocate(this % edist)
end if
end subroutine reaction_clear
!===============================================================================
@ -322,21 +265,11 @@ module ace_header
!===============================================================================
subroutine nuclide_clear(this)
class(Nuclide), intent(inout) :: this ! The Nuclide object to clear
class(Nuclide), intent(inout) :: this
integer :: i ! Loop counter
if (associated(this % nu_d_edist)) then
do i = 1, size(this % nu_d_edist)
call this % nu_d_edist(i) % clear()
end do
deallocate(this % nu_d_edist)
end if
if (associated(this % urr_data)) then
deallocate(this % urr_data)
end if
if (associated(this % urr_data)) deallocate(this % urr_data)
if (allocated(this % reactions)) then
do i = 1, size(this % reactions)

View file

@ -0,0 +1,63 @@
module angle_distribution
use constants, only: ZERO, ONE
use distribution_univariate, only: DistributionContainer
use random_lcg, only: prn
use search, only: binary_search
implicit none
private
!===============================================================================
! ANGLEDISTRIBUTION represents an angular distribution that is to be used in an
! uncorrelated angle-energy distribution. This occurs whenever the angle
! distrbution is given in File 4 in an ENDF file. The distribution of angles
! depends on the incoming energy of the neutron, so this type stores a
! distribution for each of a set of incoming energies.
!===============================================================================
type, public :: AngleDistribution
real(8), allocatable :: energy(:)
type(DistributionContainer), allocatable :: distribution(:)
contains
procedure :: sample => angle_sample
end type AngleDistribution
contains
function angle_sample(this, E) result(mu)
class(AngleDistribution), intent(in) :: this
real(8), intent(in) :: E ! incoming energy
real(8) :: mu ! sampled cosine of scattering angle
integer :: i ! index on incoming energy grid
integer :: n ! number of incoming energies
real(8) :: r ! interpolation factor on incoming energy grid
! Determine number of incoming energies
n = size(this%energy)
! Find energy bin and calculate interpolation factor -- if the energy is
! outside the range of the tabulated energies, choose the first or last bins
if (E < this%energy(1)) then
i = 1
r = ZERO
elseif (E > this%energy(n)) then
i = n - 1
r = ONE
else
i = binary_search(this%energy, n, E)
r = (E - this%energy(i))/(this%energy(i+1) - this%energy(i))
end if
! Sample between the ith and (i+1)th bin
if (r > prn()) i = i + 1
! Sample i-th distribution
mu = this%distribution(i)%obj%sample()
! Make sure mu is in range [-1,1]
if (abs(mu) > ONE) mu = sign(ONE, mu)
end function angle_sample
end module angle_distribution

View file

@ -1,7 +1,7 @@
module distribution_univariate
use constants, only: ZERO, HALF, HISTOGRAM, LINEAR_LINEAR, MAX_LINE_LEN, &
MAX_WORD_LEN
use constants, only: ZERO, ONE, HALF, HISTOGRAM, LINEAR_LINEAR, &
MAX_LINE_LEN, MAX_WORD_LEN
use error, only: fatal_error
use math, only: maxwell_spectrum, watt_spectrum
use random_lcg, only: prn
@ -78,6 +78,12 @@ module distribution_univariate
procedure :: initialize => tabular_initialize
end type Tabular
type, extends(Distribution) :: Equiprobable
real(8), allocatable :: x(:)
contains
procedure :: sample => equiprobable_sample
end type Equiprobable
contains
function discrete_sample(this) result(x)
@ -238,6 +244,25 @@ contains
this%c(:) = this%c(:)/this%c(n)
end subroutine tabular_initialize
function equiprobable_sample(this) result(x)
class(Equiprobable), intent(in) :: this
real(8) :: x
integer :: i
integer :: n
real(8) :: r
real(8) :: xl, xr
n = size(this%x)
r = prn()
i = 1 + int((n - 1)*r)
xl = this%x(i)
xr = this%x(i+1)
x = xl + ((n - 1)*r - i + ONE) * (xr - xl)
end function equiprobable_sample
subroutine distribution_from_xml(dist, node_dist)
class(Distribution), allocatable, intent(inout) :: dist
type(Node), pointer :: node_dist

View file

@ -13,6 +13,40 @@ module endf_header
integer :: n_pairs ! # of pairs of (x,y) values
real(8), allocatable :: x(:) ! values of abscissa
real(8), allocatable :: y(:) ! values of ordinate
contains
procedure :: from_ace
end type Tab1
contains
subroutine from_ace(this, xss, idx)
class(Tab1), intent(inout) :: this
real(8), intent(in) :: xss(:)
integer, intent(in) :: idx
integer :: nr, ne
! Determine number of regions
nr = nint(xss(idx))
this%n_regions = nr
! Read interpolation region data
if (nr > 0) then
allocate(this%nbt(nr))
allocate(this%int(nr))
this%nbt(:) = nint(xss(idx + 1 : idx + nr))
this%int(:) = nint(xss(idx + nr + 1 : idx + 2*nr))
end if
! Determine number of pairs
ne = int(XSS(idx + 2*nr + 1))
this%n_pairs = ne
! Read (x,y) pairs
allocate(this%x(ne))
allocate(this%y(ne))
this%x(:) = xss(idx + 2*nr + 2 : idx + 2*nr + 1 + ne)
this%y(:) = xss(idx + 2*nr + 2 + ne : idx + 2*nr + 1 + 2*ne)
end subroutine from_ace
end module endf_header

429
src/energy_distribution.F90 Normal file
View file

@ -0,0 +1,429 @@
module energy_distribution
use constants, only: ZERO, ONE, TWO, PI, HISTOGRAM, LINEAR_LINEAR
use endf_header, only: Tab1
use interpolation, only: interpolate_tab1
use math, only: maxwell_spectrum, watt_spectrum
use random_lcg, only: prn
use search, only: binary_search
!===============================================================================
! ENERGYDISTRIBUTION (abstract) defines an energy distribution that is a
! function of the incident energy of a projectile. Each derived type must
! implement a sample() function that returns a sampled outgoing energy given an
! incoming energy
!===============================================================================
type, abstract :: EnergyDistribution
contains
procedure(iSampleEnergy), deferred :: sample
end type EnergyDistribution
abstract interface
function iSampleEnergy(this, E_in) result(E_out)
import EnergyDistribution
class(EnergyDistribution), intent(in) :: this
real(8), intent(in) :: E_in
real(8) :: E_out
end function iSampleEnergy
end interface
type :: EnergyDistributionContainer
class(EnergyDistribution), allocatable :: obj
end type EnergyDistributionContainer
!===============================================================================
! Derived classes
!===============================================================================
!===============================================================================
! TABULAREQUIPROBABLE represents an energy distribution with tabular
! equiprobable energy bins as given in ACE law 1. This is an older
! representation that has largely been replaced with ACE laws 4, 44, and 61.
!===============================================================================
type, extends(EnergyDistribution) :: TabularEquiprobable
integer :: n_region ! number of interpolation regions
integer, allocatable :: breakpoints(:) ! breakpoints of interpolation regions
integer, allocatable :: interpolation(:) ! interpolation region codes
real(8), allocatable :: energy_in(:) ! incoming energies
real(8), allocatable :: energy_out(:,:) ! table of outgoing energies for
! each incoming energy
contains
procedure :: sample => equiprobable_sample
end type TabularEquiprobable
!===============================================================================
! LEVELINELASTIC gives the energy distribution for level inelastic scattering by
! neutrons as in ENDF MT=51--90.
!===============================================================================
type, extends(EnergyDistribution) :: LevelInelastic
real(8) :: threshold
real(8) :: mass_ratio
contains
procedure :: sample => level_inelastic_sample
end type LevelInelastic
!===============================================================================
! CONTINUOUSTABULAR gives an energy distribution represented as a tabular
! distribution with histogram or linear-linear interpolation. This corresponds
! to ACE law 4, which NJOY produces for a number of ENDF energy distributions.
!===============================================================================
type CTTable
integer :: interpolation
integer :: n_discrete
real(8), allocatable :: e_out(:)
real(8), allocatable :: p(:)
real(8), allocatable :: c(:)
end type CTTable
type, extends(EnergyDistribution) :: ContinuousTabular
integer :: n_region
integer, allocatable :: breakpoints(:)
integer, allocatable :: interpolation(:)
real(8), allocatable :: energy_in(:)
type(CTTable), allocatable :: energy_out(:)
contains
procedure :: sample => continuous_sample
end type ContinuousTabular
!===============================================================================
! MAXWELLENERGY gives the energy distribution of neutrons emitted from a Maxwell
! fission spectrum. This corresponds to ACE law 7 and ENDF File 5, LF=7.
!===============================================================================
type, extends(EnergyDistribution) :: MaxwellEnergy
type(Tab1) :: theta ! incoming-energy-dependent parameter
real(8) :: u ! restriction energy
contains
procedure :: sample => maxwellenergy_sample
end type MaxwellEnergy
!===============================================================================
! EVAPORATION represents an evaporation spectrum corresponding to ACE law 9 and
! ENDF File 5, LF=9.
!===============================================================================
type, extends(EnergyDistribution) :: Evaporation
type(Tab1) :: theta
real(8) :: u
contains
procedure :: sample => evaporation_sample
end type Evaporation
!===============================================================================
! WATTENERGY gives the energy distribution of neutrons emitted from a Watt
! fission spectrum. This corresponds to ACE law 11 and ENDF File 5, LF=11.
!===============================================================================
type, extends(EnergyDistribution) :: WattEnergy
type(Tab1) :: a
type(Tab1) :: b
real(8) :: u
contains
procedure :: sample => watt_sample
end type WattEnergy
!===============================================================================
! NBODYPHASESPACE gives the energy distribution for particles emitted from
! neutron and charged-particle reactions. This corresponds to ACE law 66 and
! ENDF File 6, LAW=6.
!===============================================================================
type, extends(EnergyDistribution) :: NBodyPhaseSpace
integer :: n_bodies
real(8) :: mass_ratio
real(8) :: A
real(8) :: Q
contains
procedure :: sample => nbody_sample
end type NBodyPhaseSpace
contains
function equiprobable_sample(this, E_in) result(E_out)
class(TabularEquiprobable), intent(in) :: this
real(8), intent(in) :: E_in ! incoming energy
real(8) :: E_out ! sampled outgoing energy
integer :: i, k, l ! indices
integer :: n_energy_in ! number of incoming energies
integer :: n_energy_out ! number of outgoing energies
real(8) :: r ! interpolation factor on incoming energy
real(8) :: E_i_1, E_i_K ! endpoints on outgoing grid i
real(8) :: E_i1_1, E_i1_K ! endpoints on outgoing grid i+1
real(8) :: E_1, E_K ! endpoints interpolated between i and i+1
real(8) :: E_l_k, E_l_k1 ! adjacent E on outgoing grid l
! Determine number of incoming/outgoing energies
n_energy_in = size(this%energy_in)
n_energy_out = size(this%energy_out, 1)
! Determine index on incoming energy grid and interpolation factor
i = binary_search(this%energy_in, size(this%energy_in), E_in)
r = (E_in - this%energy_in(i)) / &
(this%energy_in(i+1) - this%energy_in(i))
! Sample outgoing energy bin
k = 1 + int(n_energy_out * prn())
! Determine E_1 and E_K
E_i_1 = this%energy_out(1, i)
E_i_K = this%energy_out(n_energy_out, i)
E_i1_1 = this%energy_out(1, i+1)
E_i1_K = this%energy_out(n_energy_out, i+1)
E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
E_K = E_i_K + r*(E_i1_K - E_i_K)
! Randomly select between the outgoing table for incoming energy E_i and
! E_(i+1)
if (prn() < r) then
l = i + 1
else
l = i
end if
! Determine E_l_k and E_l_k+1
E_l_k = this%energy_out(k, l)
E_l_k1 = this%energy_out(k+1, l)
! Determine E' (denoted here as E_out)
E_out = E_l_k + prn()*(E_l_k1 - E_l_k)
! Now interpolate between incident energy bins i and i + 1
if (l == i) then
E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1)
else
E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1)
end if
end function equiprobable_sample
function level_inelastic_sample(this, E_in) result(E_out)
class(LevelInelastic), intent(in) :: this
real(8), intent(in) :: E_in
real(8) :: E_out
E_out = this%mass_ratio*(E_in - this%threshold)
end function level_inelastic_sample
function continuous_sample(this, E_in) result(E_out)
class(ContinuousTabular), intent(in) :: this
real(8), intent(in) :: E_in ! incoming energy
real(8) :: E_out ! sampled outgoing energy
integer :: i, k, l ! indices
integer :: n_energy_in ! number of incoming energies
integer :: n_energy_out ! number of outgoing energies
real(8) :: r ! interpolation factor on incoming energy
real(8) :: r1 ! random number on [0,1)
real(8) :: frac ! interpolation factor on outgoing energy
real(8) :: E_i_1, E_i_K ! endpoints on outgoing grid i
real(8) :: E_i1_1, E_i1_K ! endpoints on outgoing grid i+1
real(8) :: E_1, E_K ! endpoints interpolated between i and i+1
real(8) :: E_l_k, E_l_k1 ! adjacent E on outgoing grid l
real(8) :: p_l_k, p_l_k1 ! adjacent p on outgoing grid l
real(8) :: c_k, c_k1 ! cumulative probability
logical :: histogram_interp ! whether histogram interpolation is used
! Read number of interpolation regions and incoming energies
if (this%n_region == 1) then
histogram_interp = (this%interpolation(1) == 1)
else
histogram_interp = .false.
end if
! Find energy bin and calculate interpolation factor -- if the energy is
! outside the range of the tabulated energies, choose the first or last bins
n_energy_in = size(this%energy_in)
if (E_in < this%energy_in(1)) then
i = 1
r = ZERO
elseif (E_in > this%energy_in(n_energy_in)) then
i = n_energy_in - 1
r = ONE
else
i = binary_search(this%energy_in, n_energy_in, E_in)
r = (E_in - this%energy_in(i)) / &
(this%energy_in(i+1) - this%energy_in(i))
end if
! Sample between the ith and (i+1)th bin
if (histogram_interp) then
l = i
else
if (r > prn()) then
l = i + 1
else
l = i
end if
end if
! Interpolation for energy E1 and EK
n_energy_out = size(this%energy_out(i)%e_out)
E_i_1 = this%energy_out(i)%e_out(1)
E_i_K = this%energy_out(i)%e_out(n_energy_out)
n_energy_out = size(this%energy_out(i+1)%e_out)
E_i1_1 = this%energy_out(i+1)%e_out(1)
E_i1_K = this%energy_out(i+1)%e_out(n_energy_out)
E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
E_K = E_i_K + r*(E_i1_K - E_i_K)
! Determine outgoing energy bin
n_energy_out = size(this%energy_out(l)%e_out)
r1 = prn()
c_k = this%energy_out(l)%c(1)
do k = 1, n_energy_out - 1
c_k1 = this%energy_out(l)%c(k+1)
if (r1 < c_k1) exit
c_k = c_k1
end do
! Check to make sure k is <= NP - 1
k = min(k, n_energy_out - 1)
E_l_k = this%energy_out(l)%e_out(k)
p_l_k = this%energy_out(l)%p(k)
if (this%energy_out(l)%interpolation == HISTOGRAM) then
! Histogram interpolation
if (p_l_k > ZERO) then
E_out = E_l_k + (r1 - c_k)/p_l_k
else
E_out = E_l_k
end if
elseif (this%energy_out(l)%interpolation == LINEAR_LINEAR) then
! Linear-linear interpolation
E_l_k1 = this%energy_out(l)%e_out(k+1)
p_l_k1 = this%energy_out(l)%p(k+1)
frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k)
if (frac == ZERO) then
E_out = E_l_k + (r1 - c_k)/p_l_k
else
E_out = E_l_k + (sqrt(max(ZERO, p_l_k*p_l_k + &
TWO*frac*(r1 - c_k))) - p_l_k)/frac
end if
end if
! Now interpolate between incident energy bins i and i + 1
if (.not. histogram_interp) then
if (l == i) then
E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1)
else
E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1)
end if
end if
end function continuous_sample
function maxwellenergy_sample(this, E_in) result(E_out)
class(MaxwellEnergy), intent(in) :: this
real(8), intent(in) :: E_in ! incoming energy
real(8) :: E_out ! sampled outgoing energy
real(8) :: theta ! Maxwell distribution parameter
! Get temperature corresponding to incoming energy
theta = interpolate_tab1(this%theta, E_in)
do
! Sample maxwell fission spectrum
E_out = maxwell_spectrum(theta)
! Accept energy based on restriction energy
if (E_out <= E_in - this%u) exit
end do
end function maxwellenergy_sample
function evaporation_sample(this, E_in) result(E_out)
class(Evaporation), intent(in) :: this
real(8), intent(in) :: E_in ! incoming energy
real(8) :: E_out ! sampled outgoing energy
real(8) :: theta ! evaporation spectrum parameter
real(8) :: x, y, v
! Get temperature corresponding to incoming energy
theta = interpolate_tab1(this%theta, E_in)
y = (E_in - this%U)/theta
v = 1 - exp(-y)
! Sample outgoing energy based on evaporation spectrum probability
! density function
do
x = -log((ONE - v*prn())*(ONE - v*prn()))
if (x <= y) exit
end do
E_out = x*theta
end function evaporation_sample
function watt_sample(this, E_in) result(E_out)
class(WattEnergy), intent(in) :: this
real(8), intent(in) :: E_in ! incoming energy
real(8) :: E_out ! sampled outgoing energy
real(8) :: a, b ! Watt spectrum parameters
! Determine Watt parameter 'a' from tabulated function
a = interpolate_tab1(this%a, E_in)
! Determine Watt parameter 'b' from tabulated function
b = interpolate_tab1(this%b, E_in)
do
! Sample energy-dependent Watt fission spectrum
E_out = watt_spectrum(a, b)
! Accept energy based on restriction energy
if (E_out <= E_in - this%u) exit
end do
end function watt_sample
function nbody_sample(this, E_in) result(E_out)
class(NBodyPhaseSpace), intent(in) :: this
real(8), intent(in) :: E_in ! incoming energy
real(8) :: E_out ! sampled outgoing energy
real(8) :: Ap ! total mass of particles in neutron masses
real(8) :: E_max ! maximum possible COM energy
real(8) :: x, y, v
real(8) :: r1, r2, r3, r4, r5, r6
! Determine E_max parameter
Ap = this%mass_ratio
E_max = (Ap - ONE)/Ap * (this%A/(this%A + ONE)*E_in + this%Q)
! x is essentially a Maxwellian distribution
x = maxwell_spectrum(ONE)
select case (this%n_bodies)
case (3)
y = maxwell_spectrum(ONE)
case (4)
r1 = prn()
r2 = prn()
r3 = prn()
y = -log(r1*r2*r3)
case (5)
r1 = prn()
r2 = prn()
r3 = prn()
r4 = prn()
r5 = prn()
r6 = prn()
y = -log(r1*r2*r3*r4) - log(r5) * cos(PI/TWO*r6)**2
end select
! Now determine v and E_out
v = x/(x+y)
E_out = E_max * v
end function nbody_sample
end module energy_distribution

View file

@ -22,7 +22,6 @@ module global
#endif
implicit none
save
! ============================================================================
! GEOMETRY-RELATED VARIABLES

View file

@ -2,7 +2,6 @@ module interpolation
use constants
use endf_header, only: Tab1
use error, only: fatal_error
use search, only: binary_search
use string, only: to_str

View file

@ -322,14 +322,8 @@ contains
integer :: i ! loop index over nuclides
integer :: unit_ ! unit to write to
integer :: size_total ! memory used by nuclide (bytes)
integer :: size_angle_total ! total memory used for angle dist. (bytes)
integer :: size_energy_total ! total memory used for energy dist. (bytes)
integer :: size_xs ! memory used for cross-sections (bytes)
integer :: size_angle ! memory used for an angle distribution (bytes)
integer :: size_energy ! memory used for a energy distributions (bytes)
integer :: size_urr ! memory used for probability tables (bytes)
character(11) :: law ! secondary energy distribution law
type(UrrData), pointer :: urr
! set default unit for writing information
@ -340,8 +334,6 @@ contains
end if
! Initialize totals
size_angle_total = 0
size_energy_total = 0
size_urr = 0
size_xs = 0
@ -356,33 +348,15 @@ contains
write(unit_,*) ' # of reactions = ' // trim(to_str(nuc % n_reaction))
! Information on each reaction
write(unit_,*) ' Reaction Q-value COM Law IE size(angle) size(energy)'
write(unit_,*) ' Reaction Q-value COM IE'
do i = 1, nuc % n_reaction
associate (rxn => nuc % reactions(i))
! Determine size of angle distribution
if (rxn % has_angle_dist) then
size_angle = rxn % adist % n_energy * 16 + size(rxn % adist % data) * 8
else
size_angle = 0
end if
! Determine size of energy distribution and law
if (rxn % has_energy_dist) then
size_energy = size(rxn % edist % data) * 8
law = to_str(rxn % edist % law)
else
size_energy = 0
law = 'None'
end if
write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,A4,1X,I6,1X,I11,1X,I11)') &
write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,I6)') &
reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, &
law(1:4), rxn % threshold, size_angle, size_energy
rxn % threshold
! Accumulate data size
size_xs = size_xs + (nuc % n_grid - rxn%threshold + 1) * 8
size_angle_total = size_angle_total + size_angle
size_energy_total = size_energy_total + size_energy
end associate
end do
@ -408,19 +382,11 @@ contains
size_urr = urr % n_energy * (urr % n_prob * 6 + 1) * 8
end if
! Calculate total memory
size_total = size_xs + size_angle_total + size_energy_total + size_urr
! Write memory used
write(unit_,*) ' Memory Requirements'
write(unit_,*) ' Cross sections = ' // trim(to_str(size_xs)) // ' bytes'
write(unit_,*) ' Secondary angle distributions = ' // &
trim(to_str(size_angle_total)) // ' bytes'
write(unit_,*) ' Secondary energy distributions = ' // &
trim(to_str(size_energy_total)) // ' bytes'
write(unit_,*) ' Probability Tables = ' // &
trim(to_str(size_urr)) // ' bytes'
write(unit_,*) ' Total = ' // trim(to_str(size_total)) // ' bytes'
! Blank line at end of nuclide
write(unit_,*)

File diff suppressed because it is too large Load diff

View file

@ -19,8 +19,8 @@ element tallies {
(element id { xsd:int } | attribute id { xsd:int }) &
(element name { xsd:string { maxLength="52" } } |
attribute name { xsd:string { maxLength="52" } })? &
(element estimator { ( "analog" | "tracklength" ) } |
attribute estimator { ( "analog" | "tracklength" ) })? &
(element estimator { ( "analog" | "tracklength" | "collision" ) } |
attribute estimator { ( "analog" | "tracklength" | "collision" ) })? &
element filter {
(element type { ( "cell" | "cellborn" | "material" | "universe" |
"surface" | "distribcell" | "mesh" | "energy" | "energyout" | "mu" |

View file

@ -120,12 +120,14 @@
<choice>
<value>analog</value>
<value>tracklength</value>
<value>collision</value>
</choice>
</element>
<attribute name="estimator">
<choice>
<value>analog</value>
<value>tracklength</value>
<value>collision</value>
</choice>
</attribute>
</choice>

View file

@ -1,7 +1,6 @@
module search
use constants
use error, only: fatal_error
implicit none

View file

@ -0,0 +1,148 @@
module secondary_correlated
use constants, only: ZERO, ONE, TWO, HISTOGRAM, LINEAR_LINEAR
use distribution_univariate, only: DistributionContainer
use secondary_header, only: AngleEnergy
use random_lcg, only: prn
use search, only: binary_search
!===============================================================================
! CORRELATEDANGLEENERGY represents a correlated angle-energy distribution. This
! corresponds to ACE law 61 and ENDF File 6, LAW=1, LANG/=2.
!===============================================================================
type AngleEnergyTable
integer :: interpolation
integer :: n_discrete
real(8), allocatable :: e_out(:)
real(8), allocatable :: p(:)
real(8), allocatable :: c(:)
type(DistributionContainer), allocatable :: angle(:)
end type AngleEnergyTable
type, extends(AngleEnergy) :: CorrelatedAngleEnergy
integer :: n_region ! number of interpolation regions
integer, allocatable :: breakpoints(:) ! breakpoints of interpolation regions
integer, allocatable :: interpolation(:) ! interpolation region codes
real(8), allocatable :: energy_in(:) ! incoming energies
type(AngleEnergyTable), allocatable :: table(:) ! outgoing E/mu distributions
contains
procedure :: sample => correlated_sample
end type CorrelatedAngleEnergy
contains
subroutine correlated_sample(this, E_in, E_out, mu)
class(CorrelatedAngleEnergy), intent(in) :: this
real(8), intent(in) :: E_in ! incoming energy
real(8), intent(out) :: E_out ! sampled outgoing energy
real(8), intent(out) :: mu ! sapmled scattering cosine
integer :: i, k, l ! indices
integer :: n_energy_in ! number of incoming energies
integer :: n_energy_out ! number of outgoing energies
real(8) :: r ! interpolation factor on incoming energy
real(8) :: r1 ! random number on [0,1)
real(8) :: frac ! interpolation factor on outgoing energy
real(8) :: E_i_1, E_i_K ! endpoints on outgoing grid i
real(8) :: E_i1_1, E_i1_K ! endpoints on outgoing grid i+1
real(8) :: E_1, E_K ! endpoints interpolated between i and i+1
real(8) :: E_l_k, E_l_k1 ! adjacent E on outgoing grid l
real(8) :: p_l_k, p_l_k1 ! adjacent p on outgoing grid l
real(8) :: c_k, c_k1 ! cumulative probability
! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! Before the secondary distribution refactor, an isotropic polar cosine was
! always sampled but then overwritten with the polar cosine sampled from the
! correlated distribution. To preserve the random number stream, we keep
! this dummy sampling here but can remove it later (will change answers)
mu = TWO*prn() - ONE
! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! find energy bin and calculate interpolation factor -- if the energy is
! outside the range of the tabulated energies, choose the first or last bins
n_energy_in = size(this%energy_in)
if (E_in < this%energy_in(1)) then
i = 1
r = ZERO
elseif (E_in > this%energy_in(n_energy_in)) then
i = n_energy_in - 1
r = ONE
else
i = binary_search(this%energy_in, n_energy_in, E_in)
r = (E_in - this%energy_in(i)) / &
(this%energy_in(i+1) - this%energy_in(i))
end if
! Sample between the ith and (i+1)th bin
if (r > prn()) then
l = i + 1
else
l = i
end if
! interpolation for energy E1 and EK
n_energy_out = size(this%table(i)%e_out)
E_i_1 = this%table(i)%e_out(1)
E_i_K = this%table(i)%e_out(n_energy_out)
n_energy_out = size(this%table(i+1)%e_out)
E_i1_1 = this%table(i+1)%e_out(1)
E_i1_K = this%table(i+1)%e_out(n_energy_out)
E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
E_K = E_i_K + r*(E_i1_K - E_i_K)
! determine outgoing energy bin
n_energy_out = size(this%table(l)%e_out)
r1 = prn()
c_k = this%table(l)%c(1)
do k = 1, n_energy_out - 1
c_k1 = this%table(l)%c(k+1)
if (r1 < c_k1) exit
c_k = c_k1
end do
! check to make sure k is <= NP - 1
k = min(k, n_energy_out - 1)
E_l_k = this%table(l)%e_out(k)
p_l_k = this%table(l)%p(k)
if (this%table(l)%interpolation == HISTOGRAM) then
! Histogram interpolation
if (p_l_k > ZERO) then
E_out = E_l_k + (r1 - c_k)/p_l_k
else
E_out = E_l_k
end if
elseif (this%table(l)%interpolation == LINEAR_LINEAR) then
! Linear-linear interpolation
E_l_k1 = this%table(l)%e_out(k+1)
p_l_k1 = this%table(l)%p(k+1)
frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k)
if (frac == ZERO) then
E_out = E_l_k + (r1 - c_k)/p_l_k
else
E_out = E_l_k + (sqrt(max(ZERO, p_l_k*p_l_k + &
TWO*frac*(r1 - c_k))) - p_l_k)/frac
end if
end if
! Now interpolate between incident energy bins i and i + 1
if (l == i) then
E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1)
else
E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1)
end if
! Find correlated angular distribution for closest outgoing energy bin
if (r1 - c_k < c_k1 - r1) then
mu = this%table(l)%angle(k)%obj%sample()
else
mu = this%table(l)%angle(k + 1)%obj%sample()
end if
end subroutine correlated_sample
end module secondary_correlated

78
src/secondary_header.F90 Normal file
View file

@ -0,0 +1,78 @@
module secondary_header
use endf_header, only: Tab1
use interpolation, only: interpolate_tab1
use random_lcg, only: prn
!===============================================================================
! ANGLEENERGY (abstract) defines a correlated or uncorrelated angle-energy
! distribution that is a function of incoming energy. Each derived type must
! implement a sample() subroutine that returns an outgoing energy and scattering
! cosine given an incoming energy.
!===============================================================================
type, abstract :: AngleEnergy
contains
procedure(iSampleAngleEnergy), deferred :: sample
end type AngleEnergy
abstract interface
subroutine iSampleAngleEnergy(this, E_in, E_out, mu)
import AngleEnergy
class(AngleEnergy), intent(in) :: this
real(8), intent(in) :: E_in
real(8), intent(out) :: E_out
real(8), intent(out) :: mu
end subroutine iSampleAngleEnergy
end interface
type :: AngleEnergyContainer
class(AngleEnergy), allocatable :: obj
end type AngleEnergyContainer
!===============================================================================
! SECONDARYDISTRIBUTION stores multiple angle-energy distributions, each of
! which has a given probability of occurring for a given incoming energy. In
! general, most secondary distributions only have one angle-energy distribution,
! but for some cases (e.g., (n,2n) in certain nuclides) multiple distinct
! distributions exist.
!===============================================================================
type :: SecondaryDistribution
type(Tab1), allocatable :: applicability(:)
type(AngleEnergyContainer), allocatable :: distribution(:)
contains
procedure :: sample => secondary_sample
end type SecondaryDistribution
contains
subroutine secondary_sample(this, E_in, E_out, mu)
class(SecondaryDistribution), intent(in) :: this
real(8), intent(in) :: E_in ! incoming energy
real(8), intent(out) :: E_out ! sampled outgoing energy
real(8), intent(out) :: mu ! sampled scattering cosine
integer :: n ! number of angle-energy distributions
real(8) :: p_valid ! probability that given distribution is valid
n = size(this%applicability)
if (n > 1) then
do i = 1, n
! Determine probability that i-th energy distribution is sampled
p_valid = interpolate_tab1(this%applicability(i), E_in)
! If i-th distribution is sampled, sample energy from the distribution
if (prn() <= p_valid) then
call this%distribution(i)%obj%sample(E_in, E_out, mu)
exit
end if
end do
else
! If only one distribution is present, go ahead and sample it
call this%distribution(1)%obj%sample(E_in, E_out, mu)
end if
end subroutine secondary_sample
end module secondary_header

164
src/secondary_kalbach.F90 Normal file
View file

@ -0,0 +1,164 @@
module secondary_kalbach
use constants, only: ZERO, ONE, TWO, HISTOGRAM, LINEAR_LINEAR
use secondary_header, only: AngleEnergy
use random_lcg, only: prn
use search, only: binary_search
!===============================================================================
! KalbachMann represents a correlated angle-energy distribution with the angular
! distribution represented using Kalbach-Mann systematics. This corresponds to
! ACE law 44 and ENDF File 6, LAW=1, LANG=2.
!===============================================================================
type KalbachMannTable
integer :: n_discrete
integer :: interpolation
real(8), allocatable :: e_out(:)
real(8), allocatable :: p(:)
real(8), allocatable :: c(:)
real(8), allocatable :: r(:)
real(8), allocatable :: a(:)
end type KalbachMannTable
type, extends(AngleEnergy) :: KalbachMann
integer :: n_region ! number of interpolation regions
integer, allocatable :: breakpoints(:) ! breakpoints of interpolation regions
integer, allocatable :: interpolation(:) ! interpolation region codes
real(8), allocatable :: energy_in(:) ! incoming energies
type(KalbachMannTable), allocatable :: table(:) ! outgoing E/mu parameters
contains
procedure :: sample => kalbachmann_sample
end type KalbachMann
contains
subroutine kalbachmann_sample(this, E_in, E_out, mu)
class(KalbachMann), intent(in) :: this
real(8), intent(in) :: E_in ! incoming energy
real(8), intent(out) :: E_out ! sampled outgoing energy
real(8), intent(out) :: mu ! sampled scattering cosine
integer :: i, k, l ! indices
integer :: n_energy_in ! number of incoming energies
integer :: n_energy_out ! number of outgoing energies
real(8) :: r ! interpolation factor on incoming energy
real(8) :: r1 ! random number on [0,1)
real(8) :: frac ! interpolation factor on outgoing energy
real(8) :: E_i_1, E_i_K ! endpoints on outgoing grid i
real(8) :: E_i1_1, E_i1_K ! endpoints on outgoing grid i+1
real(8) :: E_1, E_K ! endpoints interpolated between i and i+1
real(8) :: E_l_k, E_l_k1 ! adjacent E on outgoing grid l
real(8) :: p_l_k, p_l_k1 ! adjacent p on outgoing grid l
real(8) :: c_k, c_k1 ! cumulative probability
real(8) :: km_r, km_a ! Kalbach-Mann parameters
real(8) :: T
! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! Before the secondary distribution refactor, an isotropic polar cosine was
! always sampled but then overwritten with the polar cosine sampled from the
! correlated distribution. To preserve the random number stream, we keep
! this dummy sampling here but can remove it later (will change answers)
mu = TWO*prn() - ONE
! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! find energy bin and calculate interpolation factor -- if the energy is
! outside the range of the tabulated energies, choose the first or last bins
n_energy_in = size(this%energy_in)
if (E_in < this%energy_in(1)) then
i = 1
r = ZERO
elseif (E_in > this%energy_in(n_energy_in)) then
i = n_energy_in - 1
r = ONE
else
i = binary_search(this%energy_in, n_energy_in, E_in)
r = (E_in - this%energy_in(i)) / &
(this%energy_in(i+1) - this%energy_in(i))
end if
! Sample between the ith and (i+1)th bin
if (r > prn()) then
l = i + 1
else
l = i
end if
! interpolation for energy E1 and EK
n_energy_out = size(this%table(i)%e_out)
E_i_1 = this%table(i)%e_out(1)
E_i_K = this%table(i)%e_out(n_energy_out)
n_energy_out = size(this%table(i+1)%e_out)
E_i1_1 = this%table(i+1)%e_out(1)
E_i1_K = this%table(i+1)%e_out(n_energy_out)
E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
E_K = E_i_K + r*(E_i1_K - E_i_K)
! determine outgoing energy bin
n_energy_out = size(this%table(l)%e_out)
r1 = prn()
c_k = this%table(l)%c(1)
do k = 1, n_energy_out - 1
c_k1 = this%table(l)%c(k+1)
if (r1 < c_k1) exit
c_k = c_k1
end do
! check to make sure k is <= NP - 1
k = min(k, n_energy_out - 1)
E_l_k = this%table(l)%e_out(k)
p_l_k = this%table(l)%p(k)
if (this%table(l)%interpolation == HISTOGRAM) then
! Histogram interpolation
if (p_l_k > ZERO) then
E_out = E_l_k + (r1 - c_k)/p_l_k
else
E_out = E_l_k
end if
! Determine Kalbach-Mann parameters
km_r = this%table(l)%r(k)
km_a = this%table(l)%a(k)
elseif (this%table(l)%interpolation == LINEAR_LINEAR) then
! Linear-linear interpolation
E_l_k1 = this%table(l)%e_out(k+1)
p_l_k1 = this%table(l)%p(k+1)
frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k)
if (frac == ZERO) then
E_out = E_l_k + (r1 - c_k)/p_l_k
else
E_out = E_l_k + (sqrt(max(ZERO, p_l_k*p_l_k + &
TWO*frac*(r1 - c_k))) - p_l_k)/frac
end if
! Determine Kalbach-Mann parameters
km_r = this%table(l)%r(k) + (E_out - E_l_k)/(E_l_k1 - E_l_k) * &
(this%table(l)%r(k+1) - this%table(l)%r(k))
km_a = this%table(l)%a(k) + (E_out - E_l_k)/(E_l_k1 - E_l_k) * &
(this%table(l)%a(k+1) - this%table(l)%a(k))
end if
! Now interpolate between incident energy bins i and i + 1
if (l == i) then
E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1)
else
E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1)
end if
! Sampled correlated angle from Kalbach-Mann parameters
if (prn() > km_r) then
T = (TWO*prn() - ONE) * sinh(km_a)
mu = log(T + sqrt(T*T + ONE))/km_a
else
r1 = prn()
mu = log(r1*exp(km_a) + (ONE - r1)*exp(-km_a))/km_a
end if
end subroutine kalbachmann_sample
end module secondary_kalbach

View file

@ -0,0 +1,48 @@
module secondary_uncorrelated
use angle_distribution, only: AngleDistribution
use constants, only: ONE, TWO
use energy_distribution, only: EnergyDistribution
use secondary_header, only: AngleEnergy
use random_lcg, only: prn
!===============================================================================
! UNCORRELATEDANGLEENERGY represents an uncorrelated angle-energy
! distribution. This corresponds to when an energy distribution is given in ENDF
! File 5/6 and an angular distribution is given in ENDF File 4.
!===============================================================================
type, extends(AngleEnergy) :: UncorrelatedAngleEnergy
logical :: fission = .false.
type(AngleDistribution) :: angle
class(EnergyDistribution), allocatable :: energy
contains
procedure :: sample => uncorrelated_sample
end type UncorrelatedAngleEnergy
contains
subroutine uncorrelated_sample(this, E_in, E_out, mu)
class(UncorrelatedAngleEnergy), intent(in) :: this
real(8), intent(in) :: E_in ! incoming energy
real(8), intent(out) :: E_out ! sampled outgoing energy
real(8), intent(out) :: mu ! sampled scattering cosine
! Sample cosine of scattering angle
if (this%fission) then
! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! For fission, the angle is not used, so just assign a dummy value
mu = ONE
! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
elseif (allocated(this%angle%energy)) then
mu = this%angle%sample(E_in)
else
! no angle distribution given => assume isotropic for all energies
mu = TWO*prn() - ONE
end if
! Sample outgoing energy
E_out = this%energy%sample(E_in)
end subroutine uncorrelated_sample
end module secondary_uncorrelated

View file

@ -45,7 +45,7 @@ contains
call write_attribute_string(file_id, "n_batches", &
"description", "Total number of batches")
! Write eigenvalue information
! Write eigenvalue information
if (run_mode == MODE_EIGENVALUE) then
! write number of inactive/active batches and generations/batch
call write_dataset(file_id, "n_inactive", n_inactive)
@ -174,8 +174,10 @@ contains
case (CELL_FILL)
call write_dataset(cell_group, "fill_type", "universe")
call write_dataset(cell_group, "fill", universes(c%fill)%id)
if (size(c%offset) > 0) then
call write_dataset(cell_group, "offset", c%offset)
if (allocated(c%offset)) then
if (size(c%offset) > 0) then
call write_dataset(cell_group, "offset", c%offset)
end if
end if
if (allocated(c%translation)) then
@ -362,8 +364,10 @@ contains
call write_dataset(lattice_group, "outer", lat%outer)
! Write distribcell offsets if present
if (size(lat%offset) > 0) then
call write_dataset(lattice_group, "offsets", lat%offset)
if (allocated(lat%offset)) then
if (size(lat%offset) > 0) then
call write_dataset(lattice_group, "offsets", lat%offset)
end if
end if
select type (lat)

View file

@ -1,10 +0,0 @@
#!/bin/bash
# This simple script ensures that all binary
# output files have been deleted in all the
# folders. This can occur if a previous error
# occurred and the test suite was rerun without
# deleting left over binary files. This will
# cause an assertion error in some of the
# tests.
find . \( -name "*.h5" -o -name "*.ppm" \) -exec rm -f {} \;

View file

@ -8,7 +8,7 @@ import shutil
import re
import glob
import socket
from subprocess import call
from subprocess import call, check_output
from collections import OrderedDict
from optparse import OptionParser
@ -42,9 +42,9 @@ parser.add_option("-s", "--script", action="store_true", dest="script",
# Default compiler paths
FC='gfortran'
MPI_DIR='/opt/mpich/3.1.3-gnu'
HDF5_DIR='/opt/hdf5/1.8.15-gnu'
PHDF5_DIR='/opt/phdf5/1.8.15-gnu'
MPI_DIR='/opt/mpich/3.2-gnu'
HDF5_DIR='/opt/hdf5/1.8.16-gnu'
PHDF5_DIR='/opt/phdf5/1.8.16-gnu'
# Script mode for extra capability
script_mode = False
@ -73,11 +73,13 @@ set(CTEST_UPDATE_COMMAND "git")
set(CTEST_CONFIGURE_COMMAND "${{CMAKE_COMMAND}} -H${{CTEST_SOURCE_DIRECTORY}} -B${{CTEST_BINARY_DIRECTORY}} ${{CTEST_BUILD_OPTIONS}}")
set(CTEST_MEMORYCHECK_COMMAND "{valgrind_cmd}")
set(CTEST_MEMORYCHECK_COMMAND_OPTIONS "--tool=memcheck --leak-check=yes --show-reachable=yes --num-callers=20 --track-fds=yes")
set(CTEST_MEMORYCHECK_SUPPRESSIONS_FILE ${{CTEST_SOURCE_DIRECTORY}}/../tests/valgrind.supp)
#set(CTEST_MEMORYCHECK_SUPPRESSIONS_FILE ${{CTEST_SOURCE_DIRECTORY}}/../tests/valgrind.supp)
set(MEM_CHECK {mem_check})
if(MEM_CHECK)
set(ENV{{MEM_CHECK}} ${{MEM_CHECK}})
endif()
set(CTEST_COVERAGE_COMMAND "{gcov_cmd}")
set(CTEST_COVERAGE_COMMAND "gcov")
set(COVERAGE {coverage})
set(ENV{{COVERAGE}} ${{COVERAGE}})
@ -87,9 +89,11 @@ ctest_start("{dashboard}")
ctest_configure(RETURN_VALUE res)
{update}
ctest_build(RETURN_VALUE res)
if(NOT MEM_CHECK)
ctest_test({tests} PARALLEL_LEVEL {n_procs}, RETURN_VALUE res)
endif()
if(MEM_CHECK)
ctest_memcheck({tests}, RETURN_VALUE res)
ctest_memcheck({tests} RETURN_VALUE res)
endif(MEM_CHECK)
if(COVERAGE)
ctest_coverage(RETURN_VALUE res)
@ -105,6 +109,32 @@ endif()
# Define test data structure
tests = OrderedDict()
def cleanup(path):
"""Remove generated output files."""
for dirpath, dirnames, filenames in os.walk(path):
for fname in filenames:
for ext in ['.h5', '.ppm', '.voxel']:
if fname.endswith(ext):
os.remove(os.path.join(dirpath, fname))
def which(program):
def is_exe(fpath):
return os.path.isfile(fpath) and os.access(fpath, os.X_OK)
fpath, fname = os.path.split(program)
if fpath:
if is_exe(program):
return program
else:
for path in os.environ["PATH"].split(os.pathsep):
path = path.strip('"')
exe_file = os.path.join(path, program)
if is_exe(exe_file):
return exe_file
return None
class Test(object):
def __init__(self, name, debug=False, optimize=False, mpi=False, openmp=False,
phdf5=False, valgrind=False, coverage=False):
@ -119,8 +149,6 @@ class Test(object):
self.success = True
self.msg = None
self.skipped = False
self.valgrind_cmd = ""
self.gcov_cmd = ""
self.cmake = ['cmake', '-H..', '-Bbuild',
'-DPYTHON_EXECUTABLE=' + sys.executable]
@ -231,42 +259,6 @@ class Test(object):
self.success = False
self.msg = 'Failed on testing.'
# Checks to see if file exists in PWD or PATH
def check_compiler(self):
result = False
if os.path.isfile(self.fc):
result = True
for path in os.environ["PATH"].split(":"):
if os.path.isfile(os.path.join(path, self.fc)):
result = True
if not result:
self.msg = 'Compiler not found: {0}'.\
format((os.path.join(path, self.fc)))
self.success = False
# Get valgrind command from user's environment
def find_valgrind(self):
result = False
for path in os.environ["PATH"].split(":"):
if os.path.isfile(os.path.join(path, 'valgrind')):
self.valgrind_cmd = os.path.join(path, 'valgrind')
result = True
break
if not result:
self.msg = 'valgrind not found.'
self.success = False
# Get coverage command from user's environment
def find_coverage(self):
result = False
for path in os.environ["PATH"].split(":"):
if os.path.isfile(os.path.join(path, 'gcov')):
self.gcov_cmd = os.path.join(path, 'gcov')
result = True
break
if not result:
self.msg = 'gcov not found.'
self.success = False
# Simple function to add a test to the global tests dictionary
def add_test(name, debug=False, optimize=False, mpi=False, openmp=False,\
@ -342,7 +334,7 @@ else:
# Setup CTest script vars. Not used in non-script mode
pwd = os.getcwd()
ctest_vars = {
'source_dir': os.path.join(pwd, '..'),
'source_dir': os.path.join(pwd, os.pardir),
'build_dir': os.path.join(pwd, 'build'),
'host_name': socket.gethostname(),
'dashboard': dash,
@ -363,10 +355,10 @@ else:
# Set up default valgrind tests (subset of all tests)
# Currently takes too long to run all the tests with valgrind
# Only used in script mode
valgrind_default_tests = "basic|cmfd_feed|confidence_intervals|\
density_atombcm|eigenvalue_genperbatch|energy_grid|entropy|\
filter_cell|lattice_multiple|output|plot_background|reflective_plane|\
rotation|salphabeta_multiple|score_absorption|seed|source_energy_mono|\
valgrind_default_tests = "cmfd_feed|confidence_intervals|\
density|eigenvalue_genperbatch|energy_grid|entropy|\
lattice_multiple|output|plotreflective_plane|\
rotation|salphabetascore_absorption|seed|source_energy_mono|\
sourcepoint_batch|statepoint_interval|survival_biasing|\
tally_assumesep|translation|uniform_fs|universe|void"
@ -383,7 +375,7 @@ if len(list(tests.keys())) == 0:
# Begin testing
shutil.rmtree('build', ignore_errors=True)
call(['./cleanup']) # removes all binary and hdf5 output files from tests
cleanup('.')
for key in iter(tests):
test = tests[key]
@ -395,29 +387,34 @@ for key in iter(tests):
sys.stdout.flush()
# Verify fortran compiler exists
test.check_compiler()
if not test.success:
if which(test.fc) is None:
self.msg = 'Compiler not found: {0}'.format(test.fc)
self.success = False
continue
# Get valgrind command
# Verify valgrind command exists
if test.valgrind:
test.find_valgrind()
if not test.success:
continue
valgrind_cmd = which('valgrind')
if valgrind_cmd is None:
self.msg = 'No valgrind executable found.'
self.success = False
continue
else:
valgrind_cmd = ''
# Get coverage command
# Verify gcov/lcov exist
if test.coverage:
test.find_coverage()
if not test.success:
continue
if which('gcov') is None:
self.msg = 'No {} executable found.'.format(exe)
self.success = False
continue
# Set test specific CTest script vars. Not used in non-script mode
ctest_vars.update({'build_name' : test.get_build_name()})
ctest_vars.update({'build_opts' : test.get_build_opts()})
ctest_vars.update({'mem_check' : test.valgrind})
ctest_vars.update({'coverage' : test.coverage})
ctest_vars.update({'valgrind_cmd' : test.valgrind_cmd})
ctest_vars.update({'gcov_cmd' : test.gcov_cmd})
ctest_vars.update({'build_name': test.get_build_name()})
ctest_vars.update({'build_opts': test.get_build_opts()})
ctest_vars.update({'mem_check': test.valgrind})
ctest_vars.update({'coverage': test.coverage})
ctest_vars.update({'valgrind_cmd': valgrind_cmd})
# Check for user custom tests
# INCLUDE is a CTest command that allows for a subset
@ -458,7 +455,7 @@ for key in iter(tests):
test.run_ctests()
# Leave build directory
os.chdir('..')
os.chdir(os.pardir)
# Copy over log file
if script_mode:
@ -471,11 +468,37 @@ for key in iter(tests):
logfilename = logfilename + '_{0}.log'.format(test.name)
shutil.copy(logfile[0], logfilename)
# For coverage builds, use lcov to generate HTML output
if test.coverage:
if which('lcov') is None or which('genhtml') is None:
print('No lcov/genhtml command found. '
'Could not generate coverage report.')
else:
shutil.rmtree('coverage', ignore_errors=True)
call(['lcov', '--directory', '.', '--capture',
'--output-file', 'coverage.info'])
call(['genhtml', '--output-directory', 'coverage', 'coverage.info'])
os.remove('coverage.info')
if test.valgrind:
# Copy memcheck output to memcheck directory
shutil.rmtree('memcheck', ignore_errors=True)
os.mkdir('memcheck')
memcheck_out = glob.glob('build/Testing/Temporary/MemoryChecker.*.log')
for fname in memcheck_out:
shutil.copy(fname, 'memcheck/')
# Remove generated XML files
xml_files = check_output(['git', 'ls-files', '.', '--exclude-standard',
'--others']).split()
for f in xml_files:
os.remove(f)
# Clear build directory and remove binary and hdf5 files
shutil.rmtree('build', ignore_errors=True)
if script_mode:
os.remove('ctestscript.run')
call(['./cleanup'])
cleanup('.')
# Print out summary of results
print('\n' + '='*54)

View file

@ -1,8 +0,0 @@
<?xml version="1.0"?>
<geometry>
<!-- Sphere with radius 10 -->
<surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
<cell id="1" material="1" region="-1" />
</geometry>

View file

@ -1,9 +0,0 @@
<?xml version="1.0"?>
<materials>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U-235" xs="71c" ao="1.0" />
</material>
</materials>

View file

@ -1,2 +0,0 @@
k-combined:
3.021779E-01 3.813358E-03

View file

@ -0,0 +1,14 @@
<?xml version="1.0"?>
<geometry>
<surface id="1" type="sphere" coeffs="0 0 0 3"/>
<surface id="2" type="sphere" coeffs="0 0 0 6"/>
<surface id="3" type="sphere" coeffs="0 0 0 9"/>
<surface id="4" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
<cell id="1" material="1" region="-1" />
<cell id="2" material="2" region="1 -2" />
<cell id="3" material="3" region="2 -3" />
<cell id="4" material="4" region="3 -4" />
</geometry>

View file

@ -0,0 +1,26 @@
<?xml version="1.0"?>
<materials>
<material id="1">
<density value="0.1" units="atom/b-cm" />
<nuclide name="U-235" xs="71c" ao="1.0" />
</material>
<material id="2">
<density value="4.5e22" units="atom/cm3" />
<nuclide name="U-235" xs="71c" ao="1.0" />
</material>
<material id="3">
<density value="12.3e3" units="kg/m3" />
<nuclide name="U-235" xs="71c" ao="1.0" />
</material>
<material id="4">
<density units="sum" />
<nuclide name="U-235" xs="71c" ao="0.3e-2" />
<nuclide name="U-238" xs="71c" ao="0.5e-1" />
<nuclide name="H-1" xs="71c" ao="0.1e-2" />
</material>
</materials>

View file

@ -0,0 +1,2 @@
k-combined:
1.088237E+00 1.999252E-02

View file

View file

@ -1,8 +0,0 @@
<?xml version="1.0"?>
<geometry>
<!-- Sphere with radius 10 -->
<surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
<cell id="1" material="1" region="-1" />
</geometry>

View file

@ -1,9 +0,0 @@
<?xml version="1.0"?>
<materials>
<material id="1">
<density value="0.1" units="atom/b-cm" />
<nuclide name="U-235" xs="71c" ao="1.0" />
</material>
</materials>

View file

@ -1,2 +0,0 @@
k-combined:
1.752274E+00 4.032481E-02

View file

@ -1,11 +0,0 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness
if __name__ == '__main__':
harness = TestHarness('statepoint.10.*')
harness.main()

View file

@ -1,8 +0,0 @@
<?xml version="1.0"?>
<geometry>
<!-- Sphere with radius 10 -->
<surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
<cell id="1" material="1" region="-1" />
</geometry>

View file

@ -1,9 +0,0 @@
<?xml version="1.0"?>
<materials>
<material id="1">
<density value="4.5e22" units="atom/cm3" />
<nuclide name="U-235" xs="71c" ao="1.0" />
</material>
</materials>

View file

@ -1,2 +0,0 @@
k-combined:
1.092376E+00 1.759788E-02

View file

@ -1,16 +0,0 @@
<?xml version="1.0"?>
<settings>
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
</eigenvalue>
<source>
<space type="box">
<parameters>-4 -4 -4 4 4 4</parameters>
</space>
</source>
</settings>

View file

@ -1,11 +0,0 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness
if __name__ == '__main__':
harness = TestHarness('statepoint.10.*')
harness.main()

View file

@ -1,8 +0,0 @@
<?xml version="1.0"?>
<geometry>
<!-- Sphere with radius 10 -->
<surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
<cell id="1" material="1" region="-1" />
</geometry>

View file

@ -1,9 +0,0 @@
<?xml version="1.0"?>
<materials>
<material id="1">
<density value="12.3e3" units="kg/m3" />
<nuclide name="U-235" xs="71c" ao="1.0" />
</material>
</materials>

View file

@ -1,2 +0,0 @@
k-combined:
7.994522E-01 1.065745E-02

View file

@ -1,11 +0,0 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness
if __name__ == '__main__':
harness = TestHarness('statepoint.10.*')
harness.main()

View file

@ -1,8 +0,0 @@
<?xml version="1.0"?>
<geometry>
<!-- Sphere with radius 10 -->
<surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
<cell id="1" material="1" region="-1" />
</geometry>

View file

@ -1,11 +0,0 @@
<?xml version="1.0"?>
<materials>
<material id="1">
<density units="sum" />
<nuclide name="U-235" xs="71c" ao="0.3e-2" />
<nuclide name="U-238" xs="71c" ao="0.5e-1" />
<nuclide name="H-1" xs="71c" ao="0.1e-2" />
</material>
</materials>

View file

@ -1,2 +0,0 @@
k-combined:
3.231215E-01 6.421320E-03

View file

@ -1,16 +0,0 @@
<?xml version="1.0"?>
<settings>
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
</eigenvalue>
<source>
<space type="box">
<parameters>-4 -4 -4 4 4 4</parameters>
</space>
</source>
</settings>

View file

@ -1,11 +0,0 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness
if __name__ == '__main__':
harness = TestHarness('statepoint.10.*')
harness.main()

View file

@ -0,0 +1,5 @@
<?xml version="1.0"?>
<geometry>
<surface id="1" type="sphere" coeffs="0 0 0 100" boundary="vacuum"/>
<cell id="1" material="1" region="-1" />
</geometry>

View file

@ -0,0 +1,11 @@
<?xml version="1.0"?>
<materials>
<default_xs>71c</default_xs>
<material id="1">
<density value="20" units="g/cc" />
<nuclide name="U-233" ao="1.0" />
<nuclide name="H-2" ao="1.0" />
<nuclide name="Na-23" ao="1.0" />
<nuclide name="Ta-181" ao="1.0" />
</material>
</materials>

View file

@ -0,0 +1,2 @@
k-combined:
2.130076E+00 1.938907E-03

View file

@ -1,16 +1,11 @@
<?xml version="1.0"?>
<settings>
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
</eigenvalue>
<source>
<space type="box">
<parameters>-4 -4 -4 4 4 4</parameters>
</space>
<space type="point" parameters="0. 0. 0." />
</source>
</settings>

View file

@ -0,0 +1,30 @@
#!/usr/bin/env python
"""The purpose of this test is to provide coverage of energy distributions that
are not covered in other tests. It has a single material with the following
nuclides:
U-233: Only nuclide that has a Watt fission spectrum
H-2: Only nuclide that has an N-body phase space distribution, in this case for
(n,2n)
Na-23: Has an evaporation spectrum and also has reactions that have multiple
angle-energy distributions, so it provides coverage for both of those
situations.
Ta-181: One of a few nuclides that has reactions with Kalbach-Mann distributions
that use linear-linear interpolation.
"""
import glob
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness
if __name__ == '__main__':
harness = TestHarness('statepoint.10.*')
harness.main()

View file

@ -1 +0,0 @@
57d6fd9cb5180c38efd2729a5dea0708cbd5fd0bf7dcf0c9d5c9cef5d818aeab5a926d03e70dedcf1b60d5740938fb3ba80e6ccdb09c661d159c0893da3bd593

View file

@ -1,76 +0,0 @@
k-combined:
9.903196E-01 4.279617E-02
tally 1:
4.215917E+01
3.561920E+02
4.174788E+01
3.505184E+02
4.603223E+01
4.242918E+02
4.496760E+01
4.075599E+02
4.088099E+01
3.376516E+02
tally 2:
4.157239E+01
3.482158E+02
4.227810E+01
3.613293E+02
4.376107E+01
3.835007E+02
4.644205E+01
4.327195E+02
4.191554E+01
3.522147E+02
tally 3:
4.215917E+01
3.561920E+02
4.174788E+01
3.505184E+02
4.603223E+01
4.242918E+02
4.496402E+01
4.075053E+02
4.088458E+01
3.377000E+02
tally 4:
1.531988E+01
4.816326E+01
9.274393E+00
1.821174E+01
1.595868E+01
5.124238E+01
1.299895E+00
6.417145E-01
1.510024E+01
4.604170E+01
8.533361E+00
1.462765E+01
1.658141E+01
5.595629E+01
1.427417E+00
6.621807E-01
1.683102E+01
5.741400E+01
9.845257E+00
2.028406E+01
1.773179E+01
6.477077E+01
1.536972E+00
6.111079E-01
1.586070E+01
5.360975E+01
9.928220E+00
2.089005E+01
1.737609E+01
6.161847E+01
1.700608E+00
8.439708E-01
1.607027E+01
5.490113E+01
7.569336E+00
1.280955E+01
1.606086E+01
5.308665E+01
9.898901E-01
3.143027E-01

View file

@ -1,59 +0,0 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
import openmc
class FilterAzimuthalTestHarness(PyAPITestHarness):
def _build_inputs(self):
filt1 = openmc.Filter(type='azimuthal',
bins=(-3.1416, -1.8850, -0.6283, 0.6283, 1.8850,
3.1416))
tally1 = openmc.Tally(tally_id=1)
tally1.add_filter(filt1)
tally1.add_score('flux')
tally1.estimator = 'tracklength'
tally2 = openmc.Tally(tally_id=2)
tally2.add_filter(filt1)
tally2.add_score('flux')
tally2.estimator = 'analog'
filt3 = openmc.Filter(type='azimuthal', bins=(5,))
tally3 = openmc.Tally(tally_id=3)
tally3.add_filter(filt3)
tally3.add_score('flux')
tally3.estimator = 'tracklength'
mesh = openmc.Mesh(mesh_id=1)
mesh.lower_left = [-182.07, -182.07]
mesh.upper_right = [182.07, 182.07]
mesh.dimension = [2, 2]
filt_mesh = openmc.Filter(type='mesh', bins=(1,))
tally4 = openmc.Tally(tally_id=4)
tally4.add_filter(filt3)
tally4.add_filter(filt_mesh)
tally4.add_score('flux')
tally4.estimator = 'tracklength'
self._input_set.tallies = openmc.TalliesFile()
self._input_set.tallies.add_tally(tally1)
self._input_set.tallies.add_tally(tally2)
self._input_set.tallies.add_tally(tally3)
self._input_set.tallies.add_tally(tally4)
self._input_set.tallies.add_mesh(mesh)
super(FilterAzimuthalTestHarness, self)._build_inputs()
def _cleanup(self):
super(FilterAzimuthalTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'tallies.xml')
if os.path.exists(f): os.remove(f)
if __name__ == '__main__':
harness = FilterAzimuthalTestHarness('statepoint.10.*', True)
harness.main()

View file

@ -1 +0,0 @@
f8359184c02fbab5dca5368689a84924066ab1fb09cae575588ceddd696d5461db577498df9959365d89fe933e9b338390e44e362c603c6f2aa5bcf4acc14b20

View file

@ -1,11 +0,0 @@
k-combined:
9.903196E-01 4.279617E-02
tally 1:
0.000000E+00
0.000000E+00
1.767552E+01
6.295417E+01
3.863588E+00
3.013300E+00
5.356594E+01
5.839391E+02

View file

@ -1,29 +0,0 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
import openmc
class FilterCellTestHarness(PyAPITestHarness):
def _build_inputs(self):
filt = openmc.Filter(type='cell', bins=(10, 21, 22, 23))
tally = openmc.Tally(tally_id=1)
tally.add_filter(filt)
tally.add_score('total')
self._input_set.tallies = openmc.TalliesFile()
self._input_set.tallies.add_tally(tally)
super(FilterCellTestHarness, self)._build_inputs()
def _cleanup(self):
super(FilterCellTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'tallies.xml')
if os.path.exists(f): os.remove(f)
if __name__ == '__main__':
harness = FilterCellTestHarness('statepoint.10.*', True)
harness.main()

View file

@ -1 +0,0 @@
8ae662f8881ce8cdec550069c6233c2c91e9a10f7200af6892cf6f2d77712ccfa17895dbd2eee02e6daf3d665c6ed84b29e17d89ff519e70c37b36d75a431d53

View file

@ -1,11 +0,0 @@
k-combined:
9.903196E-01 4.279617E-02
tally 1:
0.000000E+00
0.000000E+00
8.921179E+01
1.601939E+03
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00

View file

@ -1,29 +0,0 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
import openmc
class FilterCellbornTestHarness(PyAPITestHarness):
def _build_inputs(self):
filt = openmc.Filter(type='cellborn', bins=(10, 21, 22, 23))
tally = openmc.Tally(tally_id=1)
tally.add_filter(filt)
tally.add_score('total')
self._input_set.tallies = openmc.TalliesFile()
self._input_set.tallies.add_tally(tally)
super(FilterCellbornTestHarness, self)._build_inputs()
def _cleanup(self):
super(FilterCellbornTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'tallies.xml')
if os.path.exists(f): os.remove(f)
if __name__ == '__main__':
harness = FilterCellbornTestHarness('statepoint.10.*', True)
harness.main()

View file

@ -1 +0,0 @@
a7c8ce7ffbc3a7b965d8a3077a4d9132130561afef19047b279b2d23198e248b09664856a092a32394894e19fef7708cebad99b3839d735c4e98ae0c9af58cb7

View file

@ -1,15 +0,0 @@
k-combined:
9.903196E-01 4.279617E-02
tally 1:
8.141852E-04
1.337187E-07
4.849156E-03
4.744020E-06
4.460252E-03
4.015453E-06
1.028479E-02
2.136252E-05
5.002274E-03
5.056965E-06
1.974747E-03
7.882970E-07

View file

@ -1,30 +0,0 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
import openmc
class FilterDelayedgroupTestHarness(PyAPITestHarness):
def _build_inputs(self):
filt = openmc.Filter(type='delayedgroup',
bins=(1, 2, 3, 4, 5, 6))
tally = openmc.Tally(tally_id=1)
tally.add_filter(filt)
tally.add_score('delayed-nu-fission')
self._input_set.tallies = openmc.TalliesFile()
self._input_set.tallies.add_tally(tally)
super(FilterDelayedgroupTestHarness, self)._build_inputs()
def _cleanup(self):
super(FilterDelayedgroupTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'tallies.xml')
if os.path.exists(f): os.remove(f)
if __name__ == '__main__':
harness = FilterDelayedgroupTestHarness('statepoint.10.*', True)
harness.main()

View file

@ -1 +0,0 @@
51d3e2c43f36712a7b26c5fa26e0e2ca6fb9af205af04f0f8cd44c6b100e36382417c2c63d711e4677ce3c1958d15072727d5fd32424a3f6eb08d1f3b1c7db5a

View file

@ -1,11 +0,0 @@
k-combined:
9.903196E-01 4.279617E-02
tally 1:
2.844008E+01
1.619630E+02
4.425619E+01
3.938244E+02
5.527425E+01
6.120383E+02
9.799897E+00
1.957877E+01

View file

@ -1,30 +0,0 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
import openmc
class FilterEnergyTestHarness(PyAPITestHarness):
def _build_inputs(self):
filt = openmc.Filter(type='energy',
bins=(0.0, 0.253e-6, 1.0e-3, 1.0, 20.0))
tally = openmc.Tally(tally_id=1)
tally.add_filter(filt)
tally.add_score('total')
self._input_set.tallies = openmc.TalliesFile()
self._input_set.tallies.add_tally(tally)
super(FilterEnergyTestHarness, self)._build_inputs()
def _cleanup(self):
super(FilterEnergyTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'tallies.xml')
if os.path.exists(f): os.remove(f)
if __name__ == '__main__':
harness = FilterEnergyTestHarness('statepoint.10.*', True)
harness.main()

View file

@ -1 +0,0 @@
f0810606c5f947a9fe03bcfc87de3883ce46f59d8603e02ed30f853ebf301b2dc6bdcd109889801ada9e6e0b7be4932efeca97d4beea875af8c8e3ecb7511444

View file

@ -1,11 +0,0 @@
k-combined:
9.903196E-01 4.279617E-02
tally 1:
2.842000E+01
1.620214E+02
4.361000E+01
3.810139E+02
5.297000E+01
5.616595E+02
6.530000E+00
8.828900E+00

View file

@ -1,30 +0,0 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
import openmc
class FilterEnergyoutTestHarness(PyAPITestHarness):
def _build_inputs(self):
filt = openmc.Filter(type='energyout',
bins=(0.0, 0.253e-6, 1.0e-3, 1.0, 20.0))
tally = openmc.Tally(tally_id=1)
tally.add_filter(filt)
tally.add_score('scatter')
self._input_set.tallies = openmc.TalliesFile()
self._input_set.tallies.add_tally(tally)
super(FilterEnergyoutTestHarness, self)._build_inputs()
def _cleanup(self):
super(FilterEnergyoutTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'tallies.xml')
if os.path.exists(f): os.remove(f)
if __name__ == '__main__':
harness = FilterEnergyoutTestHarness('statepoint.10.*', True)
harness.main()

View file

@ -1 +0,0 @@
c4d4334d44956d6dc9abe854a5e9403d7f8a87ffb04a15a3d128e8d18eb4111f46ca277b751e1b0e836d69527502f9abba115a4b2fc64c38da63a9d57968d860

View file

@ -1,67 +0,0 @@
k-combined:
9.903196E-01 4.279617E-02
tally 1:
2.576000E+01
1.331666E+02
0.000000E+00
0.000000E+00
7.000000E-02
1.300000E-03
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.050675E+00
2.274991E-01
0.000000E+00
0.000000E+00
2.070821E+00
8.886068E-01
2.660000E+00
1.422000E+00
0.000000E+00
0.000000E+00
3.897000E+01
3.042635E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
4.352932E-01
4.705717E-02
0.000000E+00
0.000000E+00
1.018668E+00
2.090017E-01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
4.570000E+00
4.182700E+00
0.000000E+00
0.000000E+00
4.968000E+01
4.940534E+02
6.537406E-02
1.230788E-03
0.000000E+00
0.000000E+00
8.678070E-02
2.482037E-03
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
3.290000E+00
2.178900E+00
1.610879E-01
5.883677E-03
6.530000E+00
8.828900E+00
3.151783E-01
2.052521E-02

View file

@ -1,34 +0,0 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
import openmc
class FilterGroupTransferTestHarness(PyAPITestHarness):
def _build_inputs(self):
filt1 = openmc.Filter(type='energy',
bins=(0.0, 0.253e-6, 1.0e-3, 1.0, 20.0))
filt2 = openmc.Filter(type='energyout',
bins=(0.0, 0.253e-6, 1.0e-3, 1.0, 20.0))
tally = openmc.Tally(tally_id=1)
tally.add_filter(filt1)
tally.add_filter(filt2)
tally.add_score('scatter')
tally.add_score('nu-fission')
self._input_set.tallies = openmc.TalliesFile()
self._input_set.tallies.add_tally(tally)
super(FilterGroupTransferTestHarness, self)._build_inputs()
def _cleanup(self):
super(FilterGroupTransferTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'tallies.xml')
if os.path.exists(f): os.remove(f)
if __name__ == '__main__':
harness = FilterGroupTransferTestHarness('statepoint.10.*', True)
harness.main()

View file

@ -1 +0,0 @@
7689b2c88391128377b7f9bfcda347a42f77d69d194186629fa965ecd3fc51be0bfd1ac92fb9d7551128d8b6ed5241ead4fb94b27ae29d80230863e78fbbcb68

View file

@ -1,11 +0,0 @@
k-combined:
9.903196E-01 4.279617E-02
tally 1:
2.868239E+01
1.648549E+02
6.779424E+00
9.202676E+00
6.446222E+01
8.387204E+02
3.367496E+01
2.349072E+02

View file

@ -1,29 +0,0 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
import openmc
class FilterMaterialTestHarness(PyAPITestHarness):
def _build_inputs(self):
filt = openmc.Filter(type='material', bins=(1, 2, 3, 4))
tally = openmc.Tally(tally_id=1)
tally.add_filter(filt)
tally.add_score('total')
self._input_set.tallies = openmc.TalliesFile()
self._input_set.tallies.add_tally(tally)
super(FilterMaterialTestHarness, self)._build_inputs()
def _cleanup(self):
super(FilterMaterialTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'tallies.xml')
if os.path.exists(f): os.remove(f)
if __name__ == '__main__':
harness = FilterMaterialTestHarness('statepoint.10.*', True)
harness.main()

View file

@ -1 +0,0 @@
ecc649936e2cc364b079944f47e18fb81ec7290017b4bd5837e5aa1e24e1146df77897f44c7c2a88500e3f525566b51777cd9b84ec6a636f5883e411e4c1f75c

View file

@ -1,121 +0,0 @@
k-combined:
9.903196E-01 4.279617E-02
tally 1:
1.241000E+01
3.088870E+01
1.241000E+01
3.088870E+01
1.364000E+01
3.727140E+01
1.364000E+01
3.727140E+01
3.251000E+01
2.118597E+02
3.251000E+01
2.118597E+02
7.297000E+01
1.066904E+03
7.297000E+01
1.066904E+03
tally 2:
9.880000E+00
1.964520E+01
9.880000E+00
1.964520E+01
1.022000E+01
2.099620E+01
1.022000E+01
2.099620E+01
1.479000E+01
4.397670E+01
1.479000E+01
4.397670E+01
3.470000E+01
2.412094E+02
3.470000E+01
2.412094E+02
6.194000E+01
7.687326E+02
6.194000E+01
7.687326E+02
tally 3:
3.560000E+00
2.681800E+00
3.560000E+00
2.681800E+00
1.930000E+00
7.915000E-01
1.930000E+00
7.915000E-01
3.870000E+00
3.109100E+00
3.870000E+00
3.109100E+00
3.500000E-01
3.630000E-02
3.500000E-01
3.630000E-02
3.680000E+00
2.840200E+00
3.680000E+00
2.840200E+00
2.050000E+00
8.735000E-01
2.050000E+00
8.735000E-01
3.910000E+00
3.085100E+00
3.910000E+00
3.085100E+00
3.900000E-01
3.610000E-02
3.900000E-01
3.610000E-02
5.130000E+00
5.422100E+00
5.130000E+00
5.422100E+00
3.100000E+00
1.959200E+00
3.100000E+00
1.959200E+00
5.840000E+00
6.914600E+00
5.840000E+00
6.914600E+00
5.400000E-01
8.980000E-02
5.400000E-01
8.980000E-02
1.215000E+01
3.061010E+01
1.215000E+01
3.061010E+01
7.220000E+00
1.081680E+01
7.220000E+00
1.081680E+01
1.355000E+01
3.699090E+01
1.355000E+01
3.699090E+01
1.360000E+00
5.098000E-01
1.360000E+00
5.098000E-01
2.199000E+01
9.837430E+01
2.199000E+01
9.837430E+01
1.243000E+01
3.167470E+01
1.243000E+01
3.167470E+01
2.451000E+01
1.233915E+02
2.451000E+01
1.233915E+02
2.460000E+00
1.687000E+00
2.460000E+00
1.687000E+00

View file

@ -1,53 +0,0 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
import openmc
class FilterMuTestHarness(PyAPITestHarness):
def _build_inputs(self):
filt1 = openmc.Filter(type='mu',
bins=(-1.0, -0.5, 0.0, 0.5, 1.0))
tally1 = openmc.Tally(tally_id=1)
tally1.add_filter(filt1)
tally1.add_score('scatter')
tally1.add_score('nu-scatter')
filt2 = openmc.Filter(type='mu', bins=(5,))
tally2 = openmc.Tally(tally_id=2)
tally2.add_filter(filt2)
tally2.add_score('scatter')
tally2.add_score('nu-scatter')
mesh = openmc.Mesh(mesh_id=1)
mesh.lower_left = [-182.07, -182.07]
mesh.upper_right = [182.07, 182.07]
mesh.dimension = [2, 2]
filt_mesh = openmc.Filter(type='mesh', bins=(1,))
tally3 = openmc.Tally(tally_id=3)
tally3.add_filter(filt2)
tally3.add_filter(filt_mesh)
tally3.add_score('scatter')
tally3.add_score('nu-scatter')
self._input_set.tallies = openmc.TalliesFile()
self._input_set.tallies.add_tally(tally1)
self._input_set.tallies.add_tally(tally2)
self._input_set.tallies.add_tally(tally3)
self._input_set.tallies.add_mesh(mesh)
super(FilterMuTestHarness, self)._build_inputs()
def _cleanup(self):
super(FilterMuTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'tallies.xml')
if os.path.exists(f): os.remove(f)
if __name__ == '__main__':
harness = FilterMuTestHarness('statepoint.10.*', True)
harness.main()

View file

@ -1 +0,0 @@
301824991a022884215609f39797a61933faf7ccacf81ad6bb883af08857563e8bd74ab946fc4fd072860168d77f76d0c76d1467375158072dce431fc6a1c449

View file

@ -1,76 +0,0 @@
k-combined:
9.903196E-01 4.279617E-02
tally 1:
2.127061E+01
9.220793E+01
5.602776E+01
6.373945E+02
6.367492E+01
8.138443E+02
5.529942E+01
6.140264E+02
1.951517E+01
7.668661E+01
tally 2:
2.075936E+01
8.757254E+01
5.524881E+01
6.153139E+02
6.475252E+01
8.402281E+02
5.446664E+01
5.961174E+02
2.074180E+01
8.681580E+01
tally 3:
2.128073E+01
9.230382E+01
5.601764E+01
6.371703E+02
6.367492E+01
8.138443E+02
5.529942E+01
6.140264E+02
1.951517E+01
7.668661E+01
tally 4:
8.088647E+00
1.396899E+01
3.960907E+00
3.249150E+00
8.430714E+00
1.435355E+01
7.192159E-01
1.641710E-01
1.974619E+01
8.105078E+01
1.212452E+01
3.016420E+01
2.228348E+01
1.050847E+02
1.748809E+00
9.501796E-01
2.257423E+01
1.038902E+02
1.351331E+01
3.969787E+01
2.507638E+01
1.283664E+02
2.193118E+00
1.424580E+00
2.192232E+01
9.859711E+01
1.096779E+01
2.506373E+01
2.074138E+01
8.670015E+01
1.469145E+00
8.072204E-01
6.850719E+00
9.425536E+00
4.584038E+00
4.399762E+00
7.176883E+00
1.090693E+01
8.244944E-01
1.794291E-01

View file

@ -1,59 +0,0 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
import openmc
class FilterPolarTestHarness(PyAPITestHarness):
def _build_inputs(self):
filt1 = openmc.Filter(type='polar',
bins=(0.0, 0.6283, 1.2566, 1.8850, 2.5132,
3.1416))
tally1 = openmc.Tally(tally_id=1)
tally1.add_filter(filt1)
tally1.add_score('flux')
tally1.estimator = 'tracklength'
tally2 = openmc.Tally(tally_id=2)
tally2.add_filter(filt1)
tally2.add_score('flux')
tally2.estimator = 'analog'
filt3 = openmc.Filter(type='polar', bins=(5,))
tally3 = openmc.Tally(tally_id=3)
tally3.add_filter(filt3)
tally3.add_score('flux')
tally3.estimator = 'tracklength'
mesh = openmc.Mesh(mesh_id=1)
mesh.lower_left = [-182.07, -182.07]
mesh.upper_right = [182.07, 182.07]
mesh.dimension = [2, 2]
filt_mesh = openmc.Filter(type='mesh', bins=(1,))
tally4 = openmc.Tally(tally_id=4)
tally4.add_filter(filt3)
tally4.add_filter(filt_mesh)
tally4.add_score('flux')
tally4.estimator = 'tracklength'
self._input_set.tallies = openmc.TalliesFile()
self._input_set.tallies.add_tally(tally1)
self._input_set.tallies.add_tally(tally2)
self._input_set.tallies.add_tally(tally3)
self._input_set.tallies.add_tally(tally4)
self._input_set.tallies.add_mesh(mesh)
super(FilterPolarTestHarness, self)._build_inputs()
def _cleanup(self):
super(FilterPolarTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'tallies.xml')
if os.path.exists(f): os.remove(f)
if __name__ == '__main__':
harness = FilterPolarTestHarness('statepoint.10.*', True)
harness.main()

View file

@ -1 +0,0 @@
164804414f48a818c93e197f2901ce6ae375d88071a03e89c920dbc4462e7a2c8d2c85acf6560fcd6eb3d7c0c53d3b426ab1cc4b7721266fe8adec3e7231149e

View file

@ -1,11 +0,0 @@
k-combined:
9.903196E-01 4.279617E-02
tally 1:
7.510505E+01
1.143811E+03
8.792943E+00
1.575416E+01
4.214462E+01
3.642975E+02
4.335157E+00
3.864423E+00

View file

@ -1,29 +0,0 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
import openmc
class FilterUniverseTestHarness(PyAPITestHarness):
def _build_inputs(self):
filt = openmc.Filter(type='universe', bins=(1, 2, 3, 4))
tally = openmc.Tally(tally_id=1)
tally.add_filter(filt)
tally.add_score('total')
self._input_set.tallies = openmc.TalliesFile()
self._input_set.tallies.add_tally(tally)
super(FilterUniverseTestHarness, self)._build_inputs()
def _cleanup(self):
super(FilterUniverseTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'tallies.xml')
if os.path.exists(f): os.remove(f)
if __name__ == '__main__':
harness = FilterUniverseTestHarness('statepoint.10.*', True)
harness.main()

View file

@ -1,181 +0,0 @@
<?xml version="1.0"?>
<geometry>
<surface id="1" type="z-cylinder" coeffs="0. 0. 0.41" />
<surface id="2" type="z-cylinder" coeffs="0. 0. 0.475" />
<surface id="3" type="z-cylinder" coeffs="0. 0. 0.56" />
<surface id="4" type="z-cylinder" coeffs="0. 0. 0.62" />
<surface id="5" type="z-cylinder" coeffs="0. 0. 187.6" />
<surface id="6" type="z-cylinder" coeffs="0. 0. 209.0" />
<surface id="7" type="z-cylinder" coeffs="0. 0. 229.0" />
<surface id="8" type="z-cylinder" coeffs="0. 0. 249.0" boundary="vacuum" />
<surface id="31" type="z-plane" coeffs="-229.0" boundary="vacuum" />
<surface id="32" type="z-plane" coeffs="-199.0" />
<surface id="33" type="z-plane" coeffs="-193.0" />
<surface id="34" type="z-plane" coeffs="-183.0" />
<surface id="35" type="z-plane" coeffs="0.0" />
<surface id="36" type="z-plane" coeffs="183.0" />
<surface id="37" type="z-plane" coeffs="203.0" />
<surface id="38" type="z-plane" coeffs="215.0" />
<surface id="39" type="z-plane" coeffs="223.0" boundary="vacuum" />
<!-- All geometry on base universe -->
<cell id="1" fill="200" region=" -6 34 -35" /> <!-- Lower core -->
<cell id="2" fill="201" region=" -6 35 -36" /> <!-- Upper core -->
<cell id="3" material="8" region=" -7 31 -32" /> <!-- Lower core plate region -->
<cell id="4" material="9" region=" -5 32 -33" /> <!-- Bottom nozzle region -->
<cell id="5" material="12" region=" -5 33 -34" /> <!-- Bottom FA region -->
<cell id="6" material="11" region=" -5 36 -37" /> <!-- Top FA region -->
<cell id="7" material="10" region=" -5 37 -38" /> <!-- Top nozzle region -->
<cell id="8" material="7" region=" -7 38 -39" /> <!-- Upper plate region -->
<cell id="9" material="4" region="6 -7 32 -38" /> <!-- Downcomer -->
<cell id="10" material="5" region="7 -8 31 -39" /> <!-- RPV -->
<cell id="11" material="6" region="5 -6 32 -34" /> <!-- Bottom of radial reflector -->
<cell id="12" material="7" region="5 -6 36 -38" /> <!-- Top of radial reflector -->
<!-- Fuel pin, cladding, cold water -->
<cell id="21" universe="1" material="1" region="-1" />
<cell id="22" universe="1" material="2" region="1 -2" />
<cell id="23" universe="1" material="3" region="2" />
<!-- Instrumentation guide tube -->
<cell id="24" universe="2" material="3" region="-3" />
<cell id="25" universe="2" material="2" region="3 -4" />
<cell id="26" universe="2" material="3" region="4" />
<!-- Fuel pin, cladding, hot water -->
<cell id="27" universe="3" material="1" region="-1" />
<cell id="28" universe="3" material="2" region="1 -2" />
<cell id="29" universe="3" material="4" region="2" />
<!-- Instrumentation guide tube -->
<cell id="30" universe="4" material="4" region="-3" />
<cell id="31" universe="4" material="2" region="3 -4" />
<cell id="32" universe="4" material="4" region="4" />
<!-- cell for water assembly (cold) -->
<cell id="50" universe="5" material="4" region="34 -35" />
<!-- containing cell for fuel assembly -->
<cell id="60" universe="6" fill="100" region="34 -35" />
<!-- cell for water assembly (hot) -->
<cell id="70" universe="7" material="3" region="35 -36" />
<!-- containing cell for fuel assembly -->
<cell id="80" universe="8" fill="101" region="35 -36" />
<!-- Fuel Assembly (Lower Half) -->
<lattice id="100">
<dimension>17 17</dimension>
<lower_left>-10.71 -10.71</lower_left>
<pitch>1.26 1.26</pitch>
<universes>
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1
1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1
1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
</universes>
</lattice>
<!-- Fuel Assembly (Upper Half) -->
<lattice id="101">
<dimension>17 17</dimension>
<lower_left>-10.71 -10.71</lower_left>
<pitch>1.26 1.26</pitch>
<universes>
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3
3 3 3 4 3 3 3 3 3 3 3 3 3 4 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 4 3 3 3 3 3 3 3 3 3 4 3 3 3
3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
</universes>
</lattice>
<!-- Core Lattice (Lower Half) -->
<lattice id="200">
<dimension>21 21</dimension>
<lower_left>-224.91 -224.91</lower_left>
<pitch>21.42 21.42</pitch>
<universes>
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
5 5 5 5 5 5 5 6 6 6 6 6 6 6 5 5 5 5 5 5 5
5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5
5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5
5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5
5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5
5 5 5 5 5 5 5 6 6 6 6 6 6 6 5 5 5 5 5 5 5
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
</universes>
</lattice>
<!-- Core Lattice (Upper Half) -->
<lattice id="201">
<dimension>21 21</dimension>
<lower_left>-224.91 -224.91</lower_left>
<pitch>21.42 21.42</pitch>
<universes>
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
7 7 7 7 7 7 7 8 8 8 8 8 8 8 7 7 7 7 7 7 7
7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7
7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7
7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7
7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7
7 7 7 7 7 7 7 8 8 8 8 8 8 8 7 7 7 7 7 7 7
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
</universes>
</lattice>
</geometry>

View file

@ -1,272 +0,0 @@
<?xml version="1.0"?>
<materials>
<default_xs>71c</default_xs>
<!-- Fuel composition -->
<material id="1">
<density value="10.062" units="g/cm3" />
<nuclide name="U-234" ao="4.9476e-6" />
<nuclide name="U-235" ao="4.8218e-4" />
<nuclide name="U-236" ao="9.0402e-5" />
<nuclide name="U-238" ao="2.1504e-2" />
<nuclide name="Np-237" ao="7.3733e-6" />
<nuclide name="Pu-238" ao="1.5148e-6" />
<nuclide name="Pu-239" ao="1.3955e-4" />
<nuclide name="Pu-240" ao="3.4405e-5" />
<nuclide name="Pu-241" ao="2.1439e-5" />
<nuclide name="Pu-242" ao="3.7422e-6" />
<nuclide name="Am-241" ao="4.5041e-7" />
<nuclide name="Am-242m" ao="9.2301e-9" />
<nuclide name="Am-243" ao="4.7878e-7" />
<nuclide name="Cm-242" ao="1.0485e-7" />
<nuclide name="Cm-243" ao="1.4268e-9" />
<nuclide name="Cm-244" ao="8.8756e-8" />
<nuclide name="Cm-245" ao="3.5285e-9" />
<nuclide name="Mo-95" ao="2.6497e-5" />
<nuclide name="Tc-99" ao="3.2772e-5" />
<nuclide name="Ru-101" ao="3.0742e-5" />
<nuclide name="Ru-103" ao="2.3505e-6" />
<nuclide name="Ag-109" ao="2.0009e-6" />
<nuclide name="Xe-135" ao="1.0801e-8" />
<nuclide name="Cs-133" ao="3.4612e-5" />
<nuclide name="Nd-143" ao="2.6078e-5" />
<nuclide name="Nd-145" ao="1.9898e-5" />
<nuclide name="Sm-147" ao="1.6128e-6" />
<nuclide name="Sm-149" ao="1.1627e-7" />
<nuclide name="Sm-150" ao="7.1727e-6" />
<nuclide name="Sm-151" ao="5.4947e-7" />
<nuclide name="Sm-152" ao="3.0221e-6" />
<nuclide name="Eu-153" ao="2.6209e-6" />
<nuclide name="Gd-155" ao="1.5369e-9" />
<nuclide name="O-16" ao="4.5737e-2" />
</material>
<!-- Cladding composition -->
<material id="2">
<density value="5.77" units="g/cm3" />
<nuclide name="Zr-90" ao="0.5145" />
<nuclide name="Zr-91" ao="0.1122" />
<nuclide name="Zr-92" ao="0.1715" />
<nuclide name="Zr-94" ao="0.1738" />
<nuclide name="Zr-96" ao="0.0280" />
</material>
<!-- Cold borated water -->
<material id="3">
<density value="0.07416" units="atom/b-cm" />
<nuclide name="H-1" ao="2.0" />
<nuclide name="O-16" ao="1.0" />
<nuclide name="B-10" ao="6.490e-4" />
<nuclide name="B-11" ao="2.689e-3" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Hot borated water -->
<material id="4">
<density value="0.06614" units="atom/b-cm" />
<nuclide name="H-1" ao="2.0" />
<nuclide name="O-16" ao="1.0" />
<nuclide name="B-10" ao="6.490e-4" />
<nuclide name="B-11" ao="2.689e-3" />
<sab name="HH2O" xs="71t" />
</material>
<!-- RPV Composition -->
<material id="5">
<density value="7.9" units="g/cm3" />
<nuclide name="Fe-54" wo="0.05437098" />
<nuclide name="Fe-56" wo="0.88500663" />
<nuclide name="Fe-57" wo="0.0208008" />
<nuclide name="Fe-58" wo="0.00282159" />
<nuclide name="Ni-58" wo="0.0067198" />
<nuclide name="Ni-60" wo="0.0026776" />
<nuclide name="Ni-61" wo="0.0001183" />
<nuclide name="Ni-62" wo="0.0003835" />
<nuclide name="Ni-64" wo="0.0001008" />
<nuclide name="Mn-55" wo="0.01" />
<nuclide name="Mo-92" wo="0.000849" />
<nuclide name="Mo-94" wo="0.0005418" />
<nuclide name="Mo-95" wo="0.0009438" />
<nuclide name="Mo-96" wo="0.0010002" />
<nuclide name="Mo-97" wo="0.0005796" />
<nuclide name="Mo-98" wo="0.0014814" />
<nuclide name="Mo-100" wo="0.0006042" />
<nuclide name="Si-28" wo="0.00367464" />
<nuclide name="Si-29" wo="0.00019336" />
<nuclide name="Si-30" wo="0.000132" />
<nuclide name="Cr-50" wo="0.00010435" />
<nuclide name="Cr-52" wo="0.002092475" />
<nuclide name="Cr-53" wo="0.00024185" />
<nuclide name="Cr-54" wo="6.1325e-05" />
<nuclide name="C-Nat" wo="0.0025" />
<nuclide name="Cu-63" wo="0.0013696" />
<nuclide name="Cu-65" wo="0.0006304" />
</material>
<!-- Lower radial reflector -->
<material id="6">
<density value="4.32" units="g/cm3" />
<nuclide name="H-1" wo="0.0095661" />
<nuclide name="O-16" wo="0.0759107" />
<nuclide name="B-10" wo="3.08409e-5" />
<nuclide name="B-11" wo="1.40499e-4" />
<nuclide name="Fe-54" wo="0.035620772088" />
<nuclide name="Fe-56" wo="0.579805982228" />
<nuclide name="Fe-57" wo="0.01362750048" />
<nuclide name="Fe-58" wo="0.001848545204" />
<nuclide name="Ni-58" wo="0.055298376566" />
<nuclide name="Ni-60" wo="0.022034425592" />
<nuclide name="Ni-61" wo="0.000973510811" />
<nuclide name="Ni-62" wo="0.003155886695" />
<nuclide name="Ni-64" wo="0.000829500336" />
<nuclide name="Mn-55" wo="0.0182870" />
<nuclide name="Si-28" wo="0.00839976771" />
<nuclide name="Si-29" wo="0.00044199679" />
<nuclide name="Si-30" wo="0.0003017355" />
<nuclide name="Cr-50" wo="0.007251360806" />
<nuclide name="Cr-52" wo="0.145407678031" />
<nuclide name="Cr-53" wo="0.016806340306" />
<nuclide name="Cr-54" wo="0.004261520857" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Upper radial reflector / Top plate region -->
<material id="7">
<density value="4.28" units="g/cm3" />
<nuclide name="H-1" wo="0.0086117" />
<nuclide name="O-16" wo="0.0683369" />
<nuclide name="B-10" wo="2.77638e-5" />
<nuclide name="B-11" wo="1.26481e-4" />
<nuclide name="Fe-54" wo="0.035953677186" />
<nuclide name="Fe-56" wo="0.585224740891" />
<nuclide name="Fe-57" wo="0.01375486056" />
<nuclide name="Fe-58" wo="0.001865821363" />
<nuclide name="Ni-58" wo="0.055815129186" />
<nuclide name="Ni-60" wo="0.022240333032" />
<nuclide name="Ni-61" wo="0.000982608081" />
<nuclide name="Ni-62" wo="0.003185377845" />
<nuclide name="Ni-64" wo="0.000837251856" />
<nuclide name="Mn-55" wo="0.0184579" />
<nuclide name="Si-28" wo="0.00847831314" />
<nuclide name="Si-29" wo="0.00044612986" />
<nuclide name="Si-30" wo="0.000304557" />
<nuclide name="Cr-50" wo="0.00731912987" />
<nuclide name="Cr-52" wo="0.146766614995" />
<nuclide name="Cr-53" wo="0.01696340737" />
<nuclide name="Cr-54" wo="0.004301347765" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Bottom plate region -->
<material id="8">
<density value="7.184" units="g/cm3" />
<nuclide name="H-1" wo="0.0011505" />
<nuclide name="O-16" wo="0.0091296" />
<nuclide name="B-10" wo="3.70915e-6" />
<nuclide name="B-11" wo="1.68974e-5" />
<nuclide name="Fe-54" wo="0.03855611055" />
<nuclide name="Fe-56" wo="0.627585036425" />
<nuclide name="Fe-57" wo="0.014750478" />
<nuclide name="Fe-58" wo="0.002000875025" />
<nuclide name="Ni-58" wo="0.059855207342" />
<nuclide name="Ni-60" wo="0.023850159704" />
<nuclide name="Ni-61" wo="0.001053732407" />
<nuclide name="Ni-62" wo="0.003415945715" />
<nuclide name="Ni-64" wo="0.000897854832" />
<nuclide name="Mn-55" wo="0.0197940" />
<nuclide name="Si-28" wo="0.00909197802" />
<nuclide name="Si-29" wo="0.00047842098" />
<nuclide name="Si-30" wo="0.000326601" />
<nuclide name="Cr-50" wo="0.007848910646" />
<nuclide name="Cr-52" wo="0.157390026871" />
<nuclide name="Cr-53" wo="0.018191270146" />
<nuclide name="Cr-54" wo="0.004612692337" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Bottom nozzle region -->
<material id="9">
<density value="2.53" units="g/cm3" />
<nuclide name="H-1" wo="0.0245014" />
<nuclide name="O-16" wo="0.1944274" />
<nuclide name="B-10" wo="7.89917e-5" />
<nuclide name="B-11" wo="3.59854e-4" />
<nuclide name="Fe-54" wo="0.030411411144" />
<nuclide name="Fe-56" wo="0.495012237964" />
<nuclide name="Fe-57" wo="0.01163454624" />
<nuclide name="Fe-58" wo="0.001578204652" />
<nuclide name="Ni-58" wo="0.047211231662" />
<nuclide name="Ni-60" wo="0.018811987544" />
<nuclide name="Ni-61" wo="0.000831139127" />
<nuclide name="Ni-62" wo="0.002694352115" />
<nuclide name="Ni-64" wo="0.000708189552" />
<nuclide name="Mn-55" wo="0.0156126" />
<nuclide name="Si-28" wo="0.007171335558" />
<nuclide name="Si-29" wo="0.000377356542" />
<nuclide name="Si-30" wo="0.0002576079" />
<nuclide name="Cr-50" wo="0.006190885148" />
<nuclide name="Cr-52" wo="0.124142524198" />
<nuclide name="Cr-53" wo="0.014348496148" />
<nuclide name="Cr-54" wo="0.003638294506" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Top nozzle region -->
<material id="10">
<density value="1.746" units="g/cm3" />
<nuclide name="H-1" wo="0.0358870" />
<nuclide name="O-16" wo="0.2847761" />
<nuclide name="B-10" wo="1.15699e-4" />
<nuclide name="B-11" wo="5.27075e-4" />
<nuclide name="Fe-54" wo="0.02644016154" />
<nuclide name="Fe-56" wo="0.43037146399" />
<nuclide name="Fe-57" wo="0.0101152584" />
<nuclide name="Fe-58" wo="0.00137211607" />
<nuclide name="Ni-58" wo="0.04104621835" />
<nuclide name="Ni-60" wo="0.0163554502" />
<nuclide name="Ni-61" wo="0.000722605975" />
<nuclide name="Ni-62" wo="0.002342513875" />
<nuclide name="Ni-64" wo="0.0006157116" />
<nuclide name="Mn-55" wo="0.0135739" />
<nuclide name="Si-28" wo="0.006234853554" />
<nuclide name="Si-29" wo="0.000328078746" />
<nuclide name="Si-30" wo="0.0002239677" />
<nuclide name="Cr-50" wo="0.005382452306" />
<nuclide name="Cr-52" wo="0.107931450781" />
<nuclide name="Cr-53" wo="0.012474806806" />
<nuclide name="Cr-54" wo="0.003163190107" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Top of Fuel Assemblies -->
<material id="11">
<density value="3.044" units="g/cm3" />
<nuclide name="H-1" wo="0.0162913" />
<nuclide name="O-16" wo="0.1292776" />
<nuclide name="B-10" wo="5.25228e-5" />
<nuclide name="B-11" wo="2.39272e-4" />
<nuclide name="Zr-90" wo="0.43313403903" />
<nuclide name="Zr-91" wo="0.09549277374" />
<nuclide name="Zr-92" wo="0.14759527104" />
<nuclide name="Zr-94" wo="0.15280552077" />
<nuclide name="Zr-96" wo="0.02511169542" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Bottom of Fuel Assemblies -->
<material id="12">
<density value="1.762" units="g/cm3" />
<nuclide name="H-1" wo="0.0292856" />
<nuclide name="O-16" wo="0.2323919" />
<nuclide name="B-10" wo="9.44159e-5" />
<nuclide name="B-11" wo="4.30120e-4" />
<nuclide name="Zr-90" wo="0.3741373658" />
<nuclide name="Zr-91" wo="0.0824858164" />
<nuclide name="Zr-92" wo="0.1274914944" />
<nuclide name="Zr-94" wo="0.1319920622" />
<nuclide name="Zr-96" wo="0.0216912612" />
<sab name="HH2O" xs="71t" />
</material>
</materials>

View file

@ -1,101 +0,0 @@
k-combined:
0.000000E+00 0.000000E+00
tally 1:
2.247257E+01
1.683779E+02
1.014000E+01
3.427342E+01
8.628000E+00
2.481430E+01
8.632000E+00
2.483728E+01
5.102293E-01
8.710841E-02
5.087118E-01
8.657086E-02
9.212024E+00
2.829472E+01
8.628000E+00
2.481430E+01
1.512000E+00
7.620560E-01
1.816851E+00
1.102658E+00
1.337996E+02
5.985519E+03
2.247257E+01
1.683779E+02
1.512960E-01
2.623972E-02
-3.775020E-01
1.055377E-01
1.916133E-01
4.680798E-02
2.754367E-02
3.320008E-04
2.028374E-02
1.319357E-02
8.974271E-03
1.681081E-03
1.658978E-01
1.520448E-02
2.878360E-01
5.645480E-02
1.014000E+01
3.427342E+01
4.798897E-02
1.551226E-03
-1.818770E-01
1.492633E-02
6.340651E-02
9.011305E-03
3.395308E-02
4.612818E-04
2.640250E-02
6.434787E-04
-8.242639E-03
9.516540E-04
8.378601E-02
2.645988E-03
9.567484E-02
7.262477E-03
8.628000E+00
2.481430E+01
4.712248E-02
1.140942E-03
-6.431930E-02
4.290580E-03
9.251642E-02
8.134201E-03
1.020119E-04
1.154184E-04
2.994164E-02
3.079076E-04
2.128844E-02
2.046549E-04
-1.637972E-02
1.459209E-04
4.629047E-02
7.823267E-04
8.632000E+00
2.483728E+01
4.651997E-02
1.133839E-03
-6.416955E-02
4.279418E-03
9.280565E-02
8.095106E-03
-2.078094E-04
1.151292E-04
3.005568E-02
3.104764E-04
2.199519E-02
2.179172E-04
-1.660645E-02
1.451345E-04
4.607553E-02
7.673412E-04
1.014000E+01
3.427342E+01
7.652723E-03
3.578992E-05

View file

@ -1,19 +0,0 @@
<?xml version="1.0"?>
<settings>
<eigenvalue>
<batches>5</batches>
<inactive>2</inactive>
<particles>500</particles>
</eigenvalue>
<source>
<space type="box">
<parameters>
-160 -160 -183
160 160 183
</parameters>
</space>
</source>
</settings>

View file

@ -1,13 +0,0 @@
<?xml version="1.0"?>
<tallies>
<tally id="1">
<filter type="cell" bins="21" />
<scores>
flux total scatter nu-scatter scatter-2 nu-scatter-2 transport n1n
absorption nu-fission kappa-fission flux-y2 total-y2 scatter-y2
nu-scatter-y2 events delayed-nu-fission
</scores>
</tally>
</tallies>

Some files were not shown because too many files have changed in this diff Show more