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

Updating local branch with latest commits in develop.

Browse source 
Merge branch 'develop' of https://github.com/mit-crpg/openmc into fix-init-args-capi
This commit is contained in:
Jose Salcedo Perez 2018-08-03 15:20:55 +00:00
commit cd0e9667f4
84 changed files with 6329 additions and 5295 deletions
Download patch file Download diff file Expand all files Collapse all files

184
CMakeLists.txt
View file

@ -1,4 +1,4 @@
cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
cmake_minimum_required(VERSION 3.3 FATAL_ERROR)
project(openmc Fortran C CXX)
# Setup output directories
@ -9,14 +9,6 @@ set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin)
# Set module path
set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake/Modules)
#===============================================================================
# Architecture specific definitions
#===============================================================================
if (${UNIX})
add_definitions(-DUNIX)
endif()
#===============================================================================
# Command line options
#===============================================================================
@ -27,10 +19,7 @@ option(debug "Compile with debug flags" OFF)
option(optimize "Turn on all compiler optimization flags" OFF)
option(coverage "Compile with coverage analysis flags" OFF)
option(mpif08 "Use Fortran 2008 MPI interface" OFF)
# Maximum number of nested coordinates levels
set(maxcoord 10 CACHE STRING "Maximum number of nested coordinate levels")
add_definitions(-DMAX_COORD=${maxcoord})
#===============================================================================
# MPI for distributed-memory parallelism
@ -39,17 +28,12 @@ add_definitions(-DMAX_COORD=${maxcoord})
set(MPI_ENABLED FALSE)
if($ENV{FC} MATCHES "(mpi[^/]*|ftn)$")
message("-- Detected MPI wrapper: $ENV{FC}")
add_definitions(-DOPENMC_MPI)
set(MPI_ENABLED TRUE)
# Get directory containing MPI wrapper
get_filename_component(MPI_DIR $ENV{FC} DIRECTORY)
endif()
# Check for Fortran 2008 MPI interface
if(MPI_ENABLED AND mpif08)
message("-- Using Fortran 2008 MPI bindings")
add_definitions(-DOPENMC_MPIF08)
endif()
#===============================================================================
@ -77,7 +61,6 @@ if(HDF5_IS_PARALLEL)
if(NOT MPI_ENABLED)
message(FATAL_ERROR "Parallel HDF5 must be used with MPI.")
endif()
add_definitions(-DPHDF5)
message("-- Using parallel HDF5")
endif()
@ -85,19 +68,15 @@ endif()
# Set compile/link flags based on which compiler is being used
#===============================================================================
# Support for Fortran in FindOpenMP was added in CMake 3.1. To support lower
# versions, we manually add the flags. However, at some point in time, the
# manual logic can be removed in favor of the block below
#if(NOT (CMAKE_VERSION VERSION_LESS 3.1))
# if(openmp)
# find_package(OpenMP)
# if(OPENMP_FOUND)
# list(APPEND f90flags ${OpenMP_Fortran_FLAGS})
# list(APPEND ldflags ${OpenMP_Fortran_FLAGS})
# endif()
# endif()
#endif()
if(openmp)
# Requires CMake 3.1+
find_package(OpenMP)
if(OPENMP_FOUND)
list(APPEND f90flags ${OpenMP_Fortran_FLAGS})
list(APPEND cxxflags ${OpenMP_CXX_FLAGS})
list(APPEND ldflags ${OpenMP_Fortran_FLAGS})
endif()
endif()
set(CMAKE_POSITION_INDEPENDENT_CODE ON)
@ -125,10 +104,6 @@ if(CMAKE_Fortran_COMPILER_ID STREQUAL GNU)
list(REMOVE_ITEM f90flags -O2)
list(APPEND f90flags -O3)
endif()
if(openmp)
list(APPEND f90flags -fopenmp)
list(APPEND ldflags -fopenmp)
endif()
if(coverage)
list(APPEND f90flags -coverage)
list(APPEND ldflags -coverage)
@ -149,10 +124,6 @@ elseif(CMAKE_Fortran_COMPILER_ID STREQUAL Intel)
if(optimize)
list(APPEND f90flags -O3)
endif()
if(openmp)
list(APPEND f90flags -qopenmp)
list(APPEND ldflags -qopenmp)
endif()
elseif(CMAKE_Fortran_COMPILER_ID STREQUAL PGI)
# PGI Fortran compiler options
@ -187,10 +158,6 @@ elseif(CMAKE_Fortran_COMPILER_ID STREQUAL XL)
list(REMOVE_ITEM f90flags -O2)
list(APPEND f90flags -O3)
endif()
if(openmp)
list(APPEND f90flags -qsmp=omp)
list(APPEND ldflags -qsmp=omp)
endif()
elseif(CMAKE_Fortran_COMPILER_ID STREQUAL Cray)
# Cray Fortran compiler options
@ -204,7 +171,7 @@ endif()
if(CMAKE_C_COMPILER_ID STREQUAL GNU)
# GCC compiler options
list(APPEND cflags -std=c99 -O2)
list(APPEND cflags -O2)
if(debug)
list(REMOVE_ITEM cflags -O2)
list(APPEND cflags -g -Wall -pedantic -fbounds-check)
@ -222,7 +189,6 @@ if(CMAKE_C_COMPILER_ID STREQUAL GNU)
elseif(CMAKE_C_COMPILER_ID STREQUAL Intel)
# Intel compiler options
list(APPEND cflags -std=c99)
if(debug)
list(APPEND cflags -g -w3 -ftrapuv -fp-stack-check -O0)
endif()
@ -235,7 +201,6 @@ elseif(CMAKE_C_COMPILER_ID STREQUAL Intel)
elseif(CMAKE_C_COMPILER_ID MATCHES Clang)
# Clang options
list(APPEND cflags -std=c99)
if(debug)
list(APPEND cflags -g -O0 -ftrapv)
endif()
@ -257,9 +222,6 @@ if(optimize)
list(REMOVE_ITEM cxxflags -O2)
list(APPEND cxxflags -O3)
endif()
if(openmp)
list(APPEND cxxflags -fopenmp)
endif()
# Show flags being used
message(STATUS "Fortran flags: ${f90flags}")
@ -267,26 +229,12 @@ message(STATUS "C flags: ${cflags}")
message(STATUS "C++ flags: ${cxxflags}")
message(STATUS "Linker flags: ${ldflags}")
#===============================================================================
# git SHA1 hash
#===============================================================================
execute_process(COMMAND git rev-parse HEAD
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
RESULT_VARIABLE GIT_SHA1_SUCCESS
OUTPUT_VARIABLE GIT_SHA1
ERROR_QUIET OUTPUT_STRIP_TRAILING_WHITESPACE)
if(GIT_SHA1_SUCCESS EQUAL 0)
add_definitions(-DGIT_SHA1="${GIT_SHA1}")
endif()
#===============================================================================
# pugixml library
#===============================================================================
add_library(pugixml src/pugixml/pugixml_c.cpp src/pugixml/pugixml.cpp)
add_library(pugixml_fortran src/pugixml/pugixml_f.F90)
target_link_libraries(pugixml_fortran pugixml)
add_library(pugixml vendor/pugixml/pugixml.cpp)
target_include_directories(pugixml PUBLIC vendor/pugixml/)
#===============================================================================
# RPATH information
@ -295,7 +243,7 @@ target_link_libraries(pugixml_fortran pugixml)
# This block of code ensures that dynamic libraries can be found via the RPATH
# whether the executable is the original one from the build directory or the
# installed one in CMAKE_INSTALL_PREFIX. Ref:
# https://cmake.org/Wiki/CMake_RPATH_handling#Always_full_RPATH
# https://gitlab.kitware.com/cmake/community/wikis/doc/cmake/RPATH-handling
# use, i.e. don't skip the full RPATH for the build tree
set(CMAKE_SKIP_BUILD_RPATH FALSE)
@ -317,17 +265,20 @@ if("${isSystemDir}" STREQUAL "-1")
endif()
#===============================================================================
# Build faddeeva library
# faddeeva library
#===============================================================================
add_library(faddeeva STATIC src/faddeeva/Faddeeva.c)
add_library(faddeeva STATIC vendor/faddeeva/Faddeeva.c)
target_compile_options(faddeeva PRIVATE ${cflags})
set_target_properties(faddeeva PROPERTIES
C_STANDARD 99
C_STANDARD_REQUIRED ON)
#===============================================================================
# List source files. Define the libopenmc and the OpenMC executable
# libopenmc
#===============================================================================
set(program "openmc")
set(LIBOPENMC_FORTRAN_SRC
add_library(libopenmc SHARED
src/algorithm.F90
src/angle_distribution.F90
src/angleenergy_header.F90
@ -344,7 +295,6 @@ set(LIBOPENMC_FORTRAN_SRC
src/dict_header.F90
src/distribution_multivariate.F90
src/distribution_univariate.F90
src/doppler.F90
src/eigenvalue.F90
src/endf.F90
src/endf_header.F90
@ -363,7 +313,7 @@ set(LIBOPENMC_FORTRAN_SRC
src/mesh_header.F90
src/message_passing.F90
src/mgxs_data.F90
src/mgxs_header.F90
src/mgxs_interface.F90
src/multipole_header.F90
src/nuclide_header.F90
src/output.F90
@ -376,11 +326,11 @@ set(LIBOPENMC_FORTRAN_SRC
src/plot_header.F90
src/product_header.F90
src/progress_header.F90
src/pugixml/pugixml_f.F90
src/random_lcg.F90
src/reaction_header.F90
src/relaxng
src/sab_header.F90
src/scattdata_header.F90
src/secondary_correlated.F90
src/secondary_kalbach.F90
src/secondary_nbody.F90
@ -430,8 +380,6 @@ set(LIBOPENMC_FORTRAN_SRC
src/tallies/tally_header.F90
src/tallies/trigger.F90
src/tallies/trigger_header.F90
)
set(LIBOPENMC_CXX_SRC
src/cell.cpp
src/initialize.cpp
src/finalize.cpp
@ -440,53 +388,71 @@ set(LIBOPENMC_CXX_SRC
src/lattice.cpp
src/math_functions.cpp
src/message_passing.cpp
src/mgxs.cpp
src/mgxs_interface.cpp
src/plot.cpp
src/pugixml/pugixml_c.cpp
src/random_lcg.cpp
src/scattdata.cpp
src/simulation.cpp
src/state_point.cpp
src/string_functions.cpp
src/surface.cpp
src/xml_interface.cpp
src/pugixml/pugixml.cpp)
add_library(libopenmc SHARED ${LIBOPENMC_FORTRAN_SRC} ${LIBOPENMC_CXX_SRC})
src/xsdata.cpp)
set_target_properties(libopenmc PROPERTIES
OUTPUT_NAME openmc
PUBLIC_HEADER include/openmc.h)
add_executable(${program} src/main.cpp)
PUBLIC_HEADER include/openmc.h
LINKER_LANGUAGE Fortran)
#===============================================================================
# Add compiler/linker flags
#===============================================================================
set_property(TARGET ${program} libopenmc pugixml_fortran
PROPERTY LINKER_LANGUAGE Fortran)
target_include_directories(libopenmc PUBLIC include ${HDF5_INCLUDE_DIRS})
# The executable and the faddeeva package use only one language. They can be
# set via target_compile_options which accepts a list.
target_compile_options(${program} PUBLIC ${cxxflags})
target_compile_options(faddeeva PRIVATE ${cflags})
target_include_directories(libopenmc
PUBLIC include
PRIVATE ${HDF5_INCLUDE_DIRS})
# The libopenmc library has both F90 and C++ so the compile flags must be set
# file-by-file via set_source_file_properties. The compile flags must first be
# converted from lists to strings.
string(REPLACE ";" " " f90flags "${f90flags}")
string(REPLACE ";" " " cxxflags "${cxxflags}")
set_source_files_properties(${LIBOPENMC_FORTRAN_SRC} PROPERTIES COMPILE_FLAGS
${f90flags})
set_source_files_properties(${LIBOPENMC_CXX_SRC} PROPERTIES COMPILE_FLAGS
${cxxflags})
# differently depending on the language. The $<COMPILE_LANGUAGE> generator
# expression was added in CMake 3.3
target_compile_options(libopenmc PRIVATE
$<$<COMPILE_LANGUAGE:Fortran>:${f90flags}>
$<$<COMPILE_LANGUAGE:CXX>:${cxxflags}>)
# Add HDF5 library directories to link line with -L
foreach(LIBDIR ${HDF5_LIBRARY_DIRS})
list(APPEND ldflags "-L${LIBDIR}")
endforeach()
target_compile_definitions(libopenmc PRIVATE -DMAX_COORD=${maxcoord})
if (UNIX)
# Used in progress_header.F90 for calling check_isatty
target_compile_definitions(libopenmc PRIVATE -DUNIX)
endif()
if (HDF5_IS_PARALLEL)
target_compile_definitions(libopenmc PRIVATE -DPHDF5)
endif()
if (MPI_ENABLED)
target_compile_definitions(libopenmc PUBLIC -DOPENMC_MPI)
if (mpif08)
target_compile_definitions(libopenmc PRIVATE -DOPENMC_MPIF08)
endif()
endif()
# Set git SHA1 hash as a compile definition
execute_process(COMMAND git rev-parse HEAD
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
RESULT_VARIABLE GIT_SHA1_SUCCESS
OUTPUT_VARIABLE GIT_SHA1
ERROR_QUIET OUTPUT_STRIP_TRAILING_WHITESPACE)
if(GIT_SHA1_SUCCESS EQUAL 0)
target_compile_definitions(libopenmc PRIVATE -DGIT_SHA1="${GIT_SHA1}")
endif()
# target_link_libraries treats any arguments starting with - but not -l as
# linker flags. Thus, we can pass both linker flags and libraries together.
target_link_libraries(libopenmc ${ldflags} ${HDF5_LIBRARIES} pugixml_fortran
target_link_libraries(libopenmc ${ldflags} ${HDF5_LIBRARIES} pugixml
faddeeva)
target_link_libraries(${program} ${ldflags} libopenmc)
#===============================================================================
# openmc executable
#===============================================================================
add_executable(openmc src/main.cpp)
target_compile_options(openmc PRIVATE ${cxxflags})
target_link_libraries(openmc libopenmc)
#===============================================================================
# Python package
@ -502,7 +468,7 @@ add_custom_command(TARGET libopenmc POST_BUILD
# Install executable, scripts, manpage, license
#===============================================================================
install(TARGETS ${program} libopenmc
install(TARGETS openmc libopenmc
RUNTIME DESTINATION bin
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib
@ -510,4 +476,4 @@ install(TARGETS ${program} libopenmc
)
install(DIRECTORY src/relaxng DESTINATION share/openmc)
install(FILES man/man1/openmc.1 DESTINATION share/man/man1)
install(FILES LICENSE DESTINATION "share/doc/${program}" RENAME copyright)
install(FILES LICENSE DESTINATION "share/doc/openmc" RENAME copyright)

76
CODE_OF_CONDUCT.md Normal file
View file

@ -0,0 +1,76 @@
# Contributor Covenant Code of Conduct
## Our Pledge
In the interest of fostering an open and welcoming environment, we as
contributors and maintainers pledge to making participation in our project and
our community a harassment-free experience for everyone, regardless of age, body
size, disability, ethnicity, sex characteristics, gender identity and expression,
level of experience, education, socio-economic status, nationality, personal
appearance, race, religion, or sexual identity and orientation.
## Our Standards
Examples of behavior that contributes to creating a positive environment
include:
* Using welcoming and inclusive language
* Being respectful of differing viewpoints and experiences
* Gracefully accepting constructive criticism
* Focusing on what is best for the community
* Showing empathy towards other community members
Examples of unacceptable behavior by participants include:
* The use of sexualized language or imagery and unwelcome sexual attention or
advances
* Trolling, insulting/derogatory comments, and personal or political attacks
* Public or private harassment
* Publishing others' private information, such as a physical or electronic
address, without explicit permission
* Other conduct which could reasonably be considered inappropriate in a
professional setting
## Our Responsibilities
Project maintainers are responsible for clarifying the standards of acceptable
behavior and are expected to take appropriate and fair corrective action in
response to any instances of unacceptable behavior.
Project maintainers have the right and responsibility to remove, edit, or
reject comments, commits, code, wiki edits, issues, and other contributions
that are not aligned to this Code of Conduct, or to ban temporarily or
permanently any contributor for other behaviors that they deem inappropriate,
threatening, offensive, or harmful.
## Scope
This Code of Conduct applies both within project spaces and in public spaces
when an individual is representing the project or its community. Examples of
representing a project or community include using an official project e-mail
address, posting via an official social media account, or acting as an appointed
representative at an online or offline event. Representation of a project may be
further defined and clarified by project maintainers.
## Enforcement
Instances of abusive, harassing, or otherwise unacceptable behavior may be
reported by contacting the project team at openmc@anl.gov. All complaints will
be reviewed and investigated and will result in a response that is deemed
necessary and appropriate to the circumstances. The project team is obligated to
maintain confidentiality with regard to the reporter of an incident. However,
note that some project team members may have a legal obligation to report
certain forms of harassment because of their affiliation (for example, staff and
faculty at universities in the United States). Further details of specific
enforcement policies may be posted separately.
Project maintainers who do not follow or enforce the Code of Conduct in good
faith may face temporary or permanent repercussions as determined by other
members of the project's leadership.
## Attribution
This Code of Conduct is adapted from the [Contributor Covenant][homepage], version 1.4,
available at https://www.contributor-covenant.org/version/1/4/code-of-conduct.html
[homepage]: https://www.contributor-covenant.org

46
CONTRIBUTING.md Normal file
View file

@ -0,0 +1,46 @@
# Contributing to OpenMC
Welcome, and thank you for considering contributing to OpenMC! We look forward
to welcoming new members to the community and will do our best to help you get
up to speed.
## Code of Conduct
Participants in the OpenMC project are expected to follow and uphold the [Code
of Conduct](CODE_OF_CONDUCT.md). Please report any unacceptable behavior to
openmc@anl.gov.
## Resources
- [GitHub Repository](https://github.com/openmc-dev/openmc)
- [Documentation](http://openmc.readthedocs.io/en/latest)
- [User's Mailing List](openmc-users@googlegroups.com)
- [Developer's Mailing List](openmc-dev@googlegroups.com)
- [Slack Community](https://openmc.slack.com/signup) (If you don't see your
domain listed, contact openmc@anl.gov)
## How to Report Bugs
OpenMC is hosted on GitHub and all bugs are reported and tracked through the
[Issues](https://github.com/openmc-dev/openmc/issues) listed on GitHub.
## How to Suggest Enhancements
We welcome suggestions for new features or enhancements to the code and
encourage you to submit them as Issues on GitHub. However, it's important to
recognize that our development team is relatively small and does not have
unlimited time to devote to new feature suggestions. If you are interested in
working on the feature you are requesting, indicate so in the issue and the
development team will be happy to discuss it.
## How to Submit Changes
All changes to OpenMC happen through pull requests. For a full overview of the
process, see the developer's guide section on [Development
Workflow](http://openmc.readthedocs.io/en/latest/devguide/workflow.html).
## Code Style
Before you run off to make changes to the code, please have a look at our [style
guide](http://openmc.readthedocs.io/en/latest/devguide/styleguide.html), which
is used when reviewing new contributions.

2
LICENSE
View file

@ -1,4 +1,4 @@
Copyright (c) 2011-2018 Massachusetts Institute of Technology
Copyright (c) 2011-2018 Massachusetts Institute of Technology and OpenMC contributors
Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in

58
README.md Normal file
View file

@ -0,0 +1,58 @@
# OpenMC Monte Carlo Particle Transport Code
[![License](https://img.shields.io/github/license/openmc-dev/openmc.svg)](http://openmc.readthedocs.io/en/latest/license.html)
[![Travis CI build status (Linux)](https://travis-ci.org/openmc-dev/openmc.svg?branch=develop)](https://travis-ci.org/openmc-dev/openmc)
[![Code Coverage](https://coveralls.io/repos/github/openmc-dev/openmc/badge.svg?branch=develop)](https://coveralls.io/github/openmc-dev/openmc?branch=develop)
The OpenMC project aims to provide a fully-featured Monte Carlo particle
transport code based on modern methods. It is a constructive solid geometry,
continuous-energy transport code that uses HDF5 format cross sections. The
project started under the Computational Reactor Physics Group at MIT.
Complete documentation on the usage of OpenMC is hosted on Read the Docs (both
for the [latest release](http://openmc.readthedocs.io/en/stable/) and
[developmental](http://openmc.readthedocs.io/en/latest/) version). If you are
interested in the project or would like to help and contribute, please send a
message to the OpenMC User's Group [mailing
list](https://groups.google.com/forum/?fromgroups=#!forum/openmc-users).
## Installation
Detailed [installation
instructions](http://openmc.readthedocs.io/en/stable/usersguide/install.html)
can be found in the User's Guide.
## Citing
If you use OpenMC in your research, please consider giving proper attribution by
citing the following publication:
- Paul K. Romano, Nicholas E. Horelik, Bryan R. Herman, Adam G. Nelson, Benoit
Forget, and Kord Smith, "[OpenMC: A State-of-the-Art Monte Carlo Code for
Research and Development](https://doi.org/10.1016/j.anucene.2014.07.048),"
*Ann. Nucl. Energy*, **82**, 90--97 (2015).
## Troubleshooting
If you run into problems compiling, installing, or running OpenMC, first check
the [Troubleshooting
section](http://openmc.readthedocs.io/en/stable/usersguide/troubleshoot.html) in
the User's Guide. If you are not able to find a solution to your problem there,
please send a message to the User's Group [mailing
list](https://groups.google.com/forum/?fromgroups=#!forum/openmc-users).
## Reporting Bugs
OpenMC is hosted on GitHub and all bugs are reported and tracked through the
[Issues](https://github.com/openmc-dev/openmc/issues) feature on GitHub. However,
GitHub Issues should not be used for common troubleshooting purposes. If you are
having trouble installing the code or getting your model to run properly, you
should first send a message to the User's Group mailing list. If it turns out
your issue really is a bug in the code, an issue will then be created on
GitHub. If you want to request that a feature be added to the code, you may
create an Issue on github.
## License
OpenMC is distributed under the MIT/X
[license](http://openmc.readthedocs.io/en/stable/license.html).

4
docs/source/conf.py
View file

@ -71,7 +71,7 @@ master_doc = 'index'
# General information about the project.
project = u'OpenMC'
copyright = u'2011-2018, Massachusetts Institute of Technology'
copyright = u'2011-2018, Massachusetts Institute of Technology and OpenMC contributors'
# The version info for the project you're documenting, acts as replacement for
# |version| and |release|, also used in various other places throughout the
@ -208,7 +208,7 @@ htmlhelp_basename = 'openmcdoc'
# (source start file, target name, title, author, documentclass [howto/manual]).
latex_documents = [
('index', 'openmc.tex', u'OpenMC Documentation',
u'Massachusetts Institute of Technology', 'manual'),
u'OpenMC contributors', 'manual'),
]
latex_elements = {

103
docs/source/devguide/styleguide.rst
View file

@ -182,37 +182,83 @@ Always use C++-style comments (``//``) as opposed to C-style (``/**/``). (It
is more difficult to comment out a large section of code that uses C-style
comments.)
Header files should always use include guards with the following style (See
`SF.8 <http://isocpp.github.io/CppCoreGuidelines/CppCoreGuidelines#sf8-use-include-guards-for-all-h-files>`_:
.. code-block:: C++
#ifndef MODULE_NAME_H
#define MODULE_NAME_H
...
content
...
#endif // MODULE_NAME_H
Do not use C-style casting. Always use the C++-style casts ``static_cast``,
``const_cast``, or ``reinterpret_cast``. (See `ES.49 <http://isocpp.github.io/CppCoreGuidelines/CppCoreGuidelines#es49-if-you-must-use-a-cast-use-a-named-cast>`_)
Source Files
------------
Use a ``.cpp`` suffix for code files and ``.h`` for header files.
Header files should always use include guards with the following style (See
`SF.8 <http://isocpp.github.io/CppCoreGuidelines/CppCoreGuidelines#sf8-use-include-guards-for-all-h-files>`_):
.. code-block:: C++
#ifndef OPENMC_MODULE_NAME_H
#define OPENMC_MODULE_NAME_H
namespace openmc {
...
content
...
}
#endif // OPENMC_MODULE_NAME_H
Avoid hidden dependencies by always including a related header file first,
followed by C/C++ library includes, other library includes, and then local
includes. For example:
.. code-block:: C++
// foo.cpp
#include "foo.h"
#include <cstddef>
#include <iostream>
#include <vector>
#include "hdf5.h"
#include "pugixml.hpp"
#include "error.h"
#include "random_lcg.h"
Naming
------
In general, write your code in lower-case. Having code in all caps does not
enhance code readability or otherwise.
Struct and class names should be CamelCase, e.g. ``HexLattice``.
Functions (including member functions) should be lower-case with underscores,
e.g. ``get_indices``.
Local variables, global variables, and struct/class attributes should be
lower-case with underscores (e.g. ``n_cells``) except for physics symbols that
are written differently by convention (e.g. ``E`` for energy).
Local variables, global variables, and struct/class member variables should be
lower-case with underscores (e.g., ``n_cells``) except for physics symbols that
are written differently by convention (e.g., ``E`` for energy). Data members of
classes (but not structs) additionally have trailing underscores (e.g.,
``a_class_member_``).
Const variables should be in upper-case with underscores, e.g. ``SQRT_PI``.
Accessors and mutators (get and set functions) may be named like
variables. These often correspond to actual member variables, but this is not
required. For example, ``int count()`` and ``void set_count(int count)``.
Variables declared constexpr or const that have static storage duration (exist
for the duration of the program) should be upper-case with underscores,
e.g., ``SQRT_PI``.
Use C++-style declarator layout (see `NL.18
<http://isocpp.github.io/CppCoreGuidelines/CppCoreGuidelines#nl18-use-c-style-declarator-layout>`_):
pointer and reference operators in declarations should be placed adject to the
base type rather than the variable name. Avoid declaring multiple names in a
single declaration to avoid confusion:
.. code-block:: C++
T* p; // good
T& p; // good
T *p; // bad
T* p, q; // misleading
Curly braces
------------
@ -280,6 +326,13 @@ the same line or omit the braces entirely.
for (int i = 0; i < 5; i++) content();
Documentation
-------------
Classes, structs, and functions are to be annotated for the `Doxygen
<http://www.stack.nl/~dimitri/doxygen/>`_ documentation generation tool. Use the
``\`` form of Doxygen commands, e.g., ``\brief`` instead of ``@brief``.
------
Python
------
@ -288,15 +341,17 @@ Style for Python code should follow PEP8_.
Docstrings for functions and methods should follow numpydoc_ style.
Python code should work with both Python 2.7+ and Python 3.0+.
Python code should work with Python 3.4+.
Use of third-party Python packages should be limited to numpy_, scipy_, and
h5py_. Use of other third-party packages must be implemented as optional
dependencies rather than required dependencies.
Use of third-party Python packages should be limited to numpy_, scipy_,
matplotlib_, pandas_, and h5py_. Use of other third-party packages must be
implemented as optional dependencies rather than required dependencies.
.. _C++ Core Guidelines: http://isocpp.github.io/CppCoreGuidelines/CppCoreGuidelines
.. _PEP8: https://www.python.org/dev/peps/pep-0008/
.. _numpydoc: https://github.com/numpy/numpy/blob/master/doc/HOWTO_DOCUMENT.rst.txt
.. _numpy: http://www.numpy.org/
.. _scipy: http://www.scipy.org/
.. _scipy: https://www.scipy.org/
.. _matplotlib: https://matplotlib.org/
.. _pandas: https://pandas.pydata.org/
.. _h5py: http://www.h5py.org/

4
docs/source/devguide/user-input.rst
View file

@ -65,8 +65,8 @@ developer or send a message to the `developers mailing list`_.
.. _property attribute: https://docs.python.org/3.6/library/functions.html#property
.. _XML Schema Part 2: http://www.w3.org/TR/xmlschema-2/
.. _boolean: http://www.w3.org/TR/xmlschema-2/#boolean
.. _xml_interface module: https://github.com/mit-crpg/openmc/blob/develop/src/xml_interface.F90
.. _input_xml module: https://github.com/mit-crpg/openmc/blob/develop/src/input_xml.F90
.. _xml_interface module: https://github.com/openmc-dev/openmc/blob/develop/src/xml_interface.F90
.. _input_xml module: https://github.com/openmc-dev/openmc/blob/develop/src/input_xml.F90
.. _RELAX NG: http://relaxng.org/
.. _compact syntax: http://relaxng.org/compact-tutorial-20030326.html
.. _trang: http://www.thaiopensource.com/relaxng/trang.html

12
docs/source/devguide/workflow.rst
View file

@ -55,7 +55,7 @@ Now that you understand the basic development workflow, let's discuss how an
individual to contribute to development. Note that this would apply to both new
features and bug fixes. The general steps for contributing are as follows:
1. Fork the main openmc repository from `mit-crpg/openmc`_. This will create a
1. Fork the main openmc repository from `openmc-dev/openmc`_. This will create a
repository with the same name under your personal account. As such, you can
commit to it as you please without disrupting other developers.
@ -74,7 +74,7 @@ features and bug fixes. The general steps for contributing are as follows:
ensure that those changes are made on a different branch.
4. Issue a pull request from GitHub and select the *develop* branch of
mit-crpg/openmc as the target.
openmc-dev/openmc as the target.
.. image:: ../_images/pullrequest.png
@ -87,7 +87,7 @@ features and bug fixes. The general steps for contributing are as follows:
request page itself.
6. After the pull request has been thoroughly vetted, it is merged back into the
*develop* branch of mit-crpg/openmc.
*develop* branch of openmc-dev/openmc.
Private Development
-------------------
@ -99,7 +99,7 @@ create a complete copy of the OpenMC repository (not a fork from GitHub). The
private repository can then either be stored just locally or in conjunction with
a private repository on Github (this requires a `paid plan`_). Alternatively,
`Bitbucket`_ offers private repositories for free. If you want to merge some
changes you've made in your private repository back to mit-crpg/openmc
changes you've made in your private repository back to openmc-dev/openmc
repository, simply follow the steps above with an extra step of pulling a branch
from your private repository into a public fork.
@ -128,9 +128,9 @@ can interfere with virtual environments.
.. _GitHub: https://github.com/
.. _git flow: http://nvie.com/git-model
.. _valgrind: http://valgrind.org/
.. _style guide: http://mit-crpg.github.io/openmc/devguide/styleguide.html
.. _style guide: http://openmc.readthedocs.io/en/latest/devguide/styleguide.html
.. _pull request: https://help.github.com/articles/using-pull-requests
.. _mit-crpg/openmc: https://github.com/mit-crpg/openmc
.. _openmc-dev/openmc: https://github.com/openmc-dev/openmc
.. _paid plan: https://github.com/plans
.. _Bitbucket: https://bitbucket.org
.. _ctest: http://www.cmake.org/cmake/help/v2.8.12/ctest.html

1
docs/source/examples/index.rst
View file

@ -24,6 +24,7 @@ Basic Usage
triso
candu
nuclear-data
nuclear-data-resonance-covariance
------------------------------------
Multi-Group Cross Section Generation

13
docs/source/examples/nuclear-data-resonance-covariance.rst Normal file
View file

@ -0,0 +1,13 @@
.. _notebook_nuclear_data_resonance_covariance:
==================================
Nuclear Data: Resonance Covariance
==================================
.. only:: html
.. notebook:: ../../../examples/jupyter/nuclear-data-resonance-covariance.ipynb
.. only:: latex
IPython notebooks must be viewed in the online HTML documentation.

2
docs/source/license.rst
View file

@ -4,7 +4,7 @@
License Agreement
=================
Copyright © 2011-2018 Massachusetts Institute of Technology
Copyright © 2011-2018 Massachusetts Institute of Technology and OpenMC contributors
Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in

5
docs/source/pythonapi/data.rst
View file

@ -79,6 +79,11 @@ Resonance Data
openmc.data.MultiLevelBreitWigner
openmc.data.ReichMoore
openmc.data.RMatrixLimited
openmc.data.ResonanceCovariances
openmc.data.ResonanceCovarianceRange
openmc.data.SingleLevelBreitWignerCovariance
openmc.data.MultiLevelBreitWignerCovariance
openmc.data.ReichMooreCovariance
openmc.data.ParticlePair
openmc.data.SpinGroup
openmc.data.Unresolved

4
docs/source/quickinstall.rst
View file

@ -83,7 +83,7 @@ Installing from Source on Linux or Mac OS X
-------------------------------------------
All OpenMC source code is hosted on `GitHub
<https://github.com/mit-crpg/openmc>`_. If you have `git
<https://github.com/openmc-dev/openmc>`_. If you have `git
<https://git-scm.com>`_, the `gcc <https://gcc.gnu.org/>`_ compiler suite,
`CMake <http://www.cmake.org>`_, and `HDF5 <https://www.hdfgroup.org/HDF5/>`_
installed, you can download and install OpenMC be entering the following
@ -91,7 +91,7 @@ commands in a terminal:
.. code-block:: sh
git clone https://github.com/mit-crpg/openmc.git
git clone https://github.com/openmc-dev/openmc.git
cd openmc
mkdir build && cd build
cmake ..

16
docs/source/releasenotes.rst
View file

@ -71,14 +71,14 @@ Bug Fixes
- 0c6915_: Bugfix for generating thermal scattering data
- 61ecb4_: Fix bugs in Python multipole objects
.. _937469: https://github.com/mit-crpg/openmc/commit/937469
.. _a149ef: https://github.com/mit-crpg/openmc/commit/a149ef
.. _2c9b21: https://github.com/mit-crpg/openmc/commit/2c9b21
.. _8047f6: https://github.com/mit-crpg/openmc/commit/8047f6
.. _0beb4c: https://github.com/mit-crpg/openmc/commit/0beb4c
.. _f124be: https://github.com/mit-crpg/openmc/commit/f124be
.. _0c6915: https://github.com/mit-crpg/openmc/commit/0c6915
.. _61ecb4: https://github.com/mit-crpg/openmc/commit/61ecb4
.. _937469: https://github.com/openmc-dev/openmc/commit/937469
.. _a149ef: https://github.com/openmc-dev/openmc/commit/a149ef
.. _2c9b21: https://github.com/openmc-dev/openmc/commit/2c9b21
.. _8047f6: https://github.com/openmc-dev/openmc/commit/8047f6
.. _0beb4c: https://github.com/openmc-dev/openmc/commit/0beb4c
.. _f124be: https://github.com/openmc-dev/openmc/commit/f124be
.. _0c6915: https://github.com/openmc-dev/openmc/commit/0c6915
.. _61ecb4: https://github.com/openmc-dev/openmc/commit/61ecb4
------------
Contributors

4
docs/source/usersguide/beginners.rst
View file

@ -153,9 +153,9 @@ and `Volume II`_. You may also find it helpful to review the following terms:
.. _Reactor Concepts Manual: http://www.tayloredge.com/periodic/trivia/ReactorConcepts.pdf
.. _Volume I: https://www.standards.doe.gov/standards-documents/1000/1019-bhdbk-1993-v1
.. _Volume II: https://www.standards.doe.gov/standards-documents/1000/1019-bhdbk-1993-v2
.. _OpenMC source code: https://github.com/mit-crpg/openmc
.. _OpenMC source code: https://github.com/openmc-dev/openmc
.. _GitHub: https://github.com/
.. _bug reports: https://github.com/mit-crpg/openmc/issues
.. _bug reports: https://github.com/openmc-dev/openmc/issues
.. _Neutron cross section: http://en.wikipedia.org/wiki/Neutron_cross_section
.. _Effective multiplication factor: https://en.wikipedia.org/wiki/Nuclear_chain_reaction#Effective_neutron_multiplication_factor
.. _Flux: http://en.wikipedia.org/wiki/Neutron_flux

4
docs/source/usersguide/install.rst
View file

@ -181,7 +181,7 @@ with GitHub since this involves setting up ssh_ keys. With git installed and
setup, the following command will download the full source code from the GitHub
repository::
git clone https://github.com/mit-crpg/openmc.git
git clone https://github.com/openmc-dev/openmc.git
By default, the cloned repository will be set to the development branch. To
switch to the source of the latest stable release, run the following commands::
@ -189,7 +189,7 @@ switch to the source of the latest stable release, run the following commands::
cd openmc
git checkout master
.. _GitHub: https://github.com/mit-crpg/openmc
.. _GitHub: https://github.com/openmc-dev/openmc
.. _git: https://git-scm.com
.. _ssh: https://en.wikipedia.org/wiki/Secure_Shell

2
docs/source/usersguide/tallies.rst
View file

@ -213,7 +213,7 @@ The following tables show all valid scores:
+----------------------+---------------------------------------------------+
|nu-fission |Total production of neutrons due to fission. |
+----------------------+---------------------------------------------------+
|nu-scatter, |This score is similar in functionality to the |
|nu-scatter |This score is similar in functionality to the |
| |``scatter`` score except the total production of |
| |neutrons due to scattering is scored vice simply |
| |the scattering rate. This accounts for |

964
examples/jupyter/nuclear-data-resonance-covariance.ipynb Normal file
View file

File diff suppressed because one or more lines are too long

1
include/openmc.h
View file

@ -134,6 +134,7 @@ extern "C" {
extern int32_t n_lattices;
extern int32_t n_materials;
extern int32_t n_meshes;
extern int n_nuclides;
extern int64_t n_particles;
extern int32_t n_plots;
extern int32_t n_realizations;

5
man/man1/openmc.1
View file

@ -59,7 +59,8 @@ Indicates the default path to a directory containing windowed multipole data if
the user has not specified the <multipole_library> tag in
.I materials.xml\fP.
.SH LICENSE
Copyright \(co 2011-2018 Massachusetts Institute of Technology.
Copyright \(co 2011-2018 Massachusetts Institute of Technology and OpenMC
contributors.
.PP
Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in
@ -79,7 +80,7 @@ IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
.SH REPORTING BUGS
The OpenMC source code is hosted on GitHub at
https://github.com/mit-crpg/openmc. With a github account, you can submit issues
https://github.com/openmc-dev/openmc. With a github account, you can submit issues
directly on the github repository that will then be reviewed by OpenMC
developers. Alternatively, you can send a bug report to
.I openmc-users@googlegroups.com\fP.

1
openmc/data/__init__.py
View file

@ -27,5 +27,6 @@ from .urr import *
from .library import *
from .fission_energy import *
from .resonance import *
from .resonance_covariance import *
from .multipole import *
from .grid import *

64
openmc/data/endf.py
View file

@ -71,6 +71,24 @@ def float_endf(s):
return float(_ENDF_FLOAT_RE.sub(r'\1e\2', s))
def _int_endf(s):
"""Convert string to int. Used for INTG records where blank entries
indicate a 0.
Parameters
----------
s : str
Integer or spaces
Returns
-------
integer
The number or 0
"""
s = s.strip()
return int(s) if s else 0
def get_text_record(file_obj):
"""Return data from a TEXT record in an ENDF-6 file.
@ -250,6 +268,50 @@ def get_tab2_record(file_obj):
return params, Tabulated2D(breakpoints, interpolation)
def get_intg_record(file_obj):
"""
Return data from an INTG record in an ENDF-6 file. Used to store the
covariance matrix in a compact format.
Parameters
----------
file_obj : file-like object
ENDF-6 file to read from
Returns
-------
numpy.ndarray
The correlation matrix described in the INTG record
"""
# determine how many items are in list and NDIGIT
items = get_cont_record(file_obj)
ndigit = int(items[2])
npar = int(items[3]) # Number of parameters
nlines = int(items[4]) # Lines to read
NROW_RULES = {2: 18, 3: 12, 4: 11, 5: 9, 6: 8}
nrow = NROW_RULES[ndigit]
# read lines and build correlation matrix
corr = np.identity(npar)
for i in range(nlines):
line = file_obj.readline()
ii = _int_endf(line[:5]) - 1 # -1 to account for 0 indexing
jj = _int_endf(line[5:10]) - 1
factor = 10**ndigit
for j in range(nrow):
if jj+j >= ii:
break
element = _int_endf(line[11+(ndigit+1)*j:11+(ndigit+1)*(j+1)])
if element > 0:
corr[ii, jj] = (element+0.5)/factor
elif element < 0:
corr[ii, jj] = (element-0.5)/factor
# Symmetrize the correlation matrix
corr = corr + corr.T - np.diag(corr.diagonal())
return corr
def get_evaluations(filename):
"""Return a list of all evaluations within an ENDF file.
@ -288,7 +350,7 @@ class Evaluation(object):
Attributes
----------
info : dict
Miscallaneous information about the evaluation.
Miscellaneous information about the evaluation.
target : dict
Information about the target material, such as its mass, isomeric state,
whether it's stable, and whether it's fissionable.

24
openmc/data/neutron.py
View file

@ -24,6 +24,7 @@ from .njoy import make_ace
from .product import Product
from .reaction import Reaction, _get_photon_products_ace
from . import resonance as res
from . import resonance_covariance as res_cov
from .urr import ProbabilityTables
import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
@ -148,6 +149,8 @@ class IncidentNeutron(EqualityMixin):
and the values are Reaction objects.
resonances : openmc.data.Resonances or None
Resonance parameters
resonance_covariance : openmc.data.ResonanceCovariance or None
Covariance for resonance parameters
summed_reactions : collections.OrderedDict
Contains summed cross sections, e.g., the total cross section. The keys
are the MT values and the values are Reaction objects.
@ -228,6 +231,10 @@ class IncidentNeutron(EqualityMixin):
def resonances(self):
return self._resonances
@property
def resonance_covariance(self):
return self._resonance_covariance
@property
def summed_reactions(self):
return self._summed_reactions
@ -289,6 +296,12 @@ class IncidentNeutron(EqualityMixin):
cv.check_type('resonances', resonances, res.Resonances)
self._resonances = resonances
@resonance_covariance.setter
def resonance_covariance(self, resonance_covariance):
cv.check_type('resonance covariance', resonance_covariance,
res_cov.ResonanceCovariances)
self._resonance_covariance = resonance_covariance
@summed_reactions.setter
def summed_reactions(self, summed_reactions):
cv.check_type('summed reactions', summed_reactions, Mapping)
@ -744,7 +757,7 @@ class IncidentNeutron(EqualityMixin):
return data
@classmethod
def from_endf(cls, ev_or_filename):
def from_endf(cls, ev_or_filename, covariance=False):
"""Generate incident neutron continuous-energy data from an ENDF evaluation
Parameters
@ -753,6 +766,10 @@ class IncidentNeutron(EqualityMixin):
ENDF evaluation to read from. If given as a string, it is assumed to
be the filename for the ENDF file.
covariance : bool
Flag to indicate whether or not covariance data from File 32 should be
retrieved
Returns
-------
openmc.data.IncidentNeutron
@ -784,6 +801,11 @@ class IncidentNeutron(EqualityMixin):
if (2, 151) in ev.section:
data.resonances = res.Resonances.from_endf(ev)
if (32, 151) in ev.section and covariance:
data.resonance_covariance = (
res_cov.ResonanceCovariances.from_endf(ev, data.resonances)
)
# Read each reaction
for mf, mt, nc, mod in ev.reaction_list:
if mf == 3:

20
openmc/data/resonance.py
View file

@ -16,6 +16,7 @@ except ImportError:
_reconstruct = False
import openmc.checkvalue as cv
class Resonances(object):
"""Resolved and unresolved resonance data
@ -90,14 +91,14 @@ class Resonances(object):
# Determine whether discrete or continuous representation
items = get_head_record(file_obj)
n_isotope = items[4] # Number of isotopes
n_isotope = items[4] # Number of isotopes
ranges = []
for iso in range(n_isotope):
items = get_cont_record(file_obj)
abundance = items[1]
fission_widths = (items[3] == 1) # fission widths are given?
n_ranges = items[4] # number of resonance energy ranges
fission_widths = (items[3] == 1) # fission widths are given?
n_ranges = items[4] # number of resonance energy ranges
for j in range(n_ranges):
items = get_cont_record(file_obj)
@ -112,7 +113,7 @@ class Resonances(object):
# unresolved resonance region
erange = Unresolved.from_endf(file_obj, items, fission_widths)
#erange.material = self
# erange.material = self
ranges.append(erange)
return cls(ranges)
@ -162,6 +163,13 @@ class ResonanceRange(object):
self._prepared = False
self._parameter_matrix = {}
def __copy__(self):
cls = type(self)
new_copy = cls.__new__(cls)
new_copy.__dict__.update(self.__dict__)
new_copy._prepared = False
return new_copy
@classmethod
def from_endf(cls, ev, file_obj, items):
"""Create resonance range from an ENDF evaluation.
@ -437,7 +445,7 @@ class MultiLevelBreitWigner(ResonanceRange):
self._l_values = np.array(l_values)
self._competitive = np.array(competitive)
for l in l_values:
self._parameter_matrix[l] = df[df.L == l].as_matrix()
self._parameter_matrix[l] = df[df.L == l].values
self._prepared = True
@ -682,7 +690,7 @@ class ReichMoore(ResonanceRange):
self._l_values = np.array(l_values)
for (l, J) in lj_values:
self._parameter_matrix[l, J] = df[(df.L == l) &
(abs(df.J) == J)].as_matrix()
(abs(df.J) == J)].values
self._prepared = True

708
openmc/data/resonance_covariance.py Normal file
View file

@ -0,0 +1,708 @@
from collections import MutableSequence
import warnings
import io
import copy
import numpy as np
import pandas as pd
from . import endf
import openmc.checkvalue as cv
from .resonance import Resonances
def _add_file2_contributions(file32params, file2params):
"""Function for aiding in adding resonance parameters from File 2 that are
not always present in File 32. Uses already imported resonance data.
Paramaters
----------
file32params : pandas.Dataframe
Incomplete set of resonance parameters contained in File 32.
file2params : pandas.Dataframe
Resonance parameters from File 2. Ordered by energy.
Returns
-------
parameters : pandas.Dataframe
Complete set of parameters ordered by L-values and then energy
"""
# Use l-values and competitiveWidth from File 2 data
# Re-sort File 2 by energy to match File 32
file2params = file2params.sort_values(by=['energy'])
file2params.reset_index(drop=True, inplace=True)
# Sort File 32 parameters by energy as well (maintaining index)
file32params.sort_values(by=['energy'], inplace=True)
# Add in values (.values converts to array first to ignore index)
file32params['L'] = file2params['L'].values
if 'competitiveWidth' in file2params.columns:
file32params['competitiveWidth'] = file2params['competitiveWidth'].values
# Resort to File 32 order (by L then by E) for use with covariance
file32params.sort_index(inplace=True)
return file32params
class ResonanceCovariances(Resonances):
"""Resolved resonance covariance data
Parameters
----------
ranges : list of openmc.data.ResonanceCovarianceRange
Distinct energy ranges for resonance data
Attributes
----------
ranges : list of openmc.data.ResonanceCovarianceRange
Distinct energy ranges for resonance data
"""
@property
def ranges(self):
return self._ranges
@ranges.setter
def ranges(self, ranges):
cv.check_type('resonance ranges', ranges, MutableSequence)
self._ranges = cv.CheckedList(ResonanceCovarianceRange,
'resonance range', ranges)
@classmethod
def from_endf(cls, ev, resonances):
"""Generate resonance covariance data from an ENDF evaluation.
Parameters
----------
ev : openmc.data.endf.Evaluation
ENDF evaluation
resonances : openmc.data.Resonance object
openmc.data.Resonanance object generated from the same evaluation
used to import values not contained in File 32
Returns
-------
openmc.data.ResonanceCovariances
Resonance covariance data
"""
file_obj = io.StringIO(ev.section[32, 151])
# Determine whether discrete or continuous representation
items = endf.get_head_record(file_obj)
n_isotope = items[4] # Number of isotopes
ranges = []
for iso in range(n_isotope):
items = endf.get_cont_record(file_obj)
abundance = items[1]
fission_widths = (items[3] == 1) # Flag for fission widths
n_ranges = items[4] # Number of resonance energy ranges
for j in range(n_ranges):
items = endf.get_cont_record(file_obj)
# Unresolved flags - 0: only scattering radius given
# 1: resolved parameters given
# 2: unresolved parameters given
unresolved_flag = items[2]
formalism = items[3] # resonance formalism
# Throw error for unsupported formalisms
if formalism in [0, 7]:
error = 'LRF='+str(formalism)+' covariance not supported '\
'for this formalism'
raise NotImplementedError(error)
if unresolved_flag in (0, 1):
# Resolved resonance region
resonance = resonances.ranges[j]
erange = _FORMALISMS[formalism].from_endf(ev, file_obj,
items, resonance)
ranges.append(erange)
elif unresolved_flag == 2:
warn = 'Unresolved resonance not supported. Covariance '\
'values for the unresolved region not imported.'
warnings.warn(warn)
return cls(ranges)
class ResonanceCovarianceRange:
"""Resonace covariance range. Base class for different formalisms.
Parameters
----------
energy_min : float
Minimum energy of the resolved resonance range in eV
energy_max : float
Maximum energy of the resolved resonance range in eV
Attributes
----------
energy_min : float
Minimum energy of the resolved resonance range in eV
energy_max : float
Maximum energy of the resolved resonance range in eV
parameters : pandas.DataFrame
Resonance parameters
covariance : numpy.array
The covariance matrix contained within the ENDF evaluation
lcomp : int
Flag indicating format of the covariance matrix within the ENDF file
file2res : openmc.data.ResonanceRange object
Corresponding resonance range with File 2 data.
mpar : int
Number of parameters in covariance matrix for each individual resonance
formalism : str
String descriptor of formalism
"""
def __init__(self, energy_min, energy_max):
self.energy_min = energy_min
self.energy_max = energy_max
def subset(self, parameter_str, bounds):
"""Produce a subset of resonance parameters and the corresponding
covariance matrix to an IncidentNeutron object.
Parameters
----------
parameter_str : str
parameter to be discriminated
(i.e. 'energy', 'captureWidth', 'fissionWidthA'...)
bounds : np.array
[low numerical bound, high numerical bound]
Returns
-------
res_cov_range : openmc.data.ResonanceCovarianceRange
ResonanceCovarianceRange object that contains a subset of the
covariance matrix (upper triangular) as well as a subset parameters
within self.file2params
"""
# Copy range and prevent change of original
res_cov_range = copy.deepcopy(self)
parameters = self.file2res.parameters
cov = res_cov_range.covariance
mpar = res_cov_range.mpar
# Create mask
mask1 = parameters[parameter_str] >= bounds[0]
mask2 = parameters[parameter_str] <= bounds[1]
mask = mask1 & mask2
res_cov_range.parameters = parameters[mask]
indices = res_cov_range.parameters.index.values
# Build subset of covariance
sub_cov_dim = len(indices)*mpar
cov_subset_vals = []
for index1 in indices:
for i in range(mpar):
for index2 in indices:
for j in range(mpar):
if index2*mpar+j >= index1*mpar+i:
cov_subset_vals.append(cov[index1*mpar+i,
index2*mpar+j])
cov_subset = np.zeros([sub_cov_dim, sub_cov_dim])
tri_indices = np.triu_indices(sub_cov_dim)
cov_subset[tri_indices] = cov_subset_vals
res_cov_range.file2res.parameters = parameters[mask]
res_cov_range.covariance = cov_subset
return res_cov_range
def sample(self, n_samples):
"""Sample resonance parameters based on the covariances provided
within an ENDF evaluation.
Parameters
----------
n_samples : int
The number of samples to produce
Returns
-------
samples : list of openmc.data.ResonanceCovarianceRange objects
List of samples size `n_samples`
"""
warn_str = 'Sampling routine does not guarantee positive values for '\
'parameters. This can lead to undefined behavior in the '\
'reconstruction routine.'
warnings.warn(warn_str)
parameters = self.parameters
cov = self.covariance
# Symmetrizing covariance matrix
cov = cov + cov.T - np.diag(cov.diagonal())
formalism = self.formalism
mpar = self.mpar
samples = []
# Handling MLBW/SLBW sampling
if formalism == 'mlbw' or formalism == 'slbw':
params = ['energy', 'neutronWidth', 'captureWidth', 'fissionWidth',
'competitiveWidth']
param_list = params[:mpar]
mean_array = parameters[param_list].values
mean = mean_array.flatten()
par_samples = np.random.multivariate_normal(mean, cov,
size=n_samples)
spin = parameters['J'].values
l_value = parameters['L'].values
for sample in par_samples:
energy = sample[0::mpar]
gn = sample[1::mpar]
gg = sample[2::mpar]
gf = sample[3::mpar] if mpar > 3 else parameters['fissionWidth'].values
gx = sample[4::mpar] if mpar > 4 else parameters['competitiveWidth'].values
gt = gn + gg + gf + gx
records = []
for j, E in enumerate(energy):
records.append([energy[j], l_value[j], spin[j], gt[j],
gn[j], gg[j], gf[j], gx[j]])
columns = ['energy', 'L', 'J', 'totalWidth', 'neutronWidth',
'captureWidth', 'fissionWidth', 'competitiveWidth']
sample_params = pd.DataFrame.from_records(records,
columns=columns)
# Copy ResonanceRange object
res_range = copy.copy(self.file2res)
res_range.parameters = sample_params
samples.append(res_range)
# Handling RM sampling
elif formalism == 'rm':
params = ['energy', 'neutronWidth', 'captureWidth',
'fissionWidthA', 'fissionWidthB']
param_list = params[:mpar]
mean_array = parameters[param_list].values
mean = mean_array.flatten()
par_samples = np.random.multivariate_normal(mean, cov,
size=n_samples)
spin = parameters['J'].values
l_value = parameters['L'].values
for sample in par_samples:
energy = sample[0::mpar]
gn = sample[1::mpar]
gg = sample[2::mpar]
gfa = sample[3::mpar] if mpar > 3 else parameters['fissionWidthA'].values
gfb = sample[4::mpar] if mpar > 3 else parameters['fissionWidthB'].values
records = []
for j, E in enumerate(energy):
records.append([energy[j], l_value[j], spin[j], gn[j],
gg[j], gfa[j], gfb[j]])
columns = ['energy', 'L', 'J', 'neutronWidth',
'captureWidth', 'fissionWidthA', 'fissionWidthB']
sample_params = pd.DataFrame.from_records(records,
columns=columns)
# Copy ResonanceRange object
res_range = copy.copy(self.file2res)
res_range.parameters = sample_params
samples.append(res_range)
return samples
class MultiLevelBreitWignerCovariance(ResonanceCovarianceRange):
"""Multi-level Breit-Wigner resolved resonance formalism covariance data.
Parameters
----------
energy_min : float
Minimum energy of the resolved resonance range in eV
energy_max : float
Maximum energy of the resolved resonance range in eV
Attributes
----------
energy_min : float
Minimum energy of the resolved resonance range in eV
energy_max : float
Maximum energy of the resolved resonance range in eV
parameters : pandas.DataFrame
Resonance parameters
covariance : numpy.array
The covariance matrix contained within the ENDF evaluation
mpar : int
Number of parameters in covariance matrix for each individual resonance
lcomp : int
Flag indicating format of the covariance matrix within the ENDF file
file2res : openmc.data.ResonanceRange object
Corresponding resonance range with File 2 data.
formalism : str
String descriptor of formalism
"""
def __init__(self, energy_min, energy_max, parameters, covariance, mpar,
lcomp, file2res):
super().__init__(energy_min, energy_max)
self.parameters = parameters
self.covariance = covariance
self.mpar = mpar
self.lcomp = lcomp
self.file2res = copy.copy(file2res)
self.formalism = 'mlbw'
@classmethod
def from_endf(cls, ev, file_obj, items, resonance):
"""Create MLBW covariance data from an ENDF evaluation.
Parameters
----------
ev : openmc.data.endf.Evaluation
ENDF evaluation
file_obj : file-like object
ENDF file positioned at the second record of a resonance range
subsection in MF=32, MT=151
items : list
Items from the CONT record at the start of the resonance range
subsection
resonance : openmc.data.ResonanceRange object
Corresponding resonance range with File 2 data.
Returns
-------
openmc.data.MultiLevelBreitWignerCovariance
Multi-level Breit-Wigner resonance covariance parameters
"""
# Read energy-dependent scattering radius if present
energy_min, energy_max = items[0:2]
nro, naps = items[4:6]
if nro != 0:
params, ape = endf.get_tab1_record(file_obj)
# Other scatter radius parameters
items = endf.get_cont_record(file_obj)
target_spin = items[0]
lcomp = items[3] # Flag for compatibility 0, 1, 2 - 2 is compact form
nls = items[4] # number of l-values
# Build covariance matrix for General Resolved Resonance Formats
if lcomp == 1:
items = endf.get_cont_record(file_obj)
# Number of short range type resonance covariances
num_short_range = items[4]
# Number of long range type resonance covariances
num_long_range = items[5]
# Read resonance widths, J values, etc
records = []
for i in range(num_short_range):
items, values = endf.get_list_record(file_obj)
mpar = items[2]
num_res = items[5]
num_par_vals = num_res*6
res_values = values[:num_par_vals]
cov_values = values[num_par_vals:]
energy = res_values[0::6]
spin = res_values[1::6]
gt = res_values[2::6]
gn = res_values[3::6]
gg = res_values[4::6]
gf = res_values[5::6]
for i, E in enumerate(energy):
records.append([energy[i], spin[i], gt[i], gn[i],
gg[i], gf[i]])
# Build the upper-triangular covariance matrix
cov_dim = mpar*num_res
cov = np.zeros([cov_dim, cov_dim])
indices = np.triu_indices(cov_dim)
cov[indices] = cov_values
# Compact format - Resonances and individual uncertainties followed by
# compact correlations
elif lcomp == 2:
items, values = endf.get_list_record(file_obj)
mean = items
num_res = items[5]
energy = values[0::12]
spin = values[1::12]
gt = values[2::12]
gn = values[3::12]
gg = values[4::12]
gf = values[5::12]
par_unc = []
for i in range(num_res):
res_unc = values[i*12+6 : i*12+12]
# Delete 0 values (not provided, no fission width)
# DAJ/DGT always zero, DGF sometimes nonzero [1, 2, 5]
res_unc_nonzero = []
for j in range(6):
if j in [1, 2, 5] and res_unc[j] != 0.0:
res_unc_nonzero.append(res_unc[j])
elif j in [0, 3, 4]:
res_unc_nonzero.append(res_unc[j])
par_unc.extend(res_unc_nonzero)
records = []
for i, E in enumerate(energy):
records.append([energy[i], spin[i], gt[i], gn[i],
gg[i], gf[i]])
corr = endf.get_intg_record(file_obj)
cov = np.diag(par_unc).dot(corr).dot(np.diag(par_unc))
# Compatible resolved resonance format
elif lcomp == 0:
cov = np.zeros([4, 4])
records = []
cov_index = 0
for i in range(nls):
items, values = endf.get_list_record(file_obj)
num_res = items[5]
for j in range(num_res):
one_res = values[18*j:18*(j+1)]
res_values = one_res[:6]
cov_values = one_res[6:]
records.append(list(res_values))
# Populate the coviariance matrix for this resonance
# There are no covariances between resonances in lcomp=0
cov[cov_index, cov_index] = cov_values[0]
cov[cov_index+1, cov_index+1 : cov_index+2] = cov_values[1:2]
cov[cov_index+1, cov_index+3] = cov_values[4]
cov[cov_index+2, cov_index+2] = cov_values[3]
cov[cov_index+2, cov_index+3] = cov_values[5]
cov[cov_index+3, cov_index+3] = cov_values[6]
cov_index += 4
if j < num_res-1: # Pad matrix for additional values
cov = np.pad(cov, ((0, 4), (0, 4)), 'constant',
constant_values=0)
# Create pandas DataFrame with resonance data, currently
# redundant with data.IncidentNeutron.resonance
columns = ['energy', 'J', 'totalWidth', 'neutronWidth',
'captureWidth', 'fissionWidth']
parameters = pd.DataFrame.from_records(records, columns=columns)
# Determine mpar (number of parameters for each resonance in
# covariance matrix)
nparams, params = parameters.shape
covsize = cov.shape[0]
mpar = int(covsize/nparams)
# Add parameters from File 2
parameters = _add_file2_contributions(parameters,
resonance.parameters)
# Create instance of class
mlbw = cls(energy_min, energy_max, parameters, cov, mpar, lcomp,
resonance)
return mlbw
class SingleLevelBreitWignerCovariance(MultiLevelBreitWignerCovariance):
"""Single-level Breit-Wigner resolved resonance formalism covariance data.
Single-level Breit-Wigner resolved resonance data is is identified by LRF=1
in the ENDF-6 format.
Parameters
----------
energy_min : float
Minimum energy of the resolved resonance range in eV
energy_max : float
Maximum energy of the resolved resonance range in eV
Attributes
----------
energy_min : float
Minimum energy of the resolved resonance range in eV
energy_max : float
Maximum energy of the resolved resonance range in eV
parameters : pandas.DataFrame
Resonance parameters
covariance : numpy.array
The covariance matrix contained within the ENDF evaluation
mpar : int
Number of parameters in covariance matrix for each individual resonance
formalism : str
String descriptor of formalism
lcomp : int
Flag indicating format of the covariance matrix within the ENDF file
file2res : openmc.data.ResonanceRange object
Corresponding resonance range with File 2 data.
"""
def __init__(self, energy_min, energy_max, parameters, covariance, mpar,
lcomp, file2res):
super().__init__(energy_min, energy_max, parameters, covariance, mpar,
lcomp, file2res)
self.formalism = 'slbw'
class ReichMooreCovariance(ResonanceCovarianceRange):
"""Reich-Moore resolved resonance formalism covariance data.
Reich-Moore resolved resonance data is identified by LRF=3 in the ENDF-6
format.
Parameters
----------
energy_min : float
Minimum energy of the resolved resonance range in eV
energy_max : float
Maximum energy of the resolved resonance range in eV
Attributes
----------
energy_min : float
Minimum energy of the resolved resonance range in eV
energy_max : float
Maximum energy of the resolved resonance range in eV
parameters : pandas.DataFrame
Resonance parameters
covariance : numpy.array
The covariance matrix contained within the ENDF evaluation
lcomp : int
Flag indicating format of the covariance matrix within the ENDF file
mpar : int
Number of parameters in covariance matrix for each individual resonance
file2res : openmc.data.ResonanceRange object
Corresponding resonance range with File 2 data.
formalism : str
String descriptor of formalism
"""
def __init__(self, energy_min, energy_max, parameters, covariance, mpar,
lcomp, file2res):
super().__init__(energy_min, energy_max)
self.parameters = parameters
self.covariance = covariance
self.mpar = mpar
self.lcomp = lcomp
self.file2res = copy.copy(file2res)
self.formalism = 'rm'
@classmethod
def from_endf(cls, ev, file_obj, items, resonance):
"""Create Reich-Moore resonance covariance data from an ENDF
evaluation. Includes the resonance parameters contained separately in
File 32.
Parameters
----------
ev : openmc.data.endf.Evaluation
ENDF evaluation
file_obj : file-like object
ENDF file positioned at the second record of a resonance range
subsection in MF=2, MT=151
items : list
Items from the CONT record at the start of the resonance range
subsection
resonance : openmc.data.Resonance object
openmc.data.Resonanance object generated from the same evaluation
used to import values not contained in File 32
Returns
-------
openmc.data.ReichMooreCovariance
Reich-Moore resonance covariance parameters
"""
# Read energy-dependent scattering radius if present
energy_min, energy_max = items[0:2]
nro, naps = items[4:6]
if nro != 0:
params, ape = endf.get_tab1_record(file_obj)
# Other scatter radius parameters
items = endf.get_cont_record(file_obj)
target_spin = items[0]
lcomp = items[3] # Flag for compatibility 0, 1, 2 - 2 is compact form
nls = items[4] # Number of l-values
# Build covariance matrix for General Resolved Resonance Formats
if lcomp == 1:
items = endf.get_cont_record(file_obj)
# Number of short range type resonance covariances
num_short_range = items[4]
# Number of long range type resonance covariances
num_long_range = items[5]
# Read resonance widths, J values, etc
channel_radius = {}
scattering_radius = {}
records = []
for i in range(num_short_range):
items, values = endf.get_list_record(file_obj)
mpar = items[2]
num_res = items[5]
num_par_vals = num_res*6
res_values = values[:num_par_vals]
cov_values = values[num_par_vals:]
energy = res_values[0::6]
spin = res_values[1::6]
gn = res_values[2::6]
gg = res_values[3::6]
gfa = res_values[4::6]
gfb = res_values[5::6]
for i, E in enumerate(energy):
records.append([energy[i], spin[i], gn[i], gg[i],
gfa[i], gfb[i]])
# Build the upper-triangular covariance matrix
cov_dim = mpar*num_res
cov = np.zeros([cov_dim, cov_dim])
indices = np.triu_indices(cov_dim)
cov[indices] = cov_values
# Compact format - Resonances and individual uncertainties followed by
# compact correlations
elif lcomp == 2:
items, values = endf.get_list_record(file_obj)
num_res = items[5]
energy = values[0::12]
spin = values[1::12]
gn = values[2::12]
gg = values[3::12]
gfa = values[4::12]
gfb = values[5::12]
par_unc = []
for i in range(num_res):
res_unc = values[i*12+6 : i*12+12]
# Delete 0 values (not provided in evaluation)
res_unc = [x for x in res_unc if x != 0.0]
par_unc.extend(res_unc)
records = []
for i, E in enumerate(energy):
records.append([energy[i], spin[i], gn[i], gg[i],
gfa[i], gfb[i]])
corr = endf.get_intg_record(file_obj)
cov = np.diag(par_unc).dot(corr).dot(np.diag(par_unc))
# Create pandas DataFrame with resonacne data
columns = ['energy', 'J', 'neutronWidth', 'captureWidth',
'fissionWidthA', 'fissionWidthB']
parameters = pd.DataFrame.from_records(records, columns=columns)
# Determine mpar (number of parameters for each resonance in
# covariance matrix)
nparams, params = parameters.shape
covsize = cov.shape[0]
mpar = int(covsize/nparams)
# Add parameters from File 2
parameters = _add_file2_contributions(parameters,
resonance.parameters)
# Create instance of ReichMooreCovariance
rmc = cls(energy_min, energy_max, parameters, cov, mpar, lcomp,
resonance)
return rmc
_FORMALISMS = {
0: ResonanceCovarianceRange,
1: SingleLevelBreitWignerCovariance,
2: MultiLevelBreitWignerCovariance,
3: ReichMooreCovariance
# 7: RMatrixLimitedCovariance
}

60
openmc/data/thermal.py
View file

@ -24,40 +24,44 @@ from openmc.stats import Discrete, Tabular
_THERMAL_NAMES = {
'c_Al27': ('al', 'al27'),
'c_Be': ('be', 'be-metal'),
'c_BeO': ('beo'),
'c_Be_in_BeO': ('bebeo', 'be-o', 'be/o'),
'c_Al27': ('al', 'al27', 'al-27'),
'c_Be': ('be', 'be-metal', 'be-met'),
'c_BeO': ('beo',),
'c_Be_in_BeO': ('bebeo', 'be-beo', 'be-o', 'be/o'),
'c_C6H6': ('benz', 'c6h6'),
'c_C_in_SiC': ('csic',),
'c_Ca_in_CaH2': ('cah'),
'c_D_in_D2O': ('dd2o', 'hwtr', 'hw'),
'c_Fe56': ('fe', 'fe56'),
'c_C_in_SiC': ('csic', 'c-sic'),
'c_Ca_in_CaH2': ('cah',),
'c_D_in_D2O': ('dd2o', 'd-d2o', 'hwtr', 'hw'),
'c_Fe56': ('fe', 'fe56', 'fe-56'),
'c_Graphite': ('graph', 'grph', 'gr'),
'c_H_in_CaH2': ('hcah2'),
'c_H_in_CH2': ('hch2', 'poly', 'pol'),
'c_Graphite_10p': ('grph10',),
'c_Graphite_30p': ('grph30',),
'c_H_in_CaH2': ('hcah2',),
'c_H_in_CH2': ('hch2', 'poly', 'pol', 'h-poly'),
'c_H_in_CH4_liquid': ('lch4', 'lmeth'),
'c_H_in_CH4_solid': ('sch4', 'smeth'),
'c_H_in_H2O': ('hh2o', 'lwtr', 'lw'),
'c_H_in_H2O_solid': ('hice',),
'c_H_in_C5O2H8': ('lucite', 'c5o2h8'),
'c_H_in_YH2': ('hyh2'),
'c_H_in_ZrH': ('hzrh', 'h-zr', 'h/zr', 'hzr'),
'c_H_in_H2O': ('hh2o', 'h-h2o', 'lwtr', 'lw'),
'c_H_in_H2O_solid': ('hice', 'h-ice'),
'c_H_in_C5O2H8': ('lucite', 'c5o2h8', 'h-luci'),
'c_H_in_YH2': ('hyh2', 'h-yh2'),
'c_H_in_ZrH': ('hzrh', 'h-zrh', 'h-zr', 'h/zr', 'hzr'),
'c_Mg24': ('mg', 'mg24'),
'c_O_in_BeO': ('obeo', 'o-be', 'o/be'),
'c_O_in_D2O': ('od2o'),
'c_O_in_H2O_ice': ('oice'),
'c_O_in_UO2': ('ouo2', 'o2-u', 'o2/u'),
'c_ortho_D': ('orthod', 'dortho'),
'c_ortho_H': ('orthoh', 'hortho'),
'c_Si_in_SiC': ('sisic'),
'c_O_in_BeO': ('obeo', 'o-beo', 'o-be', 'o/be'),
'c_O_in_D2O': ('od2o', 'o-d2o'),
'c_O_in_H2O_ice': ('oice', 'o-ice'),
'c_O_in_UO2': ('ouo2', 'o-uo2', 'o2-u', 'o2/u'),
'c_N_in_UN': ('n-un',),
'c_ortho_D': ('orthod', 'orthoD', 'dortho'),
'c_ortho_H': ('orthoh', 'orthoH', 'hortho'),
'c_Si_in_SiC': ('sisic', 'si-sic'),
'c_SiO2_alpha': ('sio2', 'sio2a'),
'c_SiO2_beta': ('sio2b'),
'c_para_D': ('parad', 'dpara'),
'c_para_H': ('parah', 'hpara'),
'c_U_in_UO2': ('uuo2', 'u-o2', 'u/o2'),
'c_Y_in_YH2': ('yyh2'),
'c_Zr_in_ZrH': ('zrzrh', 'zr-h', 'zr/h')
'c_SiO2_beta': ('sio2b',),
'c_para_D': ('parad', 'paraD', 'dpara'),
'c_para_H': ('parah', 'paraH', 'hpara'),
'c_U_in_UN': ('u-un',),
'c_U_in_UO2': ('uuo2', 'u-uo2', 'u-o2', 'u/o2'),
'c_Y_in_YH2': ('yyh2', 'y-yh2'),
'c_Zr_in_ZrH': ('zrzrh', 'zr-zrh', 'zr-h', 'zr/h')
}

2
openmc/filter_expansion.py
View file

@ -54,7 +54,7 @@ class LegendreFilter(ExpansionFilter):
r"""Score Legendre expansion moments up to specified order.
This filter allows scores to be multiplied by Legendre polynomials of the
change in particle angle ($\mu$) up to a user-specified order.
change in particle angle (:math:`\mu`) up to a user-specified order.
Parameters
----------

10
openmc/mesh.py
View file

@ -208,12 +208,12 @@ class Mesh(IDManagerMixin):
shape = np.array(lattice.shape)
width = lattice.pitch*shape
mesh = cls(mesh_id, name)
mesh.lower_left = lattice.lower_left
mesh.upper_right = lattice.lower_left + width
mesh.dimension = shape*division
return mesh
def to_xml_element(self):
@ -333,14 +333,16 @@ class Mesh(IDManagerMixin):
if n_dim == 1:
universe_array = np.array([universes])
elif n_dim == 2:
universe_array = np.empty(self.dimension, dtype=openmc.Universe)
universe_array = np.empty(self.dimension[::-1],
dtype=openmc.Universe)
i = 0
for y in range(self.dimension[1] - 1, -1, -1):
for x in range(self.dimension[0]):
universe_array[y][x] = universes[i]
i += 1
else:
universe_array = np.empty(self.dimension, dtype=openmc.Universe)
universe_array = np.empty(self.dimension[::-1],
dtype=openmc.Universe)
i = 0
for z in range(self.dimension[2]):
for y in range(self.dimension[1] - 1, -1, -1):

80
readme.rst
View file

@ -1,80 +0,0 @@
==========================================
OpenMC Monte Carlo Particle Transport Code
==========================================
|licensebadge| |travisbadge| |coverallsbadge|
The OpenMC project aims to provide a fully-featured Monte Carlo particle
transport code based on modern methods. It is a constructive solid geometry,
continuous-energy transport code that uses HDF5 format cross sections. The
project started under the Computational Reactor Physics Group at MIT.
Complete documentation on the usage of OpenMC is hosted on Read the Docs (both
for the `latest release`_ and developmental_ version). If you are interested in
the project or would like to help and contribute, please send a message to the
OpenMC User's Group `mailing list`_.
------------
Installation
------------
Detailed `installation instructions`_ can be found in the User's Guide.
------
Citing
------
If you use OpenMC in your research, please consider giving proper attribution by
citing the following publication:
- Paul K. Romano, Nicholas E. Horelik, Bryan R. Herman, Adam G. Nelson, Benoit
Forget, and Kord Smith, "`OpenMC: A State-of-the-Art Monte Carlo Code for
Research and Development <https://doi.org/10.1016/j.anucene.2014.07.048>`_,"
*Ann. Nucl. Energy*, **82**, 90--97 (2015).
---------------
Troubleshooting
---------------
If you run into problems compiling, installing, or running OpenMC, first check
the `Troubleshooting section`_ in the User's Guide. If you are not able to find
a solution to your problem there, please send a message to the User's Group
`mailing list`_.
--------------
Reporting Bugs
--------------
OpenMC is hosted on GitHub and all bugs are reported and tracked through the
Issues_ feature on GitHub. However, GitHub Issues should not be used for common
troubleshooting purposes. If you are having trouble installing the code or
getting your model to run properly, you should first send a message to the
User's Group `mailing list`_. If it turns out your issue really is a bug in the
code, an issue will then be created on GitHub. If you want to request that a
feature be added to the code, you may create an Issue on github.
-------
License
-------
OpenMC is distributed under the MIT/X license_.
.. _latest release: http://openmc.readthedocs.io/en/stable/
.. _developmental: http://openmc.readthedocs.io/en/latest/
.. _mailing list: https://groups.google.com/forum/?fromgroups=#!forum/openmc-users
.. _installation instructions: http://openmc.readthedocs.io/en/stable/usersguide/install.html
.. _Troubleshooting section: http://openmc.readthedocs.io/en/stable/usersguide/troubleshoot.html
.. _Issues: https://github.com/mit-crpg/openmc/issues
.. _license: http://openmc.readthedocs.io/en/stable/license.html
.. |licensebadge| image:: https://img.shields.io/github/license/mit-crpg/openmc.svg
:target: http://openmc.readthedocs.io/en/latest/license.html
:alt: License
.. |travisbadge| image:: https://travis-ci.org/mit-crpg/openmc.svg?branch=develop
:target: https://travis-ci.org/mit-crpg/openmc
:alt: Travis CI build status (Linux)
.. |coverallsbadge| image:: https://coveralls.io/repos/github/mit-crpg/openmc/badge.svg?branch=develop
:target: https://coveralls.io/github/mit-crpg/openmc?branch=develop
:alt: Code Coverage

2
setup.py
View file

@ -39,7 +39,7 @@ kwargs = {
'author': 'The OpenMC Development Team',
'author_email': 'openmc-dev@googlegroups.com',
'description': 'OpenMC',
'url': 'https://github.com/mit-crpg/openmc',
'url': 'https://github.com/openmc-dev/openmc',
'classifiers': [
'Development Status :: 4 - Beta',
'Intended Audience :: Developers',

4
src/api.F90
View file

@ -129,7 +129,7 @@ contains
index_ufs_mesh = -1
keff = ONE
legendre_to_tabular = .true.
legendre_to_tabular_points = 33
legendre_to_tabular_points = C_NONE
n_batch_interval = 1
n_lost_particles = 0
n_particles = 0
@ -310,7 +310,6 @@ contains
subroutine free_memory()
use cmfd_header
use mgxs_header
use plot_header
use sab_header
use settings
@ -328,7 +327,6 @@ contains
call free_memory_simulation()
call free_memory_nuclide()
call free_memory_settings()
call free_memory_mgxs()
call free_memory_sab()
call free_memory_source()
call free_memory_mesh()

4
src/cell.h
View file

@ -7,7 +7,7 @@
#include <vector>
#include "hdf5.h"
#include "pugixml/pugixml.hpp"
#include "pugixml.hpp"
namespace openmc {
@ -50,7 +50,7 @@ public:
//! A geometry primitive that links surfaces, universes, and materials
//==============================================================================
class Cell
class Cell
{
public:
int32_t id; //!< Unique ID

4
src/cmfd_input.F90
View file

@ -3,8 +3,8 @@ module cmfd_input
use, intrinsic :: ISO_C_BINDING
use cmfd_header
use mesh_header, only: mesh_dict
use mgxs_header, only: energy_bins
use mesh_header, only: mesh_dict
use mgxs_interface, only: energy_bins, num_energy_groups
use tally
use tally_header
use timer_header

19
src/constants.F90
View file

@ -405,6 +405,25 @@ module constants
DIFF_NUCLIDE_DENSITY = 2, &
DIFF_TEMPERATURE = 3
! Mgxs::get_xs enumerated types
integer(C_INT), parameter :: &
MG_GET_XS_TOTAL = 0, &
MG_GET_XS_ABSORPTION = 1, &
MG_GET_XS_INVERSE_VELOCITY = 2, &
MG_GET_XS_DECAY_RATE = 3, &
MG_GET_XS_SCATTER = 4, &
MG_GET_XS_SCATTER_MULT = 5, &
MG_GET_XS_SCATTER_FMU_MULT = 6, &
MG_GET_XS_SCATTER_FMU = 7, &
MG_GET_XS_FISSION = 8, &
MG_GET_XS_KAPPA_FISSION = 9, &
MG_GET_XS_PROMPT_NU_FISSION = 10, &
MG_GET_XS_DELAYED_NU_FISSION = 11, &
MG_GET_XS_NU_FISSION = 12, &
MG_GET_XS_CHI_PROMPT = 13, &
MG_GET_XS_CHI_DELAYED = 14
! ============================================================================
! RANDOM NUMBER STREAM CONSTANTS

77
src/constants.h
View file

@ -1,10 +1,85 @@
//! \file constants.h
//! A collection of constants
#ifndef CONSTANTS_H
#define CONSTANTS_H
#include <cmath>
#include <array>
#include <limits>
#include <vector>
namespace openmc {
namespace openmc{
typedef std::vector<double> double_1dvec;
typedef std::vector<std::vector<double> > double_2dvec;
typedef std::vector<std::vector<std::vector<double> > > double_3dvec;
typedef std::vector<std::vector<std::vector<std::vector<double> > > > double_4dvec;
typedef std::vector<std::vector<std::vector<std::vector<std::vector<double> > > > > double_5dvec;
typedef std::vector<std::vector<std::vector<std::vector<std::vector<std::vector<double> > > > > > double_6dvec;
typedef std::vector<int> int_1dvec;
typedef std::vector<std::vector<int> > int_2dvec;
typedef std::vector<std::vector<std::vector<int> > > int_3dvec;
constexpr int MAX_SAMPLE {10000};
constexpr std::array<int, 3> VERSION {0, 10, 0};
constexpr std::array<int, 2> VERSION_PARTICLE_RESTART {2, 0};
// Maximum number of words in a single line, length of line, and length of
// single word
constexpr int MAX_WORDS {500};
constexpr int MAX_LINE_LEN {250};
constexpr int MAX_WORD_LEN {150};
constexpr int MAX_FILE_LEN {255};
// Physical Constants
constexpr double K_BOLTZMANN {8.6173303e-5}; // Boltzmann constant in eV/K
// Angular distribution type
constexpr int ANGLE_ISOTROPIC {1};
constexpr int ANGLE_32_EQUI {2};
constexpr int ANGLE_TABULAR {3};
constexpr int ANGLE_LEGENDRE {4};
constexpr int ANGLE_HISTOGRAM {5};
// MGXS Table Types
constexpr int MGXS_ISOTROPIC {1}; // Isotroically weighted data
constexpr int MGXS_ANGLE {2}; // Data by angular bins
// Flag to denote this was a macroscopic data object
constexpr double MACROSCOPIC_AWR {-2.};
// Number of mu bins to use when converting Legendres to tabular type
constexpr int DEFAULT_NMU {33};
// Temperature treatment method
constexpr int TEMPERATURE_NEAREST {1};
constexpr int TEMPERATURE_INTERPOLATION {2};
// TODO: cmath::M_PI has 3 more digits precision than the Fortran constant we
// use so for now we will reuse the Fortran constant until we are OK with
// modifying test results
constexpr double PI {3.1415926535898};
const double SQRT_PI {std::sqrt(PI)};
// Mgxs::get_xs enumerated types
constexpr int MG_GET_XS_TOTAL {0};
constexpr int MG_GET_XS_ABSORPTION {1};
constexpr int MG_GET_XS_INVERSE_VELOCITY {2};
constexpr int MG_GET_XS_DECAY_RATE {3};
constexpr int MG_GET_XS_SCATTER {4};
constexpr int MG_GET_XS_SCATTER_MULT {5};
constexpr int MG_GET_XS_SCATTER_FMU_MULT {6};
constexpr int MG_GET_XS_SCATTER_FMU {7};
constexpr int MG_GET_XS_FISSION {8};
constexpr int MG_GET_XS_KAPPA_FISSION {9};
constexpr int MG_GET_XS_PROMPT_NU_FISSION {10};
constexpr int MG_GET_XS_DELAYED_NU_FISSION {11};
constexpr int MG_GET_XS_NU_FISSION {12};
constexpr int MG_GET_XS_CHI_PROMPT {13};
constexpr int MG_GET_XS_CHI_DELAYED {14};
extern "C" double FP_COINCIDENT;
extern "C" double FP_PRECISION;

216
src/doppler.F90
View file

@ -1,216 +0,0 @@
module doppler
use constants, only: ZERO, ONE, PI, K_BOLTZMANN
implicit none
real(8), parameter :: sqrt_pi_inv = ONE / sqrt(PI)
contains
!===============================================================================
! BROADEN takes a microscopic cross section at a temperature T_1 and Doppler
! broadens it to a higher temperature T_2 based on a method originally developed
! by Cullen and Weisbin (see "Exact Doppler Broadening of Tabulated Cross
! Sections," Nucl. Sci. Eng. 60, 199-229 (1976)). The only difference here is
! the F functions are evaluated based on complementary error functions rather
! than error functions as is done in the BROADR module of NJOY.
!===============================================================================
subroutine broaden(energy, xs, A_target, T, sigmaNew)
real(8), intent(in) :: energy(:) ! energy grid
real(8), intent(in) :: xs(:) ! unbroadened cross section
integer, intent(in) :: A_target ! mass number of target
real(8), intent(in) :: T ! temperature (difference)
real(8), intent(out) :: sigmaNew(:) ! broadened cross section
integer :: i, k ! loop indices
integer :: n ! number of energy points
real(8) :: F_a(0:4) ! F(a) functions as per C&W
real(8) :: F_b(0:4) ! F(b) functions as per C&W
real(8) :: H(0:4) ! H functions as per C&W
real(8), allocatable :: x(:) ! proportional to relative velocity
real(8) :: y ! proportional to neutron velocity
real(8) :: y_sq ! y**2
real(8) :: y_inv ! 1/y
real(8) :: y_inv_sq ! 1/y**2
real(8) :: alpha ! constant equal to A/kT
real(8) :: slope ! slope of xs between adjacent points
real(8) :: Ak, Bk ! coefficients at each point
real(8) :: a, b ! values of x(k)-y and x(k+1)-y
real(8) :: sigma ! broadened cross section at one point
! Determine alpha parameter -- have to convert k to eV/K
alpha = A_target/(K_BOLTZMANN * T)
! Allocate memory for x and assign values
n = size(energy)
allocate(x(n))
x = sqrt(alpha * energy)
! Loop over incoming neutron energies
ENERGY_NEUTRON: do i = 1, n
sigma = ZERO
y = x(i)
y_sq = y*y
y_inv = ONE / y
y_inv_sq = y_inv / y
! =======================================================================
! EVALUATE FIRST TERM FROM x(k) - y = 0 to -4
k = i
a = ZERO
call calculate_F(F_a, a)
do while (a >= -4.0 .and. k > 1)
! Move to next point
F_b = F_a
k = k - 1
a = x(k) - y
! Calculate F and H functions
call calculate_F(F_a, a)
H = F_a - F_b
! Calculate A(k), B(k), and slope terms
Ak = y_inv_sq*H(2) + 2.0*y_inv*H(1) + H(0)
Bk = y_inv_sq*H(4) + 4.0*y_inv*H(3) + 6.0*H(2) + 4.0*y*H(1) + y_sq*H(0)
slope = (xs(k+1) - xs(k)) / (x(k+1)**2 - x(k)**2)
! Add contribution to broadened cross section
sigma = sigma + Ak*(xs(k) - slope*x(k)**2) + slope*Bk
end do
! =======================================================================
! EXTEND CROSS SECTION TO 0 ASSUMING 1/V SHAPE
if (k == 1 .and. a >= -4.0) then
! Since x = 0, this implies that a = -y
F_b = F_a
a = -y
! Calculate F and H functions
call calculate_F(F_a, a)
H = F_a - F_b
! Add contribution to broadened cross section
sigma = sigma + xs(k)*x(k)*(y_inv_sq*H(1) + y_inv*H(0))
end if
! =======================================================================
! EVALUATE FIRST TERM FROM x(k) - y = 0 to 4
k = i
b = ZERO
call calculate_F(F_b, b)
do while (b <= 4.0 .and. k < n)
! Move to next point
F_a = F_b
k = k + 1
b = x(k) - y
! Calculate F and H functions
call calculate_F(F_b, b)
H = F_a - F_b
! Calculate A(k), B(k), and slope terms
Ak = y_inv_sq*H(2) + 2.0*y_inv*H(1) + H(0)
Bk = y_inv_sq*H(4) + 4.0*y_inv*H(3) + 6.0*H(2) + 4.0*y*H(1) + y_sq*H(0)
slope = (xs(k) - xs(k-1)) / (x(k)**2 - x(k-1)**2)
! Add contribution to broadened cross section
sigma = sigma + Ak*(xs(k) - slope*x(k)**2) + slope*Bk
end do
! =======================================================================
! EXTEND CROSS SECTION TO INFINITY ASSUMING CONSTANT SHAPE
if (k == n .and. b <= 4.0) then
! Calculate F function at last energy point
a = x(k) - y
call calculate_F(F_a, a)
! Add contribution to broadened cross section
sigma = sigma + xs(k) * (y_inv_sq*F_a(2) + 2.0*y_inv*F_a(1) + F_a(0))
end if
! =======================================================================
! EVALUATE SECOND TERM FROM x(k) + y = 0 to +4
if (y <= 4.0) then
! Swap signs on y
y = -y
y_inv = -y_inv
k = 1
! Calculate a and b based on 0 and x(1)
a = -y
b = x(k) - y
! Calculate F and H functions
call calculate_F(F_a, a)
call calculate_F(F_b, b)
H = F_a - F_b
! Add contribution to broadened cross section
sigma = sigma - xs(k) * x(k) * (y_inv_sq*H(1) + y_inv*H(0))
! Now progress forward doing the remainder of the second term
do while (b <= 4.0)
! Move to next point
F_a = F_b
k = k + 1
b = x(k) - y
! Calculate F and H functions
call calculate_F(F_b, b)
H = F_a - F_b
! Calculate A(k), B(k), and slope terms
Ak = y_inv_sq*H(2) + 2.0*y_inv*H(1) + H(0)
Bk = y_inv_sq*H(4) + 4.0*y_inv*H(3) + 6.0*H(2) + 4.0*y*H(1) &
+ y_sq*H(0)
slope = (xs(k) - xs(k-1)) / (x(k)**2 - x(k-1)**2)
! Add contribution to broadened cross section
sigma = sigma - Ak*(xs(k) - slope*x(k)**2) - slope*Bk
end do
end if
! Set broadened cross section
sigmaNew(i) = sigma
end do ENERGY_NEUTRON
end subroutine broaden
!===============================================================================
! CALCULATE_F evaluates the function:
!
! F(n,a) = 1/sqrt(pi)*int(z^n*exp(-z^2), z = a to infinity)
!
! The five values returned in a vector correspond to the integral for n = 0
! through 4. These functions are called over and over during the Doppler
! broadening routine.
!===============================================================================
subroutine calculate_F(F, a)
real(8), intent(inout) :: F(0:4)
real(8), intent(in) :: a
#ifndef NO_F2008
F(0) = 0.5*erfc(a)
#endif
F(1) = 0.5*sqrt_pi_inv*exp(-a*a)
F(2) = 0.5*F(0) + a*F(1)
F(3) = F(1)*(1.0 + a*a)
F(4) = 0.75*F(0) + F(1)*a*(1.5 + a*a)
end subroutine calculate_F
end module doppler

56
src/endf_header.F90
View file

@ -37,7 +37,6 @@ module endf_header
contains
procedure :: from_hdf5 => polynomial_from_hdf5
procedure :: evaluate => polynomial_evaluate
procedure :: from_ace => polynomial_from_ace
end type Polynomial
!===============================================================================
@ -52,7 +51,6 @@ module endf_header
real(8), allocatable :: x(:) ! values of abscissa
real(8), allocatable :: y(:) ! values of ordinate
contains
procedure :: from_ace => tabulated1d_from_ace
procedure :: from_hdf5 => tabulated1d_from_hdf5
procedure :: evaluate => tabulated1d_evaluate
end type Tabulated1D
@ -63,24 +61,6 @@ contains
! Polynomial implementation
!===============================================================================
subroutine polynomial_from_ace(this, xss, idx)
class(Polynomial), intent(inout) :: this
real(8), intent(in) :: xss(:)
integer, intent(in) :: idx
integer :: nc ! number of coefficients (order - 1)
! Clear space
if (allocated(this % coef)) deallocate(this % coef)
! Determine number of coefficients
nc = nint(xss(idx))
! Allocate space for and read coefficients
allocate(this % coef(nc))
this % coef(:) = xss(idx + 1 : idx + nc)
end subroutine polynomial_from_ace
subroutine polynomial_from_hdf5(this, dset_id)
class(Polynomial), intent(inout) :: this
integer(HID_T), intent(in) :: dset_id
@ -111,42 +91,6 @@ contains
! Tabulated1D implementation
!===============================================================================
subroutine tabulated1d_from_ace(this, xss, idx)
class(Tabulated1D), intent(inout) :: this
real(8), intent(in) :: xss(:)
integer, intent(in) :: idx
integer :: nr, ne
! Clear space
if (allocated(this % nbt)) deallocate(this % nbt)
if (allocated(this % int)) deallocate(this % int)
if (allocated(this % x)) deallocate(this % x)
if (allocated(this % y)) deallocate(this % y)
! Determine number of regions
nr = nint(xss(idx))
this % n_regions = nr
! 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 tabulated1d_from_ace
subroutine tabulated1d_from_hdf5(this, dset_id)
class(Tabulated1D), intent(inout) :: this
integer(HID_T), intent(in) :: dset_id

9
src/error.F90
View file

@ -265,4 +265,13 @@ contains
end subroutine write_message
subroutine write_message_from_c(message, message_len, level) bind(C)
integer(C_INT), intent(in), value :: message_len
character(kind=C_CHAR), intent(in) :: message(message_len)
integer(C_INT), intent(in), value :: level
character(message_len+1) :: message_out
write(message_out, *) message
call write_message(message_out, level)
end subroutine write_message_from_c
end module error

21
src/error.h
View file

@ -11,7 +11,8 @@ namespace openmc {
extern "C" void fatal_error_from_c(const char* message, int message_len);
extern "C" void warning_from_c(const char* message, int message_len);
extern "C" void write_message_from_c(const char* message, int message_len,
int level);
inline
void fatal_error(const char *message)
@ -43,5 +44,23 @@ void warning(const std::stringstream& message)
warning(message.str());
}
inline
void write_message(const char* message, int level)
{
write_message_from_c(message, strlen(message), level);
}
inline
void write_message(const std::string& message, int level)
{
write_message_from_c(message.c_str(), message.length(), level);
}
inline
void write_message(const std::stringstream& message, int level)
{
write_message(message.str(), level);
}
} // namespace openmc
#endif // ERROR_H

2
src/geometry_header.F90
View file

@ -736,7 +736,7 @@ contains
! find which material is associated with this cell (material_index
! is the index into the materials array)
if (present(instance)) then
material_index = cells(index) % material(instance)
material_index = cells(index) % material(instance + 1)
else
material_index = cells(index) % material(1)
end if

166
src/hdf5_interface.cpp
View file

@ -447,6 +447,172 @@ read_complex(hid_t obj_id, const char* name, double _Complex* buffer, bool indep
}
void
read_nd_vector(hid_t obj_id, const char* name, std::vector<double>& result,
bool must_have)
{
if (object_exists(obj_id, name)) {
read_double(obj_id, name, result.data(), true);
} else if (must_have) {
fatal_error(std::string("Must provide " + std::string(name) + "!"));
}
}
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<double> >& result, bool must_have)
{
if (object_exists(obj_id, name)) {
int dim1 = result.size();
int dim2 = result[0].size();
double temp_arr[dim1 * dim2];
read_double(obj_id, name, temp_arr, true);
int temp_idx = 0;
for (int i = 0; i < dim1; i++) {
for (int j = 0; j < dim2; j++) {
result[i][j] = temp_arr[temp_idx++];
}
}
} else if (must_have) {
fatal_error(std::string("Must provide " + std::string(name) + "!"));
}
}
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<int> >& result, bool must_have)
{
if (object_exists(obj_id, name)) {
int dim1 = result.size();
int dim2 = result[0].size();
int temp_arr[dim1 * dim2];
read_int(obj_id, name, temp_arr, true);
int temp_idx = 0;
for (int i = 0; i < dim1; i++) {
for (int j = 0; j < dim2; j++) {
result[i][j] = temp_arr[temp_idx++];
}
}
} else if (must_have) {
fatal_error(std::string("Must provide " + std::string(name) + "!"));
}
}
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<std::vector<double> > >& result,
bool must_have)
{
if (object_exists(obj_id, name)) {
int dim1 = result.size();
int dim2 = result[0].size();
int dim3 = result[0][0].size();
double temp_arr[dim1 * dim2 * dim3];
read_double(obj_id, name, temp_arr, true);
int temp_idx = 0;
for (int i = 0; i < dim1; i++) {
for (int j = 0; j < dim2; j++) {
for (int k = 0; k < dim3; k++) {
result[i][j][k] = temp_arr[temp_idx++];
}
}
}
} else if (must_have) {
fatal_error(std::string("Must provide " + std::string(name) + "!"));
}
}
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<std::vector<int> > >& result,
bool must_have)
{
if (object_exists(obj_id, name)) {
int dim1 = result.size();
int dim2 = result[0].size();
int dim3 = result[0][0].size();
int temp_arr[dim1 * dim2 * dim3];
read_int(obj_id, name, temp_arr, true);
int temp_idx = 0;
for (int i = 0; i < dim1; i++) {
for (int j = 0; j < dim2; j++) {
for (int k = 0; k < dim3; k++) {
result[i][j][k] = temp_arr[temp_idx++];
}
}
}
} else if (must_have) {
fatal_error(std::string("Must provide " + std::string(name) + "!"));
}
}
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<std::vector<std::vector<double> > > >& result,
bool must_have)
{
if (object_exists(obj_id, name)) {
int dim1 = result.size();
int dim2 = result[0].size();
int dim3 = result[0][0].size();
int dim4 = result[0][0][0].size();
double temp_arr[dim1 * dim2 * dim3 * dim4];
read_double(obj_id, name, temp_arr, true);
int temp_idx = 0;
for (int i = 0; i < dim1; i++) {
for (int j = 0; j < dim2; j++) {
for (int k = 0; k < dim3; k++) {
for (int l = 0; l < dim4; l++) {
result[i][j][k][l] = temp_arr[temp_idx++];
}
}
}
}
} else if (must_have) {
fatal_error(std::string("Must provide " + std::string(name) + "!"));
}
}
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<std::vector<std::vector<std::vector<double> > > > >& result,
bool must_have)
{
if (object_exists(obj_id, name)) {
int dim1 = result.size();
int dim2 = result[0].size();
int dim3 = result[0][0].size();
int dim4 = result[0][0][0].size();
int dim5 = result[0][0][0][0].size();
double temp_arr[dim1 * dim2 * dim3 * dim4 * dim5];
read_double(obj_id, name, temp_arr, true);
int temp_idx = 0;
for (int i = 0; i < dim1; i++) {
for (int j = 0; j < dim2; j++) {
for (int k = 0; k < dim3; k++) {
for (int l = 0; l < dim4; l++) {
for (int m = 0; m < dim5; m++) {
result[i][j][k][l][m] = temp_arr[temp_idx++];
}
}
}
}
}
} else if (must_have) {
fatal_error(std::string("Must provide " + std::string(name) + "!"));
}
}
void
read_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score, double* results)
{

34
src/hdf5_interface.h
View file

@ -7,6 +7,7 @@
#include <array>
#include <string>
#include <sstream>
#include <vector>
#include <complex.h>
@ -55,6 +56,39 @@ extern "C" void read_string(hid_t obj_id, const char* name, size_t slen,
extern "C" void read_complex(hid_t obj_id, const char* name,
double _Complex* buffer, bool indep);
void
read_nd_vector(hid_t obj_id, const char* name, std::vector<double>& result,
bool must_have = false);
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<double> >& result,
bool must_have = false);
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<int> >& result, bool must_have = false);
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<std::vector<double> > >& result,
bool must_have = false);
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<std::vector<int> > >& result,
bool must_have = false);
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<std::vector<std::vector<double> > > >& result,
bool must_have = false);
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<std::vector<std::vector<std::vector<double> > > > >& result,
bool must_have = false);
extern "C" void read_tally_results(hid_t group_id, hsize_t n_filter,
hsize_t n_score, double* results);

12
src/input_xml.F90
View file

@ -19,7 +19,7 @@ module input_xml
use mesh_header
use message_passing
use mgxs_data, only: create_macro_xs, read_mgxs
use mgxs_header
use mgxs_interface
use nuclide_header
use output, only: title, header, print_plot
use plot_header
@ -228,7 +228,7 @@ contains
&not exist! In order to run OpenMC, you first need a set of input &
&files; at a minimum, this includes settings.xml, geometry.xml, &
&and materials.xml. Please consult the user's guide at &
&http://mit-crpg.github.io/openmc for further information.")
&http://openmc.readthedocs.io for further information.")
else
! The settings.xml file is optional if we just want to make a plot.
return
@ -1369,7 +1369,7 @@ contains
&materials.xml, settings.xml, or in the OPENMC_CROSS_SECTIONS&
& environment variable. OpenMC needs such a file to identify &
&where to find ACE cross section libraries. Please consult the&
& user's guide at http://mit-crpg.github.io/openmc for &
& user's guide at http://openmc.readthedocs.io for &
&information on how to set up ACE cross section libraries.")
else
call warning("The CROSS_SECTIONS environment variable is &
@ -1387,7 +1387,7 @@ contains
&materials.xml or in the OPENMC_MG_CROSS_SECTIONS environment &
&variable. OpenMC needs such a file to identify where to &
&find MG cross section libraries. Please consult the user's &
&guide at http://mit-crpg.github.io/openmc for information on &
&guide at http://openmc.readthedocs.io for information on &
&how to set up MG cross section libraries.")
else if (len_trim(env_variable) /= 0) then
path_cross_sections = trim(env_variable)
@ -3458,7 +3458,7 @@ contains
if (run_CE) then
awr = nuclides(mat % nuclide(j)) % awr
else
awr = nuclides_MG(mat % nuclide(j)) % obj % awr
awr = get_awr_c(mat % nuclide(j))
end if
! if given weight percent, convert all values so that they are divided
@ -3483,7 +3483,7 @@ contains
if (run_CE) then
awr = nuclides(mat % nuclide(j)) % awr
else
awr = nuclides_MG(mat % nuclide(j)) % obj % awr
awr = get_awr_c(mat % nuclide(j))
end if
x = mat % atom_density(j)
sum_percent = sum_percent + x*awr

2
src/lattice.h
View file

@ -9,7 +9,7 @@
#include "constants.h"
#include "hdf5.h"
#include "pugixml/pugixml.hpp"
#include "pugixml.hpp"
namespace openmc {

2
src/material_header.F90
View file

@ -34,7 +34,7 @@ module material_header
integer :: n_nuclides = 0 ! number of nuclides
integer, allocatable :: nuclide(:) ! index in nuclides array
real(8) :: density ! total atom density in atom/b-cm
real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm
real(C_DOUBLE), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm
real(8) :: density_gpcc ! total density in g/cm^3
! To improve performance of tallying, we store an array (direct address

2
src/math_functions.cpp
View file

@ -97,7 +97,7 @@ void calc_pn_c(int n, double x, double pnx[]) {
}
// Use recursion relation to build the higher orders
for (int l = 1; l < n; l ++) {
for (int l = 1; l < n; l++) {
pnx[l + 1] = ((2 * l + 1) * x * pnx[l] - l * pnx[l - 1]) / (l + 1);
}
}

12
src/math_functions.h
View file

@ -7,22 +7,12 @@
#include <cmath>
#include <cstdlib>
#include "constants.h"
#include "random_lcg.h"
namespace openmc {
//==============================================================================
// Module constants.
//==============================================================================
// TODO: cmath::M_PI has 3 more digits precision than the Fortran constant we
// use so for now we will reuse the Fortran constant until we are OK with
// modifying test results
extern "C" constexpr double PI {3.1415926535898};
extern "C" const double SQRT_PI {std::sqrt(PI)};
//==============================================================================
//! Calculate the percentile of the standard normal distribution with a
//! specified probability level.

712
src/mgxs.cpp Normal file
View file

@ -0,0 +1,712 @@
#include <cmath>
#include <cstdlib>
#include <algorithm>
#include <valarray>
#ifdef _OPENMP
# include <omp.h>
#endif
#include "error.h"
#include "math_functions.h"
#include "random_lcg.h"
#include "string_functions.h"
#include "mgxs.h"
namespace openmc {
// Storage for the MGXS data
std::vector<Mgxs> nuclides_MG;
std::vector<Mgxs> macro_xs;
//==============================================================================
// Mgxs base-class methods
//==============================================================================
void
Mgxs::init(const std::string& in_name, double in_awr,
const double_1dvec& in_kTs, bool in_fissionable, int in_scatter_format,
int in_num_groups, int in_num_delayed_groups, bool in_is_isotropic,
const double_1dvec& in_polar, const double_1dvec& in_azimuthal)
{
// Set the metadata
name = in_name;
awr = in_awr;
kTs = in_kTs;
fissionable = in_fissionable;
scatter_format = in_scatter_format;
num_groups = in_num_groups;
num_delayed_groups = in_num_delayed_groups;
xs.resize(in_kTs.size());
is_isotropic = in_is_isotropic;
n_pol = in_polar.size();
n_azi = in_azimuthal.size();
polar = in_polar;
azimuthal = in_azimuthal;
// Set the cross section index cache
#ifdef _OPENMP
int n_threads = omp_get_max_threads();
#else
int n_threads = 1;
#endif
cache.resize(n_threads);
// std::vector.resize() will value-initialize the members of cache[:]
}
//==============================================================================
void
Mgxs::metadata_from_hdf5(hid_t xs_id, int in_num_groups,
int in_num_delayed_groups, const double_1dvec& temperature,
double tolerance, int_1dvec& temps_to_read, int& order_dim, int& method)
{
// get name
char char_name[MAX_WORD_LEN];
get_name(xs_id, char_name);
std::string in_name {char_name};
// remove the leading '/'
in_name = in_name.substr(1);
// Get the AWR
double in_awr;
if (attribute_exists(xs_id, "atomic_weight_ratio")) {
read_attr_double(xs_id, "atomic_weight_ratio", &in_awr);
} else {
in_awr = MACROSCOPIC_AWR;
}
// Determine the available temperatures
hid_t kT_group = open_group(xs_id, "kTs");
int num_temps = get_num_datasets(kT_group);
char* dset_names[num_temps];
for (int i = 0; i < num_temps; i++) {
dset_names[i] = new char[151];
}
get_datasets(kT_group, dset_names);
double_1dvec available_temps(num_temps);
for (int i = 0; i < num_temps; i++) {
read_double(kT_group, dset_names[i], &available_temps[i], true);
// convert eV to Kelvin
available_temps[i] /= K_BOLTZMANN;
// Done with dset_names, so delete it
delete[] dset_names[i];
}
std::sort(available_temps.begin(), available_temps.end());
// If only one temperature is available, lets just use nearest temperature
// interpolation
if ((num_temps == 1) && (method == TEMPERATURE_INTERPOLATION)) {
warning("Cross sections for " + strtrim(name) + " are only available " +
"at one temperature. Reverting to the nearest temperature " +
"method.");
method = TEMPERATURE_NEAREST;
}
switch(method) {
case TEMPERATURE_NEAREST:
// Find the minimum difference
for (int i = 0; i < temperature.size(); i++) {
std::valarray<double> temp_diff(available_temps.data(),
available_temps.size());
temp_diff = std::abs(temp_diff - temperature[i]);
int i_closest = std::min_element(std::begin(temp_diff), std::end(temp_diff)) -
std::begin(temp_diff);
double temp_actual = available_temps[i_closest];
if (std::abs(temp_actual - temperature[i]) < tolerance) {
if (std::find(temps_to_read.begin(), temps_to_read.end(),
std::round(temp_actual)) == temps_to_read.end()) {
temps_to_read.push_back(std::round(temp_actual));
} else {
fatal_error("MGXS Library does not contain cross section for " +
in_name + " at or near " +
std::to_string(std::round(temperature[i])) + " K.");
}
}
}
break;
case TEMPERATURE_INTERPOLATION:
for (int i = 0; i < temperature.size(); i++) {
for (int j = 0; j < num_temps - 1; j++) {
if ((available_temps[j] <= temperature[i]) &&
(temperature[i] < available_temps[j + 1])) {
if (std::find(temps_to_read.begin(),
temps_to_read.end(),
std::round(available_temps[j])) == temps_to_read.end()) {
temps_to_read.push_back(std::round((int)available_temps[j]));
}
if (std::find(temps_to_read.begin(), temps_to_read.end(),
std::round(available_temps[j + 1])) == temps_to_read.end()) {
temps_to_read.push_back(std::round((int) available_temps[j + 1]));
}
continue;
}
}
fatal_error("MGXS Library does not contain cross sections for " +
in_name + " at temperatures that bound " +
std::to_string(std::round(temperature[i])));
}
}
std::sort(temps_to_read.begin(), temps_to_read.end());
// Get the library's temperatures
int n_temperature = temps_to_read.size();
double_1dvec in_kTs(n_temperature);
for (int i = 0; i < n_temperature; i++) {
std::string temp_str(std::to_string(temps_to_read[i]) + "K");
//read exact temperature value
read_double(kT_group, temp_str.c_str(), &in_kTs[i], true);
}
close_group(kT_group);
// Load the remaining metadata
int in_scatter_format;
if (attribute_exists(xs_id, "scatter_format")) {
std::string temp_str(MAX_WORD_LEN, ' ');
read_attr_string(xs_id, "scatter_format", MAX_WORD_LEN, &temp_str[0]);
to_lower(strtrim(temp_str));
if (temp_str.compare(0, 8, "legendre") == 0) {
in_scatter_format = ANGLE_LEGENDRE;
} else if (temp_str.compare(0, 9, "histogram") == 0) {
in_scatter_format = ANGLE_HISTOGRAM;
} else if (temp_str.compare(0, 7, "tabular") == 0) {
in_scatter_format = ANGLE_TABULAR;
} else {
fatal_error("Invalid scatter_format option!");
}
} else {
in_scatter_format = ANGLE_LEGENDRE;
}
if (attribute_exists(xs_id, "scatter_shape")) {
std::string temp_str(MAX_WORD_LEN, ' ');
read_attr_string(xs_id, "scatter_shape", MAX_WORD_LEN, &temp_str[0]);
to_lower(strtrim(temp_str));
if (temp_str.compare(0, 14, "[g][g\'][order]") != 0) {
fatal_error("Invalid scatter_shape option!");
}
}
bool in_fissionable = false;
if (attribute_exists(xs_id, "fissionable")) {
int int_fiss;
read_attr_int(xs_id, "fissionable", &int_fiss);
in_fissionable = int_fiss;
} else {
fatal_error("Fissionable element must be set!");
}
// Get the library's value for the order
if (attribute_exists(xs_id, "order")) {
read_attr_int(xs_id, "order", &order_dim);
} else {
fatal_error("Order must be provided!");
}
// Store the dimensionality of the data in order_dim.
// For Legendre data, we usually refer to it as Pn where n is the order.
// However Pn has n+1 sets of points (since you need to count the P0
// moment). Adjust for that. Histogram and Tabular formats dont need this
// adjustment.
if (in_scatter_format == ANGLE_LEGENDRE) {
++order_dim;
}
// Get the angular information
int in_n_pol;
int in_n_azi;
bool in_is_isotropic = true;
if (attribute_exists(xs_id, "representation")) {
std::string temp_str(MAX_WORD_LEN, ' ');
read_attr_string(xs_id, "representation", MAX_WORD_LEN, &temp_str[0]);
to_lower(strtrim(temp_str));
if (temp_str.compare(0, 5, "angle") == 0) {
in_is_isotropic = false;
} else if (temp_str.compare(0, 9, "isotropic") != 0) {
fatal_error("Invalid Data Representation!");
}
}
if (!in_is_isotropic) {
if (attribute_exists(xs_id, "num_polar")) {
read_attr_int(xs_id, "num_polar", &in_n_pol);
} else {
fatal_error("num_polar must be provided!");
}
if (attribute_exists(xs_id, "num_azimuthal")) {
read_attr_int(xs_id, "num_azimuthal", &in_n_azi);
} else {
fatal_error("num_azimuthal must be provided!");
}
} else {
in_n_pol = 1;
in_n_azi = 1;
}
// Set the angular bins to use equally-spaced bins
double_1dvec in_polar(in_n_pol);
double dangle = PI / in_n_pol;
for (int p = 0; p < in_n_pol; p++) {
in_polar[p] = (p + 0.5) * dangle;
}
double_1dvec in_azimuthal(in_n_azi);
dangle = 2. * PI / in_n_azi;
for (int a = 0; a < in_n_azi; a++) {
in_azimuthal[a] = (a + 0.5) * dangle - PI;
}
// Finally use this data to initialize the MGXS Object
init(in_name, in_awr, in_kTs, in_fissionable, in_scatter_format,
in_num_groups, in_num_delayed_groups, in_is_isotropic, in_polar,
in_azimuthal);
}
//==============================================================================
Mgxs::Mgxs(hid_t xs_id, int energy_groups, int delayed_groups,
const double_1dvec& temperature, double tolerance, int max_order,
bool legendre_to_tabular, int legendre_to_tabular_points, int& method)
{
// Call generic data gathering routine (will populate the metadata)
int order_data;
int_1dvec temps_to_read;
metadata_from_hdf5(xs_id, energy_groups, delayed_groups, temperature,
tolerance, temps_to_read, order_data, method);
// Set number of energy and delayed groups
int final_scatter_format = scatter_format;
if (legendre_to_tabular) {
if (scatter_format == ANGLE_LEGENDRE) final_scatter_format = ANGLE_TABULAR;
}
// Load the more specific XsData information
for (int t = 0; t < temps_to_read.size(); t++) {
xs[t] = XsData(energy_groups, delayed_groups, fissionable,
final_scatter_format, n_pol, n_azi);
// Get the temperature as a string and then open the HDF5 group
std::string temp_str = std::to_string(temps_to_read[t]) + "K";
hid_t xsdata_grp = open_group(xs_id, temp_str.c_str());
xs[t].from_hdf5(xsdata_grp, fissionable, scatter_format,
final_scatter_format, order_data, max_order,
legendre_to_tabular_points, is_isotropic, n_pol, n_azi);
close_group(xsdata_grp);
} // end temperature loop
// Make sure the scattering format is updated to the final case
scatter_format = final_scatter_format;
}
//==============================================================================
Mgxs::Mgxs(const std::string& in_name, const double_1dvec& mat_kTs,
const std::vector<Mgxs*>& micros, const double_1dvec& atom_densities,
double tolerance, int& method)
{
// Get the minimum data needed to initialize:
// Dont need awr, but lets just initialize it anyways
double in_awr = -1.;
// start with the assumption it is not fissionable
bool in_fissionable = false;
for (int m = 0; m < micros.size(); m++) {
if (micros[m]->fissionable) in_fissionable = true;
}
// Force all of the following data to be the same; these will be verified
// to be true later
int in_scatter_format = micros[0]->scatter_format;
int in_num_groups = micros[0]->num_groups;
int in_num_delayed_groups = micros[0]->num_delayed_groups;
bool in_is_isotropic = micros[0]->is_isotropic;
double_1dvec in_polar = micros[0]->polar;
double_1dvec in_azimuthal = micros[0]->azimuthal;
init(in_name, in_awr, mat_kTs, in_fissionable, in_scatter_format,
in_num_groups, in_num_delayed_groups, in_is_isotropic, in_polar,
in_azimuthal);
// Create the xs data for each temperature
for (int t = 0; t < mat_kTs.size(); t++) {
xs[t] = XsData(in_num_groups, in_num_delayed_groups, in_fissionable,
in_scatter_format, in_polar.size(), in_azimuthal.size());
// Find the right temperature index to use
double temp_desired = mat_kTs[t];
// Create the list of temperature indices and interpolation factors for
// each microscopic data at the material temperature
int_1dvec micro_t(micros.size(), 0);
double_1dvec micro_t_interp(micros.size(), 0.);
for (int m = 0; m < micros.size(); m++) {
switch(method) {
case TEMPERATURE_NEAREST:
{
// Find the nearest temperature
std::valarray<double> temp_diff(micros[m]->kTs.data(),
micros[m]->kTs.size());
temp_diff = std::abs(temp_diff - temp_desired);
micro_t[m] = std::min_element(std::begin(temp_diff),
std::end(temp_diff)) -
std::begin(temp_diff);
double temp_actual = micros[m]->kTs[micro_t[m]];
if (std::abs(temp_actual - temp_desired) >= K_BOLTZMANN * tolerance) {
fatal_error("MGXS Library does not contain cross section for " +
name + " at or near " +
std::to_string(std::round(temp_desired / K_BOLTZMANN))
+ " K.");
}
}
break;
case TEMPERATURE_INTERPOLATION:
// Get a list of bounding temperatures for each actual temperature
// present in the model
for (int k = 0; k < micros[m]->kTs.size() - 1; k++) {
if ((micros[m]->kTs[k] <= temp_desired) &&
(temp_desired < micros[m]->kTs[k + 1])) {
micro_t[m] = k;
if (k == 0) {
micro_t_interp[m] = (temp_desired - micros[m]->kTs[k]) /
(micros[m]->kTs[k + 1] - micros[m]->kTs[k]);
} else {
micro_t_interp[m] = 1.;
}
}
}
} // end switch
} // end microscopic temperature loop
// We are about to loop through each of the microscopic objects
// and incorporate the contribution of each microscopic data at
// one of the two temperature interpolants to this macroscopic quantity.
// If we are doing nearest temperature interpolation, then we don't need
// to do the 2nd temperature
int num_interp_points = 2;
if (method == TEMPERATURE_NEAREST) num_interp_points = 1;
for (int interp_point = 0; interp_point < num_interp_points; interp_point++) {
double_1dvec interp(micros.size());
double_1dvec temp_indices(micros.size());
for (int m = 0; m < micros.size(); m++) {
interp[m] = (1. - micro_t_interp[m]) * atom_densities[m];
temp_indices[m] = micro_t[m] + interp_point;
}
combine(micros, interp, micro_t, t);
} // end loop to sum all micros across the temperatures
} // end temperature (t) loop
}
//==============================================================================
void
Mgxs::combine(const std::vector<Mgxs*>& micros, const double_1dvec& scalars,
const int_1dvec& micro_ts, int this_t)
{
// Build the vector of pointers to the xs objects within micros
std::vector<XsData*> those_xs(micros.size());
for (int i = 0; i < micros.size(); i++) {
if (!xs[this_t].equiv(micros[i]->xs[micro_ts[i]])) {
fatal_error("Cannot combine the Mgxs objects!");
}
those_xs[i] = &(micros[i]->xs[micro_ts[i]]);
}
xs[this_t].combine(those_xs, scalars);
}
//==============================================================================
double
Mgxs::get_xs(int xstype, int gin, int* gout, double* mu, int* dg)
{
// This method assumes that the temperature and angle indices are set
#ifdef _OPENMP
int tid = omp_get_thread_num();
XsData* xs_t = &xs[cache[tid].t];
int a = cache[tid].a;
#else
XsData* xs_t = &xs[cache[0].t];
int a = cache[0].a;
#endif
double val;
switch(xstype) {
case MG_GET_XS_TOTAL:
val = xs_t->total[a][gin];
break;
case MG_GET_XS_NU_FISSION:
val = fissionable ? xs_t->nu_fission[a][gin] : 0.;
break;
case MG_GET_XS_ABSORPTION:
val = xs_t->absorption[a][gin];
break;
case MG_GET_XS_FISSION:
val = fissionable ? xs_t->fission[a][gin] : 0.;
break;
case MG_GET_XS_KAPPA_FISSION:
val = fissionable ? xs_t->kappa_fission[a][gin] : 0.;
break;
case MG_GET_XS_SCATTER:
case MG_GET_XS_SCATTER_MULT:
case MG_GET_XS_SCATTER_FMU_MULT:
case MG_GET_XS_SCATTER_FMU:
val = xs_t->scatter[a]->get_xs(xstype, gin, gout, mu);
break;
case MG_GET_XS_PROMPT_NU_FISSION:
val = fissionable ? xs_t->prompt_nu_fission[a][gin] : 0.;
break;
case MG_GET_XS_DELAYED_NU_FISSION:
if (fissionable) {
if (dg != nullptr) {
val = xs_t->delayed_nu_fission[a][gin][*dg];
} else {
val = 0.;
for (auto& num : xs_t->delayed_nu_fission[a][gin]) {
val += num;
}
}
} else {
val = 0.;
}
break;
case MG_GET_XS_CHI_PROMPT:
if (fissionable) {
if (gout != nullptr) {
val = xs_t->chi_prompt[a][gin][*gout];
} else {
// provide an outgoing group-wise sum
val = 0.;
for (auto& num : xs_t->chi_prompt[a][gin]) {
val += num;
}
}
} else {
val = 0.;
}
break;
case MG_GET_XS_CHI_DELAYED:
if (fissionable) {
if (gout != nullptr) {
if (dg != nullptr) {
val = xs_t->chi_delayed[a][gin][*gout][*dg];
} else {
val = xs_t->chi_delayed[a][gin][*gout][0];
}
} else {
if (dg != nullptr) {
val = 0.;
for (int i = 0; i < xs_t->chi_delayed[a][gin].size(); i++) {
val += xs_t->chi_delayed[a][gin][i][*dg];
}
} else {
val = 0.;
for (int i = 0; i < xs_t->chi_delayed[a][gin].size(); i++) {
for (auto& num : xs_t->chi_delayed[a][gin][i]) {
val += num;
}
}
}
}
} else {
val = 0.;
}
break;
case MG_GET_XS_INVERSE_VELOCITY:
val = xs_t->inverse_velocity[a][gin];
break;
case MG_GET_XS_DECAY_RATE:
if (dg != nullptr) {
val = xs_t->decay_rate[a][*dg + 1];
} else {
val = xs_t->decay_rate[a][0];
}
break;
default:
val = 0.;
}
return val;
}
//==============================================================================
void
Mgxs::sample_fission_energy(int gin, int& dg, int& gout)
{
// This method assumes that the temperature and angle indices are set
#ifdef _OPENMP
int tid = omp_get_thread_num();
#else
int tid = 0;
#endif
XsData* xs_t = &xs[cache[tid].t];
double nu_fission = xs_t->nu_fission[cache[tid].a][gin];
// Find the probability of having a prompt neutron
double prob_prompt =
xs_t->prompt_nu_fission[cache[tid].a][gin];
// sample random numbers
double xi_pd = prn() * nu_fission;
double xi_gout = prn();
// Select whether the neutron is prompt or delayed
if (xi_pd <= prob_prompt) {
// the neutron is prompt
// set the delayed group for the particle to be -1, indicating prompt
dg = -1;
// sample the outgoing energy group
gout = 0;
double prob_gout =
xs_t->chi_prompt[cache[tid].a][gin][gout];
while (prob_gout < xi_gout) {
gout++;
prob_gout += xs_t->chi_prompt[cache[tid].a][gin][gout];
}
} else {
// the neutron is delayed
// get the delayed group
dg = 0;
while (xi_pd >= prob_prompt) {
dg++;
prob_prompt +=
xs_t->delayed_nu_fission[cache[tid].a][gin][dg];
}
// adjust dg in case of round-off error
dg = std::min(dg, num_delayed_groups - 1);
// sample the outgoing energy group
gout = 0;
double prob_gout =
xs_t->chi_delayed[cache[tid].a][gin][gout][dg];
while (prob_gout < xi_gout) {
gout++;
prob_gout +=
xs_t->chi_delayed[cache[tid].a][gin][gout][dg];
}
}
}
//==============================================================================
void
Mgxs::sample_scatter(int gin, int& gout, double& mu, double& wgt)
{
// This method assumes that the temperature and angle indices are set
// Sample the data
#ifdef _OPENMP
int tid = omp_get_thread_num();
#else
int tid = 0;
#endif
xs[cache[tid].t].scatter[cache[tid].a]->sample(gin, gout, mu, wgt);
}
//==============================================================================
void
Mgxs::calculate_xs(int gin, double sqrtkT, const double uvw[3],
double& total_xs, double& abs_xs, double& nu_fiss_xs)
{
// Set our indices
#ifdef _OPENMP
int tid = omp_get_thread_num();
#else
int tid = 0;
#endif
set_temperature_index(sqrtkT);
set_angle_index(uvw);
XsData* xs_t = &xs[cache[tid].t];
total_xs = xs_t->total[cache[tid].a][gin];
abs_xs = xs_t->absorption[cache[tid].a][gin];
nu_fiss_xs = fissionable ? xs_t->nu_fission[cache[tid].a][gin] : 0.;
}
//==============================================================================
bool
Mgxs::equiv(const Mgxs& that)
{
return ((num_delayed_groups == that.num_delayed_groups) &&
(num_groups == that.num_groups) &&
(n_pol == that.n_pol) &&
(n_azi == that.n_azi) &&
(std::equal(polar.begin(), polar.end(), that.polar.begin())) &&
(std::equal(azimuthal.begin(), azimuthal.end(), that.azimuthal.begin())) &&
(scatter_format == that.scatter_format));
}
//==============================================================================
void
Mgxs::set_temperature_index(double sqrtkT)
{
// See if we need to find the new index
#ifdef _OPENMP
int tid = omp_get_thread_num();
#else
int tid = 0;
#endif
if (sqrtkT != cache[tid].sqrtkT) {
double kT = sqrtkT * sqrtkT;
// initialize vector for storage of the differences
std::valarray<double> temp_diff(kTs.data(), kTs.size());
// Find the minimum difference of kT and kTs
temp_diff = std::abs(temp_diff - kT);
cache[tid].t = std::min_element(std::begin(temp_diff), std::end(temp_diff)) -
std::begin(temp_diff);
// store this temperature as the last one used
cache[tid].sqrtkT = sqrtkT;
}
}
//==============================================================================
void
Mgxs::set_angle_index(const double uvw[3])
{
// See if we need to find the new index
#ifdef _OPENMP
int tid = omp_get_thread_num();
#else
int tid = 0;
#endif
if (!is_isotropic &&
((uvw[0] != cache[tid].u) || (uvw[1] != cache[tid].v) ||
(uvw[2] != cache[tid].w))) {
// convert uvw to polar and azimuthal angles
double my_pol = std::acos(uvw[2]);
double my_azi = std::atan2(uvw[1], uvw[0]);
// Find the location, assuming equal-bin angles
double delta_angle = PI / n_pol;
int p = std::floor(my_pol / delta_angle);
delta_angle = 2. * PI / n_azi;
int a = std::floor((my_azi + PI) / delta_angle);
cache[tid].a = n_azi * p + a;
// store this direction as the last one used
cache[tid].u = uvw[0];
cache[tid].v = uvw[1];
cache[tid].w = uvw[2];
}
}
} // namespace openmc

205
src/mgxs.h Normal file
View file

@ -0,0 +1,205 @@
//! \file mgxs.h
//! A collection of classes for Multi-Group Cross Section data
#ifndef MGXS_H
#define MGXS_H
#include <string>
#include <vector>
#include "constants.h"
#include "hdf5_interface.h"
#include "xsdata.h"
namespace openmc {
//==============================================================================
// Cache contains the cached data for an MGXS object
//==============================================================================
struct CacheData {
double sqrtkT; // last temperature corresponding to t
int t; // temperature index
int a; // angle index
// last angle that corresponds to a
double u;
double v;
double w;
};
//==============================================================================
// MGXS contains the mgxs data for a nuclide/material
//==============================================================================
class Mgxs {
private:
double_1dvec kTs; // temperature in eV (k * T)
int scatter_format; // flag for if this is legendre, histogram, or tabular
int num_delayed_groups; // number of delayed neutron groups
int num_groups; // number of energy groups
std::vector<XsData> xs; // Cross section data
// MGXS Incoming Flux Angular grid information
bool is_isotropic; // used to skip search for angle indices if isotropic
int n_pol;
int n_azi;
double_1dvec polar;
double_1dvec azimuthal;
//! \brief Initializes the Mgxs object metadata
//!
//! @param in_name Name of the object.
//! @param in_awr atomic-weight ratio.
//! @param in_kTs temperatures (in units of eV) that data is available.
//! @param in_fissionable Is this item fissionable or not.
//! @param in_scatter_format Denotes whether Legendre, Tabular, or
//! Histogram scattering is used.
//! @param in_num_groups Number of energy groups.
//! @param in_num_delayed_groups Number of delayed groups.
//! @param in_is_isotropic Is this an isotropic or angular with respect to
//! the incoming particle.
//! @param in_polar Polar angle grid.
//! @param in_azimuthal Azimuthal angle grid.
void
init(const std::string& in_name, double in_awr, const double_1dvec& in_kTs,
bool in_fissionable, int in_scatter_format, int in_num_groups,
int in_num_delayed_groups, bool in_is_isotropic,
const double_1dvec& in_polar, const double_1dvec& in_azimuthal);
//! \brief Initializes the Mgxs object metadata from the HDF5 file
//!
//! @param xs_id HDF5 group id for the cross section data.
//! @param in_num_groups Number of energy groups.
//! @param in_num_delayed_groups Number of delayed groups.
//! @param temperature Temperatures to read.
//! @param tolerance Tolerance of temperature selection method.
//! @param temps_to_read Resultant list of temperatures in the library
//! to read which correspond to the requested temperatures.
//! @param order_dim Resultant dimensionality of the scattering order.
//! @param method Method of choosing nearest temperatures.
void
metadata_from_hdf5(hid_t xs_id, int in_num_groups,
int in_num_delayed_groups, const double_1dvec& temperature,
double tolerance, int_1dvec& temps_to_read, int& order_dim,
int& method);
//! \brief Performs the actual act of combining the microscopic data for a
//! single temperature.
//!
//! @param micros Microscopic objects to combine.
//! @param scalars Scalars to multiply the microscopic data by.
//! @param micro_ts The temperature index of the microscopic objects that
//! corresponds to the temperature of interest.
//! @param this_t The temperature index of the macroscopic object.
void
combine(const std::vector<Mgxs*>& micros, const double_1dvec& scalars,
const int_1dvec& micro_ts, int this_t);
//! \brief Checks to see if this and that are able to be combined
//!
//! This comparison is used when building macroscopic cross sections
//! from microscopic cross sections.
//! @param that The other Mgxs to compare to this one.
//! @return True if they can be combined, False otherwise.
bool equiv(const Mgxs& that);
public:
std::string name; // name of dataset, e.g., UO2
double awr; // atomic weight ratio
bool fissionable; // Is this fissionable
std::vector<CacheData> cache; // index and data cache
Mgxs() = default;
//! \brief Constructor that loads the Mgxs object from the HDF5 file
//!
//! @param xs_id HDF5 group id for the cross section data.
//! @param energy_groups Number of energy groups.
//! @param delayed_groups Number of delayed groups.
//! @param temperature Temperatures to read.
//! @param tolerance Tolerance of temperature selection method.
//! @param max_order Maximum order requested by the user;
//! this is only used for Legendre scattering.
//! @param legendre_to_tabular Flag to denote if any Legendre provided
//! should be converted to a Tabular representation.
//! @param legendre_to_tabular_points If a conversion is requested, this
//! provides the number of points to use in the tabular representation.
//! @param method Method of choosing nearest temperatures.
Mgxs(hid_t xs_id, int energy_groups,
int delayed_groups, const double_1dvec& temperature, double tolerance,
int max_order, bool legendre_to_tabular,
int legendre_to_tabular_points, int& method);
//! \brief Constructor that initializes and populates all data to build a
//! macroscopic cross section from microscopic cross section.
//!
//! @param in_name Name of the object.
//! @param mat_kTs temperatures (in units of eV) that data is needed.
//! @param micros Microscopic objects to combine.
//! @param atom_densities Atom densities of those microscopic quantities.
//! @param tolerance Tolerance of temperature selection method.
//! @param method Method of choosing nearest temperatures.
Mgxs(const std::string& in_name, const double_1dvec& mat_kTs,
const std::vector<Mgxs*>& micros, const double_1dvec& atom_densities,
double tolerance, int& method);
//! \brief Provides a cross section value given certain parameters
//!
//! @param xstype Type of cross section requested, according to the
//! enumerated constants.
//! @param gin Incoming energy group.
//! @param gout Outgoing energy group; use nullptr if irrelevant, or if a
//! sum is requested.
//! @param mu Cosine of the change-in-angle, for scattering quantities;
//! use nullptr if irrelevant.
//! @param dg delayed group index; use nullptr if irrelevant.
//! @return Requested cross section value.
double
get_xs(int xstype, int gin, int* gout, double* mu, int* dg);
//! \brief Samples the fission neutron energy and if prompt or delayed.
//!
//! @param gin Incoming energy group.
//! @param dg Sampled delayed group index.
//! @param gout Sampled outgoing energy group.
void
sample_fission_energy(int gin, int& dg, int& gout);
//! \brief Samples the outgoing energy and angle from a scatter event.
//!
//! @param gin Incoming energy group.
//! @param gout Sampled outgoing energy group.
//! @param mu Sampled cosine of the change-in-angle.
//! @param wgt Weight of the particle to be adjusted.
void
sample_scatter(int gin, int& gout, double& mu, double& wgt);
//! \brief Calculates cross section quantities needed for tracking.
//!
//! @param gin Incoming energy group.
//! @param sqrtkT Temperature of the material.
//! @param uvw Incoming particle direction.
//! @param total_xs Resultant total cross section.
//! @param abs_xs Resultant absorption cross section.
//! @param nu_fiss_xs Resultant nu-fission cross section.
void
calculate_xs(int gin, double sqrtkT, const double uvw[3],
double& total_xs, double& abs_xs, double& nu_fiss_xs);
//! \brief Sets the temperature index in cache given a temperature
//!
//! @param sqrtkT Temperature of the material.
void
set_temperature_index(double sqrtkT);
//! \brief Sets the angle index in cache given a direction
//!
//! @param uvw Incoming particle direction.
void
set_angle_index(const double uvw[3]);
};
} // namespace openmc
#endif // MGXS_H

113
src/mgxs_data.F90
View file

@ -1,5 +1,7 @@
module mgxs_data
use, intrinsic :: ISO_C_BINDING
use constants
use algorithm, only: find
use dict_header, only: DictCharInt
@ -7,11 +9,11 @@ module mgxs_data
use geometry_header, only: get_temperatures, cells
use hdf5_interface
use material_header, only: Material, materials, n_materials
use mgxs_header
use mgxs_interface
use nuclide_header, only: n_nuclides
use set_header, only: SetChar
use settings
use stl_vector, only: VectorReal
use stl_vector, only: VectorReal, VectorChar
use string, only: to_lower
implicit none
@ -27,14 +29,11 @@ contains
integer :: j ! index over nuclides in material
integer :: i_nuclide ! index in nuclides array
character(20) :: name ! name of library to load
integer :: representation ! Data representation
character(MAX_LINE_LEN) :: temp_str
type(Material), pointer :: mat
type(SetChar) :: already_read
integer(HID_T) :: file_id
integer(HID_T) :: xsdata_group
logical :: file_exists
type(VectorReal), allocatable :: temps(:)
type(VectorReal), allocatable, target :: temps(:)
character(MAX_WORD_LEN) :: word
integer, allocatable :: array(:)
@ -69,9 +68,6 @@ contains
&version supported by OpenMC.")
end if
! allocate arrays for MGXS storage and cross section cache
allocate(nuclides_MG(n_nuclides))
! ==========================================================================
! READ ALL MGXS CROSS SECTION TABLES
@ -80,89 +76,29 @@ contains
mat => materials(i)
NUCLIDE_LOOP: do j = 1, mat % n_nuclides
name = mat % names(j)
name = trim(mat % names(j)) // C_NULL_CHAR
i_nuclide = mat % nuclide(j)
if (.not. already_read % contains(name)) then
i_nuclide = mat % nuclide(j)
call add_mgxs_c(file_id, name, num_energy_groups, num_delayed_groups, &
temps(i_nuclide) % size(), temps(i_nuclide) % data, &
temperature_tolerance, max_order, &
logical(legendre_to_tabular, C_BOOL), &
legendre_to_tabular_points, temperature_method)
call write_message("Loading " // trim(name) // " data...", 6)
! Check to make sure cross section set exists in the library
if (object_exists(file_id, trim(name))) then
xsdata_group = open_group(file_id, trim(name))
else
call fatal_error("Data for '" // trim(name) // "' does not exist in "&
&// trim(path_cross_sections))
end if
! First find out the data representation
if (attribute_exists(xsdata_group, "representation")) then
call read_attribute(temp_str, xsdata_group, "representation")
if (trim(temp_str) == 'isotropic') then
representation = MGXS_ISOTROPIC
else if (trim(temp_str) == 'angle') then
representation = MGXS_ANGLE
else
call fatal_error("Invalid Data Representation!")
end if
else
! Default to isotropic representation
representation = MGXS_ISOTROPIC
end if
! Now allocate accordingly
select case(representation)
case(MGXS_ISOTROPIC)
allocate(MgxsIso :: nuclides_MG(i_nuclide) % obj)
case(MGXS_ANGLE)
allocate(MgxsAngle :: nuclides_MG(i_nuclide) % obj)
end select
! Now read in the data specific to the type we just declared
call nuclides_MG(i_nuclide) % obj % from_hdf5(xsdata_group, &
num_energy_groups, num_delayed_groups, temps(i_nuclide), &
temperature_method, temperature_tolerance, max_order, &
legendre_to_tabular, legendre_to_tabular_points)
! Add name to dictionary
call already_read % add(name)
call close_group(xsdata_group)
end if
end do NUCLIDE_LOOP
mat % fissionable = query_fissionable_c(mat % n_nuclides, mat % nuclide)
end do MATERIAL_LOOP
call file_close(file_id)
! Avoid some valgrind leak errors
call already_read % clear()
! Loop around material
MATERIAL_LOOP3: do i = 1, n_materials
! Get material
mat => materials(i)
! Loop around nuclides in material
NUCLIDE_LOOP2: do j = 1, mat % n_nuclides
! Is this fissionable?
if (nuclides_MG(mat % nuclide(j)) % obj % fissionable) then
mat % fissionable = .true.
end if
if (mat % fissionable) then
exit NUCLIDE_LOOP2
end if
end do NUCLIDE_LOOP2
end do MATERIAL_LOOP3
call file_close(file_id)
end subroutine read_mgxs
!===============================================================================
@ -173,8 +109,7 @@ contains
integer :: i_mat ! index in materials array
type(Material), pointer :: mat ! current material
type(VectorReal), allocatable :: kTs(:)
allocate(macro_xs(n_materials))
character(MAX_WORD_LEN) :: name ! name of material
! Get temperatures to read for each material
call get_mat_kTs(kTs)
@ -188,19 +123,13 @@ contains
! Get the material
mat => materials(i_mat)
! Get the scattering type for the first nuclide
select type(nuc => nuclides_MG(mat % nuclide(1)) % obj)
type is (MgxsIso)
allocate(MgxsIso :: macro_xs(i_mat) % obj)
type is (MgxsAngle)
allocate(MgxsAngle :: macro_xs(i_mat) % obj)
end select
name = trim(mat % name) // C_NULL_CHAR
! Do not read materials which we do not actually use in the problem to
! reduce storage
if (allocated(kTs(i_mat) % data)) then
call macro_xs(i_mat) % obj % combine(kTs(i_mat), mat, nuclides_MG, &
num_energy_groups, num_delayed_groups, max_order, &
call create_macro_xs_c(name, mat % n_nuclides, mat % nuclide, &
kTs(i_mat) % size(), kTs(i_mat) % data, mat % atom_density, &
temperature_tolerance, temperature_method)
end if
end do

3593
src/mgxs_header.F90
View file

File diff suppressed because it is too large Load diff

164
src/mgxs_interface.F90 Normal file
View file

@ -0,0 +1,164 @@
module mgxs_interface
use, intrinsic :: ISO_C_BINDING
use hdf5_interface
implicit none
interface
subroutine add_mgxs_c(file_id, name, energy_groups, delayed_groups, &
n_temps, temps, tolerance, max_order, legendre_to_tabular, &
legendre_to_tabular_points, method) bind(C)
use ISO_C_BINDING
import HID_T
implicit none
integer(HID_T), value, intent(in) :: file_id
character(kind=C_CHAR),intent(in) :: name(*)
integer(C_INT), value, intent(in) :: energy_groups
integer(C_INT), value, intent(in) :: delayed_groups
integer(C_INT), value, intent(in) :: n_temps
real(C_DOUBLE), intent(in) :: temps(1:n_temps)
real(C_DOUBLE), value, intent(in) :: tolerance
integer(C_INT), value, intent(in) :: max_order
logical(C_BOOL),value, intent(in) :: legendre_to_tabular
integer(C_INT), value, intent(in) :: legendre_to_tabular_points
integer(C_INT), intent(inout) :: method
end subroutine add_mgxs_c
function query_fissionable_c(n_nuclides, i_nuclides) result(result) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: n_nuclides
integer(C_INT), intent(in) :: i_nuclides(1:n_nuclides)
logical(C_BOOL) :: result
end function query_fissionable_c
subroutine create_macro_xs_c(name, n_nuclides, i_nuclides, n_temps, temps, &
atom_densities, tolerance, method) bind(C)
use ISO_C_BINDING
implicit none
character(kind=C_CHAR),intent(in) :: name(*)
integer(C_INT), value, intent(in) :: n_nuclides
integer(C_INT), intent(in) :: i_nuclides(1:n_nuclides)
integer(C_INT), value, intent(in) :: n_temps
real(C_DOUBLE), intent(in) :: temps(1:n_temps)
real(C_DOUBLE), intent(in) :: atom_densities(1:n_nuclides)
real(C_DOUBLE), value, intent(in) :: tolerance
integer(C_INT), intent(inout) :: method
end subroutine create_macro_xs_c
subroutine calculate_xs_c(i_mat, gin, sqrtkT, uvw, total_xs, abs_xs, &
nu_fiss_xs) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: i_mat
integer(C_INT), value, intent(in) :: gin
real(C_DOUBLE), value, intent(in) :: sqrtkT
real(C_DOUBLE), intent(in) :: uvw(1:3)
real(C_DOUBLE), intent(inout) :: total_xs
real(C_DOUBLE), intent(inout) :: abs_xs
real(C_DOUBLE), intent(inout) :: nu_fiss_xs
end subroutine calculate_xs_c
subroutine sample_scatter_c(i_mat, gin, gout, mu, wgt, uvw) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: i_mat
integer(C_INT), value, intent(in) :: gin
integer(C_INT), intent(inout) :: gout
real(C_DOUBLE), intent(inout) :: mu
real(C_DOUBLE), intent(inout) :: wgt
real(C_DOUBLE), intent(inout) :: uvw(1:3)
end subroutine sample_scatter_c
subroutine sample_fission_energy_c(i_mat, gin, dg, gout) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: i_mat
integer(C_INT), value, intent(in) :: gin
integer(C_INT), intent(inout) :: dg
integer(C_INT), intent(inout) :: gout
end subroutine sample_fission_energy_c
subroutine get_name_c(index, name_len, name) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
integer(C_INT), value, intent(in) :: name_len
character(kind=C_CHAR), intent(inout) :: name(name_len)
end subroutine get_name_c
function get_awr_c(index) result(awr) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
real(C_DOUBLE) :: awr
end function get_awr_c
function get_nuclide_xs_c(index, xstype, gin, gout, mu, dg) result(val) &
bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
integer(C_INT), value, intent(in) :: xstype
integer(C_INT), value, intent(in) :: gin
integer(C_INT), optional, intent(in) :: gout
real(C_DOUBLE), optional, intent(in) :: mu
integer(C_INT), optional, intent(in) :: dg
real(C_DOUBLE) :: val
end function get_nuclide_xs_c
function get_macro_xs_c(index, xstype, gin, gout, mu, dg) result(val) &
bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
integer(C_INT), value, intent(in) :: xstype
integer(C_INT), value, intent(in) :: gin
integer(C_INT), optional, intent(in) :: gout
real(C_DOUBLE), optional, intent(in) :: mu
integer(C_INT), optional, intent(in) :: dg
real(C_DOUBLE) :: val
end function get_macro_xs_c
subroutine set_nuclide_angle_index_c(index, uvw) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
real(C_DOUBLE), intent(in) :: uvw(1:3)
end subroutine set_nuclide_angle_index_c
subroutine set_macro_angle_index_c(index, uvw) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
real(C_DOUBLE), intent(in) :: uvw(1:3)
end subroutine set_macro_angle_index_c
subroutine set_nuclide_temperature_index_c(index, sqrtkT) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
real(C_DOUBLE), value, intent(in) :: sqrtkT
end subroutine set_nuclide_temperature_index_c
end interface
! Number of energy groups
integer(C_INT) :: num_energy_groups
! Number of delayed groups
integer(C_INT) :: num_delayed_groups
! Energy group structure with decreasing energy
real(8), allocatable :: energy_bins(:)
! Midpoint of the energy group structure
real(8), allocatable :: energy_bin_avg(:)
! Energy group structure with increasing energy
real(8), allocatable :: rev_energy_bins(:)
end module mgxs_interface

229
src/mgxs_interface.cpp Normal file
View file

@ -0,0 +1,229 @@
#include <string>
#include "error.h"
#include "math_functions.h"
#include "mgxs_interface.h"
namespace openmc {
//==============================================================================
// Mgxs data loading interface methods
//==============================================================================
void
add_mgxs_c(hid_t file_id, const char* name, int energy_groups,
int delayed_groups, int n_temps, const double temps[], double tolerance,
int max_order, bool legendre_to_tabular, int legendre_to_tabular_points,
int& method)
{
// Convert temps to a vector for the from_hdf5 function
double_1dvec temperature(temps, temps + n_temps);
write_message("Loading " + std::string(name) + " data...", 6);
// Check to make sure cross section set exists in the library
hid_t xs_grp;
if (object_exists(file_id, name)) {
xs_grp = open_group(file_id, name);
} else {
fatal_error("Data for " + std::string(name) + " does not exist in "
+ "provided MGXS Library");
}
Mgxs mg(xs_grp, energy_groups, delayed_groups, temperature, tolerance,
max_order, legendre_to_tabular, legendre_to_tabular_points, method);
nuclides_MG.push_back(mg);
close_group(xs_grp);
}
//==============================================================================
bool
query_fissionable_c(int n_nuclides, const int i_nuclides[])
{
bool result = false;
for (int n = 0; n < n_nuclides; n++) {
if (nuclides_MG[i_nuclides[n] - 1].fissionable) result = true;
}
return result;
}
//==============================================================================
void
create_macro_xs_c(const char* mat_name, int n_nuclides, const int i_nuclides[],
int n_temps, const double temps[], const double atom_densities[],
double tolerance, int& method)
{
if (n_temps > 0) {
// // Convert temps to a vector
double_1dvec temperature(temps, temps + n_temps);
// Convert atom_densities to a vector
double_1dvec atom_densities_vec(atom_densities,
atom_densities + n_nuclides);
// Build array of pointers to nuclides_MG's Mgxs objects needed for this
// material
std::vector<Mgxs*> mgxs_ptr(n_nuclides);
for (int n = 0; n < n_nuclides; n++) {
mgxs_ptr[n] = &nuclides_MG[i_nuclides[n] - 1];
}
Mgxs macro(mat_name, temperature, mgxs_ptr, atom_densities_vec,
tolerance, method);
macro_xs.emplace_back(macro);
} else {
// Preserve the ordering of materials by including a blank entry
Mgxs macro;
macro_xs.emplace_back(macro);
}
}
//==============================================================================
// Mgxs tracking/transport/tallying interface methods
//==============================================================================
void
calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3],
double& total_xs, double& abs_xs, double& nu_fiss_xs)
{
macro_xs[i_mat - 1].calculate_xs(gin - 1, sqrtkT, uvw, total_xs, abs_xs,
nu_fiss_xs);
}
//==============================================================================
void
sample_scatter_c(int i_mat, int gin, int& gout, double& mu, double& wgt,
double uvw[3])
{
int gout_c = gout - 1;
macro_xs[i_mat - 1].sample_scatter(gin - 1, gout_c, mu, wgt);
// adjust return value for fortran indexing
gout = gout_c + 1;
// Rotate the angle
rotate_angle_c(uvw, mu, nullptr);
}
//==============================================================================
void
sample_fission_energy_c(int i_mat, int gin, int& dg, int& gout)
{
int dg_c = 0;
int gout_c = 0;
macro_xs[i_mat - 1].sample_fission_energy(gin - 1, dg_c, gout_c);
// adjust return values for fortran indexing
dg = dg_c + 1;
gout = gout_c + 1;
}
//==============================================================================
double
get_nuclide_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg)
{
int gout_c;
int* gout_c_p;
int dg_c;
int* dg_c_p;
if (gout != nullptr) {
gout_c = *gout - 1;
gout_c_p = &gout_c;
} else {
gout_c_p = gout;
}
if (dg != nullptr) {
dg_c = *dg - 1;
dg_c_p = &dg_c;
} else {
dg_c_p = dg;
}
return nuclides_MG[index - 1].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p);
}
//==============================================================================
double
get_macro_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg)
{
int gout_c;
int* gout_c_p;
int dg_c;
int* dg_c_p;
if (gout != nullptr) {
gout_c = *gout - 1;
gout_c_p = &gout_c;
} else {
gout_c_p = gout;
}
if (dg != nullptr) {
dg_c = *dg - 1;
dg_c_p = &dg_c;
} else {
dg_c_p = dg;
}
return macro_xs[index - 1].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p);
}
//==============================================================================
void
set_nuclide_angle_index_c(int index, const double uvw[3])
{
// Update the values
nuclides_MG[index - 1].set_angle_index(uvw);
}
//==============================================================================
void
set_macro_angle_index_c(int index, const double uvw[3])
{
// Update the values
macro_xs[index - 1].set_angle_index(uvw);
}
//==============================================================================
void
set_nuclide_temperature_index_c(int index, double sqrtkT)
{
// Update the values
nuclides_MG[index - 1].set_temperature_index(sqrtkT);
}
//==============================================================================
// General Mgxs methods
//==============================================================================
void
get_name_c(int index, int name_len, char* name)
{
// First blank out our input string
std::string str(name_len, ' ');
std::strcpy(name, str.c_str());
// Now get the data and copy to the C-string
str = nuclides_MG[index - 1].name;
std::strcpy(name, str.c_str());
// Finally, remove the null terminator
name[std::strlen(name)] = ' ';
}
//==============================================================================
double
get_awr_c(int index)
{
return nuclides_MG[index - 1].awr;
}
} // namespace openmc

79
src/mgxs_interface.h Normal file
View file

@ -0,0 +1,79 @@
//! \file mgxs_interface.h
//! A collection of C interfaces to the C++ Mgxs class
#ifndef MGXS_INTERFACE_H
#define MGXS_INTERFACE_H
#include "hdf5_interface.h"
#include "mgxs.h"
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
extern std::vector<Mgxs> nuclides_MG;
extern std::vector<Mgxs> macro_xs;
//==============================================================================
// Mgxs data loading interface methods
//==============================================================================
extern "C" void
add_mgxs_c(hid_t file_id, const char* name, int energy_groups,
int delayed_groups, int n_temps, const double temps[], double tolerance,
int max_order, bool legendre_to_tabular, int legendre_to_tabular_points,
int& method);
extern "C" bool
query_fissionable_c(int n_nuclides, const int i_nuclides[]);
extern "C" void
create_macro_xs_c(const char* mat_name, int n_nuclides, const int i_nuclides[],
int n_temps, const double temps[], const double atom_densities[],
double tolerance, int& method);
//==============================================================================
// Mgxs tracking/transport/tallying interface methods
//==============================================================================
extern "C" void
calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3],
double& total_xs, double& abs_xs, double& nu_fiss_xs);
extern "C" void
sample_scatter_c(int i_mat, int gin, int& gout, double& mu, double& wgt,
double uvw[3]);
extern "C" void
sample_fission_energy_c(int i_mat, int gin, int& dg, int& gout);
extern "C" double
get_nuclide_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg);
extern "C" double
get_macro_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg);
extern "C" void
set_nuclide_angle_index_c(int index, const double uvw[3]);
extern "C" void
set_macro_angle_index_c(int index, const double uvw[3]);
extern "C" void
set_nuclide_temperature_index_c(int index, double sqrtkT);
//==============================================================================
// General Mgxs methods
//==============================================================================
extern "C" void
get_name_c(int index, int name_len, char* name);
extern "C" double
get_awr_c(int index);
} // namespace openmc
#endif // MGXS_INTERFACE_H

8
src/nuclide_header.F90
View file

@ -165,10 +165,10 @@ module nuclide_header
!===============================================================================
type MaterialMacroXS
real(8) :: total ! macroscopic total xs
real(8) :: absorption ! macroscopic absorption xs
real(8) :: fission ! macroscopic fission xs
real(8) :: nu_fission ! macroscopic production xs
real(C_DOUBLE) :: total ! macroscopic total xs
real(C_DOUBLE) :: absorption ! macroscopic absorption xs
real(C_DOUBLE) :: fission ! macroscopic fission xs
real(C_DOUBLE) :: nu_fission ! macroscopic production xs
end type MaterialMacroXS
!===============================================================================

10
src/output.F90
View file

@ -12,7 +12,7 @@ module output
use math, only: t_percentile
use mesh_header, only: RegularMesh, meshes
use message_passing, only: master, n_procs
use mgxs_header, only: nuclides_MG
use mgxs_interface
use nuclide_header
use particle_header, only: LocalCoord, Particle
use plot_header
@ -76,7 +76,7 @@ contains
write(UNIT=OUTPUT_UNIT, FMT=*) &
' | The OpenMC Monte Carlo Code'
write(UNIT=OUTPUT_UNIT, FMT=*) &
' Copyright | 2011-2018 Massachusetts Institute of Technology'
' Copyright | 2011-2018 MIT and OpenMC contributors'
write(UNIT=OUTPUT_UNIT, FMT=*) &
' License | http://openmc.readthedocs.io/en/latest/license.html'
write(UNIT=OUTPUT_UNIT, FMT='(11X,"Version | ",I1,".",I2,".",I1)') &
@ -171,7 +171,7 @@ contains
write(UNIT=OUTPUT_UNIT, FMT='(1X,A,A)') "Git SHA1: ", GIT_SHA1
#endif
write(UNIT=OUTPUT_UNIT, FMT=*) "Copyright (c) 2011-2018 &
&Massachusetts Institute of Technology"
&Massachusetts Institute of Technology and OpenMC contributors"
write(UNIT=OUTPUT_UNIT, FMT=*) "MIT/X license at &
&<http://openmc.readthedocs.io/en/latest/license.html>"
end if
@ -680,6 +680,7 @@ contains
character(36) :: score_name ! names of scoring function
! to be applied at write-time
type(TallyFilterMatch), allocatable :: matches(:)
character(MAX_WORD_LEN) :: temp_name
! Skip if there are no tallies
if (n_tallies == 0) return
@ -843,8 +844,9 @@ contains
write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
trim(nuclides(i_nuclide) % name)
else
call get_name_c(i_nuclide, len(temp_name), temp_name)
write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
trim(nuclides_MG(i_nuclide) % obj % name)
trim(temp_name)
end if
end if

2
src/particle_restart.F90
View file

@ -6,7 +6,7 @@ module particle_restart
use constants
use error, only: write_message
use hdf5_interface, only: file_open, file_close, read_dataset, HID_T
use mgxs_header, only: energy_bin_avg
use mgxs_interface, only: energy_bin_avg
use nuclide_header, only: micro_xs, n_nuclides
use output, only: print_particle
use particle_header, only: Particle

17
src/physics_mg.F90
View file

@ -8,13 +8,12 @@ module physics_mg
use material_header, only: Material, materials
use math, only: rotate_angle
use mesh_header, only: meshes
use mgxs_header
use mgxs_interface
use message_passing
use nuclide_header, only: material_xs
use particle_header, only: Particle
use physics_common
use random_lcg, only: prn
use scattdata_header
use settings
use simulation_header
use string, only: to_str
@ -143,16 +142,12 @@ contains
type(Particle), intent(inout) :: p
call macro_xs(p % material) % obj % sample_scatter(p % coord(1) % uvw, &
p % last_g, p % g, &
p % mu, p % wgt)
call sample_scatter_c(p % material, p % last_g, p % g, p % mu, &
p % wgt, p % coord(1) % uvw)
! Update energy value for downstream compatability (in tallying)
p % E = energy_bin_avg(p % g)
! Convert change in angle (mu) to new direction
p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, p % mu)
! Set event component
p % event = EVENT_SCATTER
@ -178,10 +173,6 @@ contains
real(8) :: mu ! fission neutron angular cosine
real(8) :: phi ! fission neutron azimuthal angle
real(8) :: weight ! weight adjustment for ufs method
class(Mgxs), pointer :: xs
! Get Pointers
xs => macro_xs(p % material) % obj
! TODO: Heat generation from fission
@ -261,7 +252,7 @@ contains
! Sample secondary energy distribution for fission reaction and set energy
! in fission bank
call xs % sample_fission_energy(p % g, bank_array(i) % uvw, dg, gout)
call sample_fission_energy_c(p % material, p % g, dg, gout)
bank_array(i) % E = real(gout, 8)
bank_array(i) % delayed_group = dg

936
src/scattdata.cpp Normal file
View file

@ -0,0 +1,936 @@
#include <algorithm>
#include <numeric>
#include <cmath>
#include "constants.h"
#include "math_functions.h"
#include "random_lcg.h"
#include "error.h"
#include "scattdata.h"
namespace openmc {
//==============================================================================
// ScattData base-class methods
//==============================================================================
void
ScattData::base_init(int order, const int_1dvec& in_gmin,
const int_1dvec& in_gmax, const double_2dvec& in_energy,
const double_2dvec& in_mult)
{
int groups = in_energy.size();
gmin = in_gmin;
gmax = in_gmax;
energy.resize(groups);
mult.resize(groups);
dist.resize(groups);
for (int gin = 0; gin < groups; gin++) {
// Store the inputted data
energy[gin] = in_energy[gin];
mult[gin] = in_mult[gin];
// Make sure the energy is normalized
double norm = std::accumulate(energy[gin].begin(), energy[gin].end(), 0.);
if (norm != 0.) {
for (auto& n : energy[gin]) n /= norm;
}
// Initialize the distribution data
dist[gin].resize(in_gmax[gin] - in_gmin[gin] + 1);
for (auto& v : dist[gin]) {
v.resize(order);
}
}
}
//==============================================================================
void
ScattData::base_combine(int max_order,
const std::vector<ScattData*>& those_scatts, const double_1dvec& scalars,
int_1dvec& in_gmin, int_1dvec& in_gmax, double_2dvec& sparse_mult,
double_3dvec& sparse_scatter)
{
int groups = those_scatts[0] -> energy.size();
// Now allocate and zero our storage spaces
double_3dvec this_matrix = double_3dvec(groups, double_2dvec(groups,
double_1dvec(max_order, 0.)));
double_2dvec mult_numer(groups, double_1dvec(groups, 0.));
double_2dvec mult_denom(groups, double_1dvec(groups, 0.));
// Build the dense scattering and multiplicity matrices
// Get the multiplicity_matrix
// To combine from nuclidic data we need to use the final relationship
// mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) /
// sum_i(N_i*(nuscatt_{i,g,g'} / mult_{i,g,g'}))
// Developed as follows:
// mult_{gg'} = nuScatt{g,g'} / Scatt{g,g'}
// mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) / sum(N_i*scatt_{i,g,g'})
// mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) /
// sum_i(N_i*(nuscatt_{i,g,g'} / mult_{i,g,g'}))
// nuscatt_{i,g,g'} can be reconstructed from the energy and scattxs member
// variables
for (int i = 0; i < those_scatts.size(); i++) {
ScattData* that = those_scatts[i];
// Build the dense matrix for that object
double_3dvec that_matrix = that->get_matrix(max_order);
// Now add that to this for the scattering and multiplicity
for (int gin = 0; gin < groups; gin++) {
// Only spend time adding that's gmin to gmax data since the rest will
// be zeros
int i_gout = 0;
for (int gout = that->gmin[gin]; gout <= that->gmax[gin]; gout++) {
// Do the scattering matrix
for (int l = 0; l < max_order; l++) {
this_matrix[gin][gout][l] += scalars[i] * that_matrix[gin][gout][l];
}
// Incorporate that's contribution to the multiplicity matrix data
double nuscatt = that->scattxs[gin] * that->energy[gin][i_gout];
mult_numer[gin][gout] += scalars[i] * nuscatt;
if (that->mult[gin][i_gout] > 0.) {
mult_denom[gin][gout] += scalars[i] * nuscatt / that->mult[gin][i_gout];
} else {
mult_denom[gin][gout] += scalars[i];
}
i_gout++;
}
}
}
// Combine mult_numer and mult_denom into the combined multiplicity matrix
double_2dvec this_mult(groups, double_1dvec(groups, 1.));
for (int gin = 0; gin < groups; gin++) {
for (int gout = 0; gout < groups; gout++) {
if (mult_denom[gin][gout] > 0.) {
this_mult[gin][gout] = mult_numer[gin][gout] / mult_denom[gin][gout];
}
}
}
mult_numer.clear();
mult_denom.clear();
// We have the data, now we need to convert to a jagged array and then use
// the initialize function to store it on the object.
for (int gin = 0; gin < groups; gin++) {
// Find the minimum and maximum group boundaries
int gmin_;
for (gmin_ = 0; gmin_ < groups; gmin_++) {
bool non_zero = false;
for (int l = 0; l < this_matrix[gin][gmin_].size(); l++) {
if (this_matrix[gin][gmin_][l] != 0.) {
non_zero = true;
break;
}
}
if (non_zero) break;
}
int gmax_;
for (gmax_ = groups - 1; gmax_ >= 0; gmax_--) {
bool non_zero = false;
for (int l = 0; l < this_matrix[gin][gmax_].size(); l++) {
if (this_matrix[gin][gmax_][l] != 0.) {
non_zero = true;
break;
}
}
if (non_zero) break;
}
// treat the case of all values being 0
if (gmin_ > gmax_) {
gmin_ = gin;
gmax_ = gin;
}
// Store the group bounds
in_gmin[gin] = gmin_;
in_gmax[gin] = gmax_;
// Store the data in the compressed format
sparse_scatter[gin].resize(gmax_ - gmin_ + 1);
sparse_mult[gin].resize(gmax_ - gmin_ + 1);
int i_gout = 0;
for (int gout = gmin_; gout <= gmax_; gout++) {
sparse_scatter[gin][i_gout] = this_matrix[gin][gout];
sparse_mult[gin][i_gout] = this_mult[gin][gout];
i_gout++;
}
}
}
//==============================================================================
void
ScattData::sample_energy(int gin, int& gout, int& i_gout)
{
// Sample the outgoing group
double xi = prn();
i_gout = 0;
gout = gmin[gin];
double prob = energy[gin][i_gout];
while((prob < xi) && (gout < gmax[gin])) {
gout++;
i_gout++;
prob += energy[gin][i_gout];
}
}
//==============================================================================
double
ScattData::get_xs(int xstype, int gin, const int* gout, const double* mu)
{
// Set the outgoing group offset index as needed
int i_gout = 0;
if (gout != nullptr) {
// short circuit the function if gout is from a zero portion of the
// scattering matrix
if ((*gout < gmin[gin]) || (*gout > gmax[gin])) { // > gmax?
return 0.;
}
i_gout = *gout - gmin[gin];
}
double val = scattxs[gin];
switch(xstype) {
case MG_GET_XS_SCATTER:
if (gout != nullptr) val *= energy[gin][i_gout];
break;
case MG_GET_XS_SCATTER_MULT:
if (gout != nullptr) {
val *= energy[gin][i_gout] / mult[gin][i_gout];
} else {
val /= std::inner_product(mult[gin].begin(), mult[gin].end(),
energy[gin].begin(), 0.0);
}
break;
case MG_GET_XS_SCATTER_FMU_MULT:
if ((gout != nullptr) && (mu != nullptr)) {
val *= energy[gin][i_gout] * calc_f(gin, *gout, *mu);
} else {
// This is not an expected path (asking for f_mu without asking for a
// group or mu is not useful
fatal_error("Invalid call to get_xs");
}
break;
case MG_GET_XS_SCATTER_FMU:
if ((gout != nullptr) && (mu != nullptr)) {
val *= energy[gin][i_gout] * calc_f(gin, *gout, *mu) / mult[gin][i_gout];
} else {
// This is not an expected path (asking for f_mu without asking for a
// group or mu is not useful
fatal_error("Invalid call to get_xs");
}
break;
}
return val;
}
//==============================================================================
// ScattDataLegendre methods
//==============================================================================
void
ScattDataLegendre::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs)
{
int groups = coeffs.size();
int order = coeffs[0][0].size();
// make a copy of coeffs that we can use to both extract data and normalize
double_3dvec matrix = coeffs;
// Get the scattering cross section value by summing the un-normalized P0
// coefficient in the variable matrix over all outgoing groups.
scattxs.resize(groups);
for (int gin = 0; gin < groups; gin++) {
int num_groups = in_gmax[gin] - in_gmin[gin] + 1;
scattxs[gin] = 0.;
for (int i_gout = 0; i_gout < num_groups; i_gout++) {
scattxs[gin] += matrix[gin][i_gout][0];
}
}
// Build the energy transfer matrix from data in the variable matrix while
// also normalizing the variable matrix itself
// (forcing the CDF of f(mu=1) == 1)
double_2dvec in_energy;
in_energy.resize(groups);
for (int gin = 0; gin < groups; gin++) {
int num_groups = in_gmax[gin] - in_gmin[gin] + 1;
in_energy[gin].resize(num_groups);
for (int i_gout = 0; i_gout < num_groups; i_gout++) {
double norm = matrix[gin][i_gout][0];
in_energy[gin][i_gout] = norm;
if (norm != 0.) {
for (auto& n : matrix[gin][i_gout]) n /= norm;
}
}
}
// Initialize the base class attributes
ScattData::base_init(order, in_gmin, in_gmax, in_energy, in_mult);
// Set the distribution (sdata.dist) values and initialize max_val
max_val.resize(groups);
for (int gin = 0; gin < groups; gin++) {
int num_groups = gmax[gin] - gmin[gin] + 1;
for (int i_gout = 0; i_gout < num_groups; i_gout++) {
dist[gin][i_gout] = matrix[gin][i_gout];
}
max_val[gin].resize(num_groups);
for (auto& n : max_val[gin]) n = 0.;
}
// Now update the maximum value
update_max_val();
}
//==============================================================================
void
ScattDataLegendre::update_max_val()
{
int groups = max_val.size();
// Step through the polynomial with fixed number of points to identify the
// maximal value
int Nmu = 1001;
double dmu = 2. / (Nmu - 1);
for (int gin = 0; gin < groups; gin++) {
int num_groups = gmax[gin] - gmin[gin] + 1;
for (int i_gout = 0; i_gout < num_groups; i_gout++) {
for (int imu = 0; imu < Nmu; imu++) {
double mu;
if (imu == 0) {
mu = -1.;
} else if (imu == (Nmu - 1)) {
mu = 1.;
} else {
mu = -1. + (imu - 1) * dmu;
}
// Calculate probability
double f = evaluate_legendre_c(dist[gin][i_gout].size() - 1,
dist[gin][i_gout].data(), mu);
// if this is a new maximum, store it
if (f > max_val[gin][i_gout]) max_val[gin][i_gout] = f;
} // end imu loop
// Since we may not have caught the true max, add 10% margin
max_val[gin][i_gout] *= 1.1;
}
}
}
//==============================================================================
double
ScattDataLegendre::calc_f(int gin, int gout, double mu)
{
double f;
if ((gout < gmin[gin]) || (gout > gmax[gin])) {
f = 0.;
} else {
int i_gout = gout - gmin[gin];
f = evaluate_legendre_c(dist[gin][i_gout].size() - 1,
dist[gin][i_gout].data(), mu);
}
return f;
}
//==============================================================================
void
ScattDataLegendre::sample(int gin, int& gout, double& mu, double& wgt)
{
// Sample the outgoing energy using the base-class method
int i_gout;
sample_energy(gin, gout, i_gout);
// Now we can sample mu using the scattering kernel using rejection
// sampling from a rectangular bounding box
double M = max_val[gin][i_gout];
int samples = 0;
while(true) {
mu = 2. * prn() - 1.;
double f = calc_f(gin, gout, mu);
if (f > 0.) {
double u = prn() * M;
if (u <= f) break;
}
samples++;
if (samples > MAX_SAMPLE) {
fatal_error("Maximum number of Legendre expansion samples reached");
}
};
// Update the weight to reflect neutron multiplicity
wgt *= mult[gin][i_gout];
}
//==============================================================================
void
ScattDataLegendre::combine(const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars)
{
// Find the max order in the data set and make sure we can combine the sets
int max_order = 0;
for (int i = 0; i < those_scatts.size(); i++) {
// Lets also make sure these items are combineable
ScattDataLegendre* that = dynamic_cast<ScattDataLegendre*>(those_scatts[i]);
if (!that) {
fatal_error("Cannot combine the ScattData objects!");
}
int that_order = that->get_order();
if (that_order > max_order) max_order = that_order;
}
max_order++; // Add one since this is a Legendre
int groups = those_scatts[0] -> energy.size();
int_1dvec in_gmin(groups);
int_1dvec in_gmax(groups);
double_3dvec sparse_scatter(groups);
double_2dvec sparse_mult(groups);
// The rest of the steps do not depend on the type of angular representation
// so we use a base class method to sum up xs and create new energy and mult
// matrices
ScattData::base_combine(max_order, those_scatts, scalars, in_gmin, in_gmax,
sparse_mult, sparse_scatter);
// Got everything we need, store it.
init(in_gmin, in_gmax, sparse_mult, sparse_scatter);
}
//==============================================================================
double_3dvec
ScattDataLegendre::get_matrix(int max_order)
{
// Get the sizes and initialize the data to 0
int groups = energy.size();
int order_dim = max_order + 1;
double_3dvec matrix = double_3dvec(groups, double_2dvec(groups,
double_1dvec(order_dim, 0.)));
for (int gin = 0; gin < groups; gin++) {
for (int i_gout = 0; i_gout < energy[gin].size(); i_gout++) {
int gout = i_gout + gmin[gin];
for (int l = 0; l < order_dim; l++) {
matrix[gin][gout][l] = scattxs[gin] * energy[gin][i_gout] *
dist[gin][i_gout][l];
}
}
}
return matrix;
}
//==============================================================================
// ScattDataHistogram methods
//==============================================================================
void
ScattDataHistogram::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs)
{
int groups = coeffs.size();
int order = coeffs[0][0].size();
// make a copy of coeffs that we can use to both extract data and normalize
double_3dvec matrix = coeffs;
// Get the scattering cross section value by summing the distribution
// over all the histogram bins in angle and outgoing energy groups
scattxs.resize(groups);
for (int gin = 0; gin < groups; gin++) {
scattxs[gin] = 0.;
for (int i_gout = 0; i_gout < matrix[gin].size(); i_gout++) {
scattxs[gin] += std::accumulate(matrix[gin][i_gout].begin(),
matrix[gin][i_gout].end(), 0.);
}
}
// Build the energy transfer matrix from data in the variable matrix
double_2dvec in_energy;
in_energy.resize(groups);
for (int gin = 0; gin < groups; gin++) {
int num_groups = in_gmax[gin] - in_gmin[gin] + 1;
in_energy[gin].resize(num_groups);
for (int i_gout = 0; i_gout < num_groups; i_gout++) {
double norm = std::accumulate(matrix[gin][i_gout].begin(),
matrix[gin][i_gout].end(), 0.);
in_energy[gin][i_gout] = norm;
if (norm != 0.) {
for (auto& n : matrix[gin][i_gout]) n /= norm;
}
}
}
// Initialize the base class attributes
ScattData::base_init(order, in_gmin, in_gmax, in_energy, in_mult);
// Build the angular distribution mu values
mu = double_1dvec(order);
dmu = 2. / order;
mu[0] = -1.;
for (int imu = 1; imu < order; imu++) {
mu[imu] = -1. + imu * dmu;
}
// Calculate f(mu) and integrate it so we can avoid rejection sampling
fmu.resize(groups);
for (int gin = 0; gin < groups; gin++) {
int num_groups = gmax[gin] - gmin[gin] + 1;
fmu[gin].resize(num_groups);
for (int i_gout = 0; i_gout < num_groups; i_gout++) {
fmu[gin][i_gout].resize(order);
// The variable matrix contains f(mu); so directly assign it
fmu[gin][i_gout] = matrix[gin][i_gout];
// Integrate the histogram
dist[gin][i_gout][0] = dmu * matrix[gin][i_gout][0];
for (int imu = 1; imu < order; imu++) {
dist[gin][i_gout][imu] = dmu * matrix[gin][i_gout][imu] +
dist[gin][i_gout][imu - 1];
}
// Now re-normalize for integral to unity
double norm = dist[gin][i_gout][order - 1];
if (norm > 0.) {
for (int imu = 0; imu < order; imu++) {
fmu[gin][i_gout][imu] /= norm;
dist[gin][i_gout][imu] /= norm;
}
}
}
}
}
//==============================================================================
double
ScattDataHistogram::calc_f(int gin, int gout, double mu)
{
double f;
if ((gout < gmin[gin]) || (gout > gmax[gin])) {
f = 0.;
} else {
// Find mu bin
int i_gout = gout - gmin[gin];
int imu;
if (mu == 1.) {
// use size -2 to have the index one before the end
imu = this->mu.size() - 2;
} else {
imu = std::floor((mu + 1.) / dmu + 1.) - 1;
}
f = fmu[gin][i_gout][imu];
}
return f;
}
//==============================================================================
void
ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt)
{
// Sample the outgoing energy using the base-class method
int i_gout;
sample_energy(gin, gout, i_gout);
// Determine the outgoing cosine bin
double xi = prn();
int imu;
if (xi < dist[gin][i_gout][0]) {
imu = 0;
} else {
// TODO lower_bound? + 1?
imu = std::upper_bound(dist[gin][i_gout].begin(),
dist[gin][i_gout].end(), xi) -
dist[gin][i_gout].begin();
}
// Randomly select mu within the imu bin
mu = prn() * dmu + this->mu[imu];
if (mu < -1.) {
mu = -1.;
} else if (mu > 1.) {
mu = 1.;
}
// Update the weight to reflect neutron multiplicity
wgt *= mult[gin][i_gout];
}
//==============================================================================
double_3dvec
ScattDataHistogram::get_matrix(int max_order)
{
// Get the sizes and initialize the data to 0
int groups = energy.size();
// We ignore the requested order for Histogram and Tabular representations
int order_dim = get_order();
double_3dvec matrix = double_3dvec(groups, double_2dvec(groups,
double_1dvec(order_dim, 0.)));
for (int gin = 0; gin < groups; gin++) {
for (int i_gout = 0; i_gout < energy[gin].size(); i_gout++) {
int gout = i_gout + gmin[gin];
for (int l = 0; l < order_dim; l++) {
matrix[gin][gout][l] = scattxs[gin] * energy[gin][i_gout] *
fmu[gin][i_gout][l];
}
}
}
return matrix;
}
//==============================================================================
void
ScattDataHistogram::combine(const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars)
{
// Find the max order in the data set and make sure we can combine the sets
int max_order = those_scatts[0]->get_order();
for (int i = 0; i < those_scatts.size(); i++) {
// Lets also make sure these items are combineable
ScattDataHistogram* that = dynamic_cast<ScattDataHistogram*>(those_scatts[i]);
if (!that) {
fatal_error("Cannot combine the ScattData objects!");
}
if (max_order != that->get_order()) {
fatal_error("Cannot combine the ScattData objects!");
}
}
int groups = those_scatts[0] -> energy.size();
int_1dvec in_gmin(groups);
int_1dvec in_gmax(groups);
double_3dvec sparse_scatter(groups);
double_2dvec sparse_mult(groups);
// The rest of the steps do not depend on the type of angular representation
// so we use a base class method to sum up xs and create new energy and mult
// matrices
ScattData::base_combine(max_order, those_scatts, scalars, in_gmin, in_gmax,
sparse_mult, sparse_scatter);
// Got everything we need, store it.
init(in_gmin, in_gmax, sparse_mult, sparse_scatter);
}
//==============================================================================
// ScattDataTabular methods
//==============================================================================
void
ScattDataTabular::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs)
{
int groups = coeffs.size();
int order = coeffs[0][0].size();
// make a copy of coeffs that we can use to both extract data and normalize
double_3dvec matrix = coeffs;
// Build the angular distribution mu values
mu = double_1dvec(order);
dmu = 2. / (order - 1);
mu[0] = -1.;
for (int imu = 1; imu < order - 1; imu++) {
mu[imu] = -1. + imu * dmu;
}
mu[order - 1] = 1.;
// Get the scattering cross section value by integrating the distribution
// over all mu points and then combining over all outgoing groups
scattxs.resize(groups);
for (int gin = 0; gin < groups; gin++) {
scattxs[gin] = 0.;
for (int i_gout = 0; i_gout < matrix[gin].size(); i_gout++) {
for (int imu = 1; imu < order; imu++) {
scattxs[gin] += 0.5 * dmu * (matrix[gin][i_gout][imu - 1] +
matrix[gin][i_gout][imu]);
}
}
}
// Build the energy transfer matrix from data in the variable matrix
double_2dvec in_energy(groups);
for (int gin = 0; gin < groups; gin++) {
int num_groups = in_gmax[gin] - in_gmin[gin] + 1;
in_energy[gin].resize(num_groups);
for (int i_gout = 0; i_gout < num_groups; i_gout++) {
double norm = 0.;
for (int imu = 1; imu < order; imu++) {
norm += 0.5 * dmu * (matrix[gin][i_gout][imu - 1] +
matrix[gin][i_gout][imu]);
}
in_energy[gin][i_gout] = norm;
}
}
// Initialize the base class attributes
ScattData::base_init(order, in_gmin, in_gmax, in_energy, in_mult);
// Calculate f(mu) and integrate it so we can avoid rejection sampling
fmu.resize(groups);
for (int gin = 0; gin < groups; gin++) {
int num_groups = gmax[gin] - gmin[gin] + 1;
fmu[gin].resize(num_groups);
for (int i_gout = 0; i_gout < num_groups; i_gout++) {
fmu[gin][i_gout].resize(order);
// The variable matrix contains f(mu); so directly assign it
fmu[gin][i_gout] = matrix[gin][i_gout];
// Ensure positivity
for (auto& val : fmu[gin][i_gout]) {
if (val < 0.) val = 0.;
}
// Now re-normalize for numerical integration issues and to take care of
// the above negative fix-up. Also accrue the CDF
double norm = 0.;
for (int imu = 1; imu < order; imu++) {
norm += 0.5 * dmu * (fmu[gin][i_gout][imu - 1] +
fmu[gin][i_gout][imu]);
// incorporate to the CDF
dist[gin][i_gout][imu] = norm;
}
// now do the normalization
if (norm > 0.) {
for (int imu = 0; imu < order; imu++) {
fmu[gin][i_gout][imu] /= norm;
dist[gin][i_gout][imu] /= norm;
}
}
}
}
}
//==============================================================================
double
ScattDataTabular::calc_f(int gin, int gout, double mu)
{
double f;
if ((gout < gmin[gin]) || (gout > gmax[gin])) {
f = 0.;
} else {
// Find mu bin
int i_gout = gout - gmin[gin];
int imu;
if (mu == 1.) {
// use size -2 to have the index one before the end
imu = this->mu.size() - 2;
} else {
imu = std::floor((mu + 1.) / dmu + 1.) - 1;
}
double r = (mu - this->mu[imu]) / (this->mu[imu + 1] - this->mu[imu]);
f = (1. - r) * fmu[gin][i_gout][imu] + r * fmu[gin][i_gout][imu + 1];
}
return f;
}
//==============================================================================
void
ScattDataTabular::sample(int gin, int& gout, double& mu, double& wgt)
{
// Sample the outgoing energy using the base-class method
int i_gout;
sample_energy(gin, gout, i_gout);
// Determine the outgoing cosine bin
int NP = this->mu.size();
double xi = prn();
double c_k = dist[gin][i_gout][0];
int k;
for (k = 0; k < NP - 1; k++) {
double c_k1 = dist[gin][i_gout][k + 1];
if (xi < c_k1) break;
c_k = c_k1;
}
// Check to make sure k is <= NP - 1
k = std::min(k, NP - 2);
// Find the pdf values we want
double p0 = fmu[gin][i_gout][k];
double mu0 = this -> mu[k];
double p1 = fmu[gin][i_gout][k + 1];
double mu1 = this -> mu[k + 1];
if (p0 == p1) {
mu = mu0 + (xi - c_k) / p0;
} else {
double frac = (p1 - p0) / (mu1 - mu0);
mu = mu0 + (std::sqrt(std::max(0., p0 * p0 + 2. * frac * (xi - c_k)))
- p0) / frac;
}
if (mu < -1.) {
mu = -1.;
} else if (mu > 1.) {
mu = 1.;
}
// Update the weight to reflect neutron multiplicity
wgt *= mult[gin][i_gout];
}
//==============================================================================
double_3dvec
ScattDataTabular::get_matrix(int max_order)
{
// Get the sizes and initialize the data to 0
int groups = energy.size();
// We ignore the requested order for Histogram and Tabular representations
int order_dim = get_order();
double_3dvec matrix = double_3dvec(groups, double_2dvec(groups,
double_1dvec(order_dim, 0.)));
for (int gin = 0; gin < groups; gin++) {
for (int i_gout = 0; i_gout < energy[gin].size(); i_gout++) {
int gout = i_gout + gmin[gin];
for (int l = 0; l < order_dim; l++) {
matrix[gin][gout][l] = scattxs[gin] * energy[gin][i_gout] *
fmu[gin][i_gout][l];
}
}
}
return matrix;
}
//==============================================================================
void
ScattDataTabular::combine(const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars)
{
// Find the max order in the data set and make sure we can combine the sets
int max_order = those_scatts[0]->get_order();
for (int i = 0; i < those_scatts.size(); i++) {
// Lets also make sure these items are combineable
ScattDataTabular* that = dynamic_cast<ScattDataTabular*>(those_scatts[i]);
if (!that) {
fatal_error("Cannot combine the ScattData objects!");
}
if (max_order != that->get_order()) {
fatal_error("Cannot combine the ScattData objects!");
}
}
int groups = those_scatts[0] -> energy.size();
int_1dvec in_gmin(groups);
int_1dvec in_gmax(groups);
double_3dvec sparse_scatter(groups);
double_2dvec sparse_mult(groups);
// The rest of the steps do not depend on the type of angular representation
// so we use a base class method to sum up xs and create new energy and mult
// matrices
ScattData::base_combine(max_order, those_scatts, scalars, in_gmin, in_gmax,
sparse_mult, sparse_scatter);
// Got everything we need, store it.
init(in_gmin, in_gmax, sparse_mult, sparse_scatter);
}
//==============================================================================
// module-level methods
//==============================================================================
void
convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab,
int n_mu)
{
// See if the user wants us to figure out how many points to use
if (n_mu == C_NONE) {
// then we will use 2 pts if its P0, or the default if a higher order
// TODO use an error minimization algorithm that also picks n_mu
if (leg.get_order() == 0) {
n_mu = 2;
} else {
n_mu = DEFAULT_NMU;
}
}
tab.base_init(n_mu, leg.gmin, leg.gmax, leg.energy, leg.mult);
tab.scattxs = leg.scattxs;
// Build mu and dmu
tab.mu = double_1dvec(n_mu);
tab.dmu = 2. / (n_mu - 1);
tab.mu[0] = -1.;
for (int imu = 1; imu < n_mu - 1; imu++) {
tab.mu[imu] = -1. + imu * tab.dmu;
}
tab.mu[n_mu - 1] = 1.;
// Calculate f(mu) and integrate it so we can avoid rejection sampling
int groups = tab.energy.size();
tab.fmu.resize(groups);
for (int gin = 0; gin < groups; gin++) {
int num_groups = tab.gmax[gin] - tab.gmin[gin] + 1;
tab.fmu[gin].resize(num_groups);
for (int i_gout = 0; i_gout < num_groups; i_gout++) {
tab.fmu[gin][i_gout].resize(n_mu);
for (int imu = 0; imu < n_mu; imu++) {
tab.fmu[gin][i_gout][imu] =
evaluate_legendre_c(leg.dist[gin][i_gout].size() - 1,
leg.dist[gin][i_gout].data(), tab.mu[imu]);
}
// Ensure positivity
for (auto& val : tab.fmu[gin][i_gout]) {
if (val < 0.) val = 0.;
}
// Now re-normalize for numerical integration issues and to take care of
// the above negative fix-up. Also accrue the CDF
double norm = 0.;
tab.dist[gin][i_gout][0] = 0.;
for (int imu = 1; imu < n_mu; imu++) {
norm += 0.5 * tab.dmu * (tab.fmu[gin][i_gout][imu - 1] +
tab.fmu[gin][i_gout][imu]);
// incorporate to the CDF
tab.dist[gin][i_gout][imu] = norm;
}
// now do the normalization
if (norm > 0.) {
for (int imu = 0; imu < n_mu; imu++) {
tab.fmu[gin][i_gout][imu] /= norm;
tab.dist[gin][i_gout][imu] /= norm;
}
}
}
}
}
} // namespace openmc

258
src/scattdata.h Normal file
View file

@ -0,0 +1,258 @@
//! \file scattdata.h
//! A collection of multi-group scattering data classes
#ifndef SCATTDATA_H
#define SCATTDATA_H
#include <vector>
namespace openmc {
// forward declarations so we can name our friend functions
class ScattDataLegendre;
class ScattDataTabular;
//==============================================================================
// SCATTDATA contains all the data needed to describe the scattering energy and
// angular distribution data
//==============================================================================
class ScattData {
protected:
//! \brief Initializes the attributes of the base class.
void
base_init(int order, const int_1dvec& in_gmin, const int_1dvec& in_gmax,
const double_2dvec& in_energy, const double_2dvec& in_mult);
//! \brief Combines microscopic ScattDatas into a macroscopic one.
void
base_combine(int max_order, const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars, int_1dvec& in_gmin, int_1dvec& in_gmax,
double_2dvec& sparse_mult, double_3dvec& sparse_scatter);
public:
double_2dvec energy; // Normalized p0 matrix for sampling Eout
double_2dvec mult; // nu-scatter multiplication (nu-scatt/scatt)
double_3dvec dist; // Angular distribution
int_1dvec gmin; // minimum outgoing group
int_1dvec gmax; // maximum outgoing group
double_1dvec scattxs; // Isotropic Sigma_{s,g_{in}}
//! \brief Calculates the value of normalized f(mu).
//!
//! The value of f(mu) is normalized as in the integral of f(mu)dmu across
//! [-1,1] is 1.
//!
//! @param gin Incoming energy group of interest.
//! @param gout Outgoing energy group of interest.
//! @param mu Cosine of the change-in-angle of interest.
//! @return The value of f(mu).
virtual double
calc_f(int gin, int gout, double mu) = 0;
//! \brief Samples the outgoing energy and angle from the ScattData info.
//!
//! @param gin Incoming energy group.
//! @param gout Sampled outgoing energy group.
//! @param mu Sampled cosine of the change-in-angle.
//! @param wgt Weight of the particle to be adjusted.
virtual void
sample(int gin, int& gout, double& mu, double& wgt) = 0;
//! \brief Initializes the ScattData object from a given scatter and
//! multiplicity matrix.
//!
//! @param in_gmin List of minimum outgoing groups for every incoming group
//! @param in_gmax List of maximum outgoing groups for every incoming group
//! @param in_mult Input sparse multiplicity matrix
//! @param coeffs Input sparse scattering matrix
virtual void
init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs) = 0;
//! \brief Combines the microscopic data.
//!
//! @param those_scatts Microscopic objects to combine.
//! @param scalars Scalars to multiply the microscopic data by.
virtual void
combine(const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars) = 0;
//! \brief Getter for the dimensionality of the scattering order.
//!
//! If Legendre this is the "n" in "Pn"; for Tabular, this is the number
//! of points, and for Histogram this is the number of bins.
//!
//! @return The order.
virtual int
get_order() = 0;
//! \brief Builds a dense scattering matrix from the constituent parts
//!
//! @param max_order If Legendre this is the maximum value of "n" in "Pn"
//! requested; ignored otherwise.
//! @return The dense scattering matrix.
virtual double_3dvec
get_matrix(int max_order) = 0;
//! \brief Samples the outgoing energy from the ScattData info.
//!
//! @param gin Incoming energy group.
//! @param gout Sampled outgoing energy group.
//! @param i_gout Sampled outgoing energy group index.
void
sample_energy(int gin, int& gout, int& i_gout);
//! \brief Provides a cross section value given certain parameters
//!
//! @param xstype Type of cross section requested, according to the
//! enumerated constants.
//! @param gin Incoming energy group.
//! @param gout Outgoing energy group; use nullptr if irrelevant, or if a
//! sum is requested.
//! @param mu Cosine of the change-in-angle, for scattering quantities;
//! use nullptr if irrelevant.
//! @return Requested cross section value.
double
get_xs(int xstype, int gin, const int* gout, const double* mu);
};
//==============================================================================
// ScattDataLegendre represents the angular distributions as Legendre kernels
//==============================================================================
class ScattDataLegendre: public ScattData {
protected:
// Maximal value for rejection sampling from a rectangle
double_2dvec max_val;
// Friend convert_legendre_to_tabular so it has access to protected
// parameters
friend void
convert_legendre_to_tabular(ScattDataLegendre& leg,
ScattDataTabular& tab, int n_mu);
public:
void
init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs);
void
combine(const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars);
//! \brief Find the maximal value of the angular distribution to use as a
// bounding box with rejection sampling.
void
update_max_val();
double
calc_f(int gin, int gout, double mu);
void
sample(int gin, int& gout, double& mu, double& wgt);
int
get_order() {return dist[0][0].size() - 1;};
double_3dvec
get_matrix(int max_order);
};
//==============================================================================
// ScattDataHistogram represents the angular distributions as a histogram, as it
// would be if it came from a "mu" tally in OpenMC
//==============================================================================
class ScattDataHistogram: public ScattData {
protected:
double_1dvec mu; // Angle distribution mu bin boundaries
double dmu; // Quick storage of the spacing between the mu bin points
double_3dvec fmu; // The angular distribution histogram
public:
void
init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs);
void
combine(const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars);
double
calc_f(int gin, int gout, double mu);
void
sample(int gin, int& gout, double& mu, double& wgt);
int
get_order() {return dist[0][0].size();};
double_3dvec
get_matrix(int max_order);
};
//==============================================================================
// ScattDataTabular represents the angular distributions as a table of mu and
// f(mu)
//==============================================================================
class ScattDataTabular: public ScattData {
protected:
double_1dvec mu; // Angle distribution mu grid points
double dmu; // Quick storage of the spacing between the mu points
double_3dvec fmu; // The angular distribution function
// Friend convert_legendre_to_tabular so it has access to protected
// parameters
friend void
convert_legendre_to_tabular(ScattDataLegendre& leg,
ScattDataTabular& tab, int n_mu);
public:
void
init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs);
void
combine(const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars);
double
calc_f(int gin, int gout, double mu);
void
sample(int gin, int& gout, double& mu, double& wgt);
int
get_order() {return dist[0][0].size();};
double_3dvec get_matrix(int max_order);
};
//==============================================================================
// Function to convert Legendre functions to tabular
//==============================================================================
//! \brief Converts a ScattDatalegendre to a ScattDataHistogram
//!
//! @param leg The initial ScattDataLegendre object.
//! @param leg The resultant ScattDataTabular object.
//! @param n_mu The number of mu points to use when building the
//! ScattDataTabular object.
void
convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab,
int n_mu);
} // namespace openmc
#endif // SCATTDATA_H

852
src/scattdata_header.F90
View file

@ -1,852 +0,0 @@
module scattdata_header
use algorithm, only: binary_search
use constants
use error, only: fatal_error
use math
use random_lcg, only: prn
implicit none
!===============================================================================
! JAGGED1D and JAGGED2D is a type which allows for jagged 1-D or 2-D array.
!===============================================================================
type :: Jagged2D
real(8), allocatable :: data(:, :)
end type Jagged2D
type :: Jagged1D
real(8), allocatable :: data(:)
end type Jagged1D
!===============================================================================
! SCATTDATA contains all the data to describe the scattering energy and
! angular distribution
!===============================================================================
type, abstract :: ScattData
! The data attribute of the energy, mult, and dist arrays
! are not necessarily 1-indexed as they instead will be allocated
! from a minimum outgoing group to an outgoing minimum group.
! Normalized p0 matrix on its own for sampling energy
type(Jagged1D), allocatable :: energy(:) ! (Gin % data(Gout))
! Nu-scatter multiplication (i.e. nu-scatt/scatt)
type(Jagged1D), allocatable :: mult(:) ! (Gin % data(Gout))
! Angular distribution
type(Jagged2D), allocatable :: dist(:) ! (Gin % data(Order/Nmu, Gout)
integer, allocatable :: gmin(:) ! Minimum outgoing group
integer, allocatable :: gmax(:) ! Maximum outgoing group
real(8), allocatable :: scattxs(:) ! Isotropic Sigma_{s,g_{in}}
contains
procedure(scattdata_init_), deferred :: init ! Initializes ScattData
procedure(scattdata_calc_f_), deferred :: calc_f ! Calculates f, given mu
procedure(scattdata_sample_), deferred :: sample ! sample the scatter event
procedure :: get_matrix => scattdata_get_matrix ! Rebuild scattering matrix
end type ScattData
abstract interface
subroutine scattdata_init_(this, gmin, gmax, mult, coeffs)
import ScattData, Jagged1D, Jagged2D
class(ScattData), intent(inout) :: this ! Object to work with
integer, intent(in) :: gmin(:) ! Min Gout
integer, intent(in) :: gmax(:) ! Max Gout
type(Jagged1D), intent(in) :: mult(:) ! Scatter Prod'n Matrix
type(Jagged2D), intent(in) :: coeffs(:) ! Coefficients to use
end subroutine scattdata_init_
pure function scattdata_calc_f_(this, gin, gout, mu) result(f)
import ScattData
class(ScattData), intent(in) :: this ! Scattering Object to work with
integer, intent(in) :: gin ! Incoming Energy Group
integer, intent(in) :: gout ! Outgoing Energy Group
real(8), intent(in) :: mu ! Angle of interest
real(8) :: f ! Return value of f(mu)
end function scattdata_calc_f_
subroutine scattdata_sample_(this, gin, gout, mu, wgt)
import ScattData
class(ScattData), intent(in) :: this ! Scattering Object to work with
integer, intent(in) :: gin ! Incoming neutron group
integer, intent(out) :: gout ! Sampled outgoin group
real(8), intent(out) :: mu ! Sampled change in angle
real(8), intent(inout) :: wgt ! Particle weight
end subroutine scattdata_sample_
end interface
type, extends(ScattData) :: ScattDataLegendre
! Maximal value for rejection sampling from rectangle
type(Jagged1D), allocatable :: max_val(:) ! (Gin % data(Gout))
contains
procedure :: init => scattdatalegendre_init
procedure :: calc_f => scattdatalegendre_calc_f
procedure :: sample => scattdatalegendre_sample
end type ScattDataLegendre
type, extends(ScattData) :: ScattDataHistogram
real(8), allocatable :: mu(:) ! Mu bins
real(8) :: dmu ! Mu spacing
! Histogram of f(mu) (dist has CDF)
type(Jagged2D), allocatable :: fmu(:) ! (Gin % data(Order/Nmu x Gout)
contains
procedure :: init => scattdatahistogram_init
procedure :: calc_f => scattdatahistogram_calc_f
procedure :: sample => scattdatahistogram_sample
procedure :: get_matrix => scattdatahistogram_get_matrix
end type ScattDataHistogram
type, extends(ScattData) :: ScattDataTabular
real(8), allocatable :: mu(:) ! Mu bins
real(8) :: dmu ! Mu spacing
! PDF of f(mu) (dist has CDF)
type(Jagged2D), allocatable :: fmu(:) ! (Gin % data(Order/Nmu x Gout)
contains
procedure :: init => scattdatatabular_init
procedure :: calc_f => scattdatatabular_calc_f
procedure :: sample => scattdatatabular_sample
procedure :: get_matrix => scattdatatabular_get_matrix
end type ScattDataTabular
!===============================================================================
! SCATTDATACONTAINER allocatable array for storing ScattData Objects (for angle)
!===============================================================================
type ScattDataContainer
class(ScattData), allocatable :: obj
end type ScattDataContainer
contains
!===============================================================================
! SCATTDATA*_INIT builds the scattdata object
!===============================================================================
subroutine scattdata_init(this, order, gmin, gmax, energy, mult)
class(ScattData), intent(inout) :: this ! Object to work on
integer, intent(in) :: order ! Data Order
integer, intent(in) :: gmin(:) ! Min Gout
integer, intent(in) :: gmax(:) ! Max Gout
type(Jagged1D), intent(inout) :: energy(:) ! Energy Transfer Matrix
type(Jagged1D), intent(in) :: mult(:) ! Scatter Prod'n Matrix
integer :: groups, gin
real(8) :: norm
groups = size(energy, dim=1)
allocate(this % gmin(groups))
allocate(this % gmax(groups))
allocate(this % energy(groups))
allocate(this % mult(groups))
allocate(this % dist(groups))
this % gmin = gmin
this % gmax = gmax
! Set the outgoing energy PDF values
do gin = 1, groups
! Make sure energy is normalized (i.e., CDF is 1)
norm = sum(energy(gin) % data(:))
if (norm /= ZERO) energy(gin) % data(:) = energy(gin) % data(:) / norm
! Set the values
allocate(this % energy(gin) % data(gmin(gin):gmax(gin)))
this % energy(gin) % data(:) = energy(gin) % data(:)
allocate(this % mult(gin) % data(gmin(gin):gmax(gin)))
this % mult(gin) % data(gmin(gin):gmax(gin)) = &
mult(gin) % data(gmin(gin):gmax(gin))
allocate(this % dist(gin) % data(order, gmin(gin):gmax(gin)))
this % dist(gin) % data = ZERO
end do
end subroutine scattdata_init
subroutine scattdatalegendre_init(this, gmin, gmax, mult, coeffs)
class(ScattDataLegendre), intent(inout) :: this ! Object to work on
integer, intent(in) :: gmin(:) ! Min Gout
integer, intent(in) :: gmax(:) ! Max Gout
type(Jagged1D), intent(in) :: mult(:) ! Scatter Prod'n Matrix
type(Jagged2D), intent(in) :: coeffs(:) ! Coefficients to use
real(8) :: dmu, mu, f, norm
integer :: imu, Nmu, gout, gin, groups, order
type(Jagged1D), allocatable :: energy(:)
type(Jagged2D), allocatable :: matrix(:)
groups = size(coeffs)
order = size(coeffs(1) % data, dim=1)
! make a copy of coeffs that we can use to extract data and normalize
allocate(matrix(groups))
do gin = 1, groups
allocate(matrix(gin) % data(order, gmin(gin):gmax(gin)))
matrix(gin) % data = coeffs(gin) % data
end do
! Get scattxs value
allocate(this % scattxs(groups))
! Get this by summing the un-normalized P0 coefficient in matrix
! over all outgoing groups
do gin = 1, groups
this % scattxs(gin) = sum(matrix(gin) % data(1, :), dim=1)
end do
allocate(energy(groups))
! Build energy transfer probability matrix from data in matrix
! while also normalizing matrix itself (making CDF of f(mu=1)=1)
do gin = 1, groups
allocate(energy(gin) % data(gmin(gin):gmax(gin)))
energy(gin) % data = ZERO
do gout = gmin(gin), gmax(gin)
norm = matrix(gin) % data(1, gout)
energy(gin) % data(gout) = norm
if (norm /= ZERO) then
matrix(gin) % data(:, gout) = matrix(gin) % data(:, gout) / norm
end if
end do
end do
call scattdata_init(this, order, gmin, gmax, energy, mult)
allocate(this % max_val(groups))
! Set dist values from matrix and initialize max_val
do gin = 1, groups
do gout = gmin(gin), gmax(gin)
this % dist(gin) % data(:, gout) = matrix(gin) % data(:, gout)
end do
allocate(this % max_val(gin) % data(gmin(gin):gmax(gin)))
this % max_val(gin) % data(:) = ZERO
end do
! Step through the polynomial with fixed number of points to identify
! the maximal value.
Nmu = 1001
dmu = TWO / real(Nmu - 1, 8)
do gin = 1, groups
do gout = gmin(gin), gmax(gin)
do imu = 1, Nmu
! Update mu. Do first and last seperate to avoid float errors
if (imu == 1) then
mu = -ONE
else if (imu == Nmu) then
mu = ONE
else
mu = -ONE + real(imu - 1, 8) * dmu
end if
! Calculate probability
f = this % calc_f(gin,gout,mu)
! If this is a new max, store it.
if (f > this % max_val(gin) % data(gout)) &
this % max_val(gin) % data(gout) = f
end do
! Finally, since we may not have caught the exact max, add 10% margin
this % max_val(gin) % data(gout) = &
this % max_val(gin) % data(gout) * 1.1_8
end do
end do
end subroutine scattdatalegendre_init
subroutine scattdatahistogram_init(this, gmin, gmax, mult, coeffs)
class(ScattDataHistogram), intent(inout) :: this ! Object to work on
integer, intent(in) :: gmin(:) ! Min Gout
integer, intent(in) :: gmax(:) ! Max Gout
type(Jagged1D), intent(in) :: mult(:) ! Scatter Prod'n Matrix
type(Jagged2D), intent(in) :: coeffs(:) ! Coefficients to use
integer :: imu, gin, gout, groups, order
real(8) :: norm
type(Jagged1D), allocatable :: energy(:)
type(Jagged2D), allocatable :: matrix(:)
groups = size(coeffs)
order = size(coeffs(1) % data, dim=1)
! make a copy of coeffs that we can use to extract data and normalize
allocate(matrix(groups))
do gin = 1, groups
allocate(matrix(gin) % data(order, gmin(gin):gmax(gin)))
matrix(gin) % data(:, :) = coeffs(gin) % data(:, :)
end do
! Get scattxs value
allocate(this % scattxs(groups))
! Get this by summing the un-normalized angular distribution in matrix
! over all outgoing groups
do gin = 1, groups
this % scattxs(gin) = sum(matrix(gin) % data(:, :))
end do
allocate(energy(groups))
! Build energy transfer probability matrix from data in matrix
! while also normalizing matrix itself (making CDF of f(mu=1)=1)
do gin = 1, groups
allocate(energy(gin) % data(gmin(gin):gmax(gin)))
do gout = gmin(gin), gmax(gin)
norm = sum(matrix(gin) % data(:, gout))
energy(gin) % data(gout) = norm
if (norm /= ZERO) then
matrix(gin) % data(:, gout) = matrix(gin) % data(:, gout) / norm
end if
end do
end do
call scattdata_init(this, order, gmin, gmax, energy, mult)
allocate(this % mu(order))
this % dmu = TWO / real(order, 8)
this % mu(1) = -ONE
do imu = 2, order
this % mu(imu) = -ONE + real(imu - 1, 8) * this % dmu
end do
! Integrate this histogram so we can avoid rejection sampling while
! also saving the original histogram in fmu
allocate(this % fmu(groups))
do gin = 1, groups
allocate(this % fmu(gin) % data(order, gmin(gin):gmax(gin)))
do gout = gmin(gin), gmax(gin)
! Store the histogram
this % fmu(gin) % data(:, gout) = matrix(gin) % data(:, gout)
! Integrate the histogram
this % dist(gin) % data(1, gout) = &
this % dmu * matrix(gin) % data(1, gout)
do imu = 2, order
this % dist(gin) % data(imu, gout) = &
this % dmu * matrix(gin) % data(imu, gout) + &
this % dist(gin) % data(imu - 1, gout)
end do
! Normalize the integral to unity
norm = this % dist(gin) % data(order, gout)
if (norm > ZERO) then
this % fmu(gin) % data(:, gout) = &
this % fmu(gin) % data(:, gout) / norm
this % dist(gin) % data(:, gout) = &
this % dist(gin) % data(:, gout) / norm
end if
end do
end do
end subroutine scattdatahistogram_init
subroutine scattdatatabular_init(this, gmin, gmax, mult, coeffs)
class(ScattDataTabular), intent(inout) :: this ! Object to work on
integer, intent(in) :: gmin(:) ! Min Gout
integer, intent(in) :: gmax(:) ! Max Gout
type(Jagged1D), intent(in) :: mult(:) ! Scatter Prod'n Matrix
type(Jagged2D), intent(in) :: coeffs(:) ! Coefficients to use
integer :: imu, gin, gout, groups, order
real(8) :: norm
type(Jagged1D), allocatable :: energy(:)
type(Jagged2D), allocatable :: matrix(:)
groups = size(coeffs)
order = size(coeffs(1) % data, dim=1)
! make a copy of coeffs that we can use to extract data and normalize
allocate(matrix(groups))
do gin = 1, groups
allocate(matrix(gin) % data(order, gmin(gin):gmax(gin)))
matrix(gin) % data = coeffs(gin) % data
end do
! Build the angular distribution mu values
allocate(this % mu(order))
this % dmu = TWO / real(order - 1, 8)
this % mu(1) = -ONE
do imu = 2, order - 1
this % mu(imu) = -ONE + real(imu - 1, 8) * this % dmu
end do
this % mu(order) = ONE
! Get scattxs
allocate(this % scattxs(groups))
! Get this by integrating the scattering distribution over all mu points
! and then combining over all outgoing groups
! over all outgoing groups
do gin = 1, groups
norm = ZERO
do gout = gmin(gin), gmax(gin)
do imu = 2, order
norm = norm + HALF * this % dmu * &
(matrix(gin) % data(imu - 1, gout) + &
matrix(gin) % data(imu, gout))
end do
end do
this % scattxs(gin) = norm
end do
allocate(energy(groups))
! Build energy transfer probability matrix from data in matrix
do gin = 1, groups
allocate(energy(gin) % data(gmin(gin):gmax(gin)))
do gout = gmin(gin), gmax(gin)
norm = ZERO
do imu = 2, order
norm = norm + HALF * this % dmu * &
(matrix(gin) % data(imu - 1, gout) + &
matrix(gin) % data(imu, gout))
end do
energy(gin) % data(gout) = norm
end do
end do
call scattdata_init(this, order, gmin, gmax, energy, mult)
! Calculate f(mu) and integrate it so we can avoid rejection sampling
allocate(this % fmu(groups))
do gin = 1, groups
allocate(this % fmu(gin) % data(order, gmin(gin):gmax(gin)))
do gout = gmin(gin), gmax(gin)
! Coeffs contain f(mu), put in f(mu) as that is where the
! PDF lives
this % fmu(gin) % data(:, gout) = matrix(gin) % data(:, gout)
! Force positivity
do imu = 1, order
if (this % fmu(gin) % data(imu, gout) < ZERO) then
this % fmu(gin) % data(imu, gout) = ZERO
end if
end do
! Re-normalize fmu for numerical integration issues and in case
! the negative fix-up introduced un-normalized data while
! accruing the CDF
norm = ZERO
do imu = 2, order
norm = norm + HALF * this % dmu * &
(this % fmu(gin) % data(imu - 1, gout) + &
this % fmu(gin) % data(imu, gout))
this % dist(gin) % data(imu, gout) = norm
end do
if (norm > ZERO) then
this % fmu(gin) % data(:, gout) = &
this % fmu(gin) % data(:, gout) / norm
this % dist(gin) % data(:, gout) = &
this % dist(gin) % data(:, gout) / norm
end if
end do
end do
end subroutine scattdatatabular_init
!===============================================================================
! SCATTDATA_*_CALC_F Calculates the value of f given mu (and gin,gout pair)
!===============================================================================
pure function scattdatalegendre_calc_f(this, gin, gout, mu) result(f)
class(ScattDataLegendre), intent(in) :: this ! The ScattData to evaluate
integer, intent(in) :: gin ! Incoming Energy Group
integer, intent(in) :: gout ! Outgoing Energy Group
real(8), intent(in) :: mu ! Angle of interest
real(8) :: f ! Return value of f(mu)
! Plug mu in to the legendre expansion and go from there
if (gout < this % gmin(gin) .or. gout > this % gmax(gin)) then
f = ZERO
else
f = evaluate_legendre(this % dist(gin) % data(:, gout), mu)
end if
end function scattdatalegendre_calc_f
pure function scattdatahistogram_calc_f(this, gin, gout, mu) result(f)
class(ScattDataHistogram), intent(in) :: this ! The ScattData to evaluate
integer, intent(in) :: gin ! Incoming Energy Group
integer, intent(in) :: gout ! Outgoing Energy Group
real(8), intent(in) :: mu ! Angle of interest
real(8) :: f ! Return value of f(mu)
integer :: imu
if (gout < this % gmin(gin) .or. gout > this % gmax(gin)) then
f = ZERO
else
! Find mu bin
if (mu == ONE) then
imu = size(this % fmu(gin) % data, dim=1)
else
imu = floor((mu + ONE) / this % dmu + ONE)
end if
f = this % fmu(gin) % data(imu, gout)
end if
end function scattdatahistogram_calc_f
pure function scattdatatabular_calc_f(this, gin, gout, mu) result(f)
class(ScattDataTabular), intent(in) :: this ! The ScattData to evaluate
integer, intent(in) :: gin ! Incoming Energy Group
integer, intent(in) :: gout ! Outgoing Energy Group
real(8), intent(in) :: mu ! Angle of interest
real(8) :: f ! Return value of f(mu)
integer :: imu
real(8) :: r
if (gout < this % gmin(gin) .or. gout > this % gmax(gin)) then
f = ZERO
else
! Find mu bin
if (mu == ONE) then
imu = size(this % fmu(gin) % data, dim=1) - 1
else
imu = floor((mu + ONE) / this % dmu + ONE)
end if
! Now interpolate to find f(mu)
r = (mu - this % mu(imu)) / (this % mu(imu + 1) - this % mu(imu))
f = (ONE - r) * this % fmu(gin) % data(imu, gout) + &
r * this % fmu(gin) % data(imu + 1, gout)
end if
end function scattdatatabular_calc_f
!===============================================================================
! SCATTDATA*_SCATTER Samples the outgoing energy and change in angle.
!===============================================================================
subroutine scattdatalegendre_sample(this, gin, gout, mu, wgt)
class(ScattDataLegendre), intent(in) :: this ! Scattering object to use
integer, intent(in) :: gin ! Incoming neutron group
integer, intent(out) :: gout ! Sampled outgoin group
real(8), intent(out) :: mu ! Sampled change in angle
real(8), intent(inout) :: wgt ! Particle weight
real(8) :: xi ! Our random number
real(8) :: prob ! Running probability
real(8) :: u, f, M
integer :: samples
xi = prn()
gout = this % gmin(gin)
prob = this % energy(gin) % data(gout)
do while ((prob < xi) .and. (gout < this % gmax(gin)))
gout = gout + 1
prob = prob + this % energy(gin) % data(gout)
end do
! Now we can sample mu using the legendre representation of the scattering
! kernel in data(1:this % order)
! Do with rejection sampling from a rectangular bounding box
! Set maximal value
M = this % max_val(gin) % data(gout)
samples = 0
do
mu = TWO * prn() - ONE
f = this % calc_f(gin, gout, mu)
if (f > ZERO) then
u = prn() * M
if (u <= f) then
exit
end if
end if
samples = samples + 1
if (samples > MAX_SAMPLE) then
call fatal_error("Maximum number of Legendre expansion samples reached!")
end if
end do
wgt = wgt * this % mult(gin) % data(gout)
end subroutine scattdatalegendre_sample
subroutine scattdatahistogram_sample(this, gin, gout, mu, wgt)
class(ScattDataHistogram), intent(in) :: this ! Scattering object to use
integer, intent(in) :: gin ! Incoming neutron group
integer, intent(out) :: gout ! Sampled outgoin group
real(8), intent(out) :: mu ! Sampled change in angle
real(8), intent(inout) :: wgt ! Particle weight
real(8) :: xi ! Our random number
real(8) :: prob ! Running probability
integer :: imu
xi = prn()
gout = this % gmin(gin)
prob = this % energy(gin) % data(gout)
do while ((prob < xi) .and. (gout < this % gmax(gin)))
gout = gout + 1
prob = prob + this % energy(gin) % data(gout)
end do
xi = prn()
if (xi < this % dist(gin) % data(1, gout)) then
imu = 1
else
imu = binary_search(this % dist(gin) % data(:, gout), &
size(this % dist(gin) % data(:, gout)), xi) + 1
end if
! Randomly select a mu in this bin.
mu = prn() * this % dmu + this % mu(imu)
wgt = wgt * this % mult(gin) % data(gout)
end subroutine scattdatahistogram_sample
subroutine scattdatatabular_sample(this, gin, gout, mu, wgt)
class(ScattDataTabular), intent(in) :: this ! Scattering object to use
integer, intent(in) :: gin ! Incoming neutron group
integer, intent(out) :: gout ! Sampled outgoin group
real(8), intent(out) :: mu ! Sampled change in angle
real(8), intent(inout) :: wgt ! Particle weight
real(8) :: xi ! Our random number
real(8) :: prob ! Running probability
real(8) :: mu0, frac, mu1
real(8) :: c_k, c_k1, p0, p1
integer :: k, NP
xi = prn()
gout = this % gmin(gin)
prob = this % energy(gin) % data(gout)
do while ((prob < xi) .and. (gout < this % gmax(gin)))
gout = gout + 1
prob = prob + this % energy(gin) % data(gout)
end do
! determine outgoing cosine bin
NP = size(this % dist(gin) % data(:, gout))
xi = prn()
c_k = this % dist(gin) % data(1, gout)
do k = 1, NP - 1
c_k1 = this % dist(gin) % data(k + 1, gout)
if (xi < c_k1) exit
c_k = c_k1
end do
! check to make sure k is <= NP - 1
k = min(k, NP - 1)
p0 = this % fmu(gin) % data(k, gout)
mu0 = this % mu(k)
! Linear-linear interpolation to find mu value w/in bin.
p1 = this % fmu(gin) % data(k + 1, gout)
mu1 = this % mu(k + 1)
if (p0 == p1) then
mu = mu0 + (xi - c_k) / p0
else
frac = (p1 - p0) / (mu1 - mu0)
mu = mu0 + &
(sqrt(max(ZERO, p0 * p0 + TWO * frac * (xi - c_k))) - p0) / frac
end if
if (mu <= -ONE) then
mu = -ONE
else if (mu >= ONE) then
mu = ONE
end if
wgt = wgt * this % mult(gin) % data(gout)
end subroutine scattdatatabular_sample
!===============================================================================
! SCATTDATA*_GET_MATRIX Reproduces the original scattering matrix (densely)
! using ScattData's information of fmu/dist, energy, and scattxs
!===============================================================================
subroutine scattdata_get_matrix(this, req_order, matrix)
class(ScattData), intent(in) :: this ! Scattering Object to work with
integer, intent(in) :: req_order ! Requested order of matrix
type(Jagged2D), allocatable, intent(inout) :: matrix(:) ! Resultant matrix just built
integer :: order, groups, gin, gout
groups = size(this % energy)
! Set gin and gout for getting the order
order = min(req_order, size(this % dist(1) % data, dim=1))
if (allocated(matrix)) deallocate(matrix)
allocate(matrix(groups))
! Initialize to 0; this way the zero entries in the dense matrix dont
! need to be explicitly set, requiring a significant increase in the
! lines of code.
do gin = 1, groups
allocate(matrix(gin) % data(order, groups))
do gout = this % gmin(gin), this % gmax(gin)
matrix(gin) % data(:, gout) = this % scattxs(gin) * &
this % energy(gin) % data(gout) * &
this % dist(gin) % data(1:order, gout)
end do
end do
end subroutine scattdata_get_matrix
subroutine scattdatahistogram_get_matrix(this, req_order, matrix)
class(ScattDataHistogram), intent(in) :: this ! Scattering Object to work with
integer, intent(in) :: req_order ! Requested order of matrix
type(Jagged2D), allocatable, intent(inout) :: matrix(:) ! Resultant matrix just built
integer :: order, groups, gin, gout
groups = size(this % energy)
order = min(req_order, size(this % dist(1) % data, dim=1))
if (allocated(matrix)) deallocate(matrix)
allocate(matrix(groups))
! Initialize to 0; this way the zero entries in the dense matrix dont
! need to be explicitly set, requiring a significant increase in the
! lines of code.
do gin = 1, groups
allocate(matrix(gin) % data(order, groups))
do gout = this % gmin(gin), this % gmax(gin)
matrix(gin) % data(:, gout) = this % scattxs(gin) * &
this % energy(gin) % data(gout) * &
this % fmu(gin) % data(1:order, gout)
end do
end do
end subroutine scattdatahistogram_get_matrix
subroutine scattdatatabular_get_matrix(this, req_order, matrix)
class(ScattDataTabular), intent(in) :: this ! Scattering Object to work with
integer, intent(in) :: req_order ! Requested order of matrix
type(Jagged2D), allocatable, intent(inout) :: matrix(:) ! Resultant matrix just built
integer :: order, groups, gin, gout
groups = size(this % energy)
order = min(req_order, size(this % dist(1) % data, dim=1))
if (allocated(matrix)) deallocate(matrix)
allocate(matrix(groups))
! Initialize to 0; this way the zero entries in the dense matrix dont
! need to be explicitly set, requiring a significant increase in the
! lines of code.
do gin = 1, groups
allocate(matrix(gin) % data(order, groups))
do gout = this % gmin(gin), this % gmax(gin)
matrix(gin) % data(:, gout) = this % scattxs(gin) * &
this % energy(gin) % data(gout) * &
this % fmu(gin) % data(1:order, gout)
end do
end do
end subroutine scattdatatabular_get_matrix
!===============================================================================
! JAGGED_FROM_DENSE_*D Creates a jagged array from a sparse dense matrix.
! The user can supply a key which indicates the values to remove, but the
! default is ZERO
!===============================================================================
subroutine jagged_from_dense_1D(dense, jagged, lo_bounds_, hi_bounds_, key_)
real(8), intent(in) :: dense(:, :)
type(Jagged1D), allocatable, intent(inout) :: jagged(:)
real(8), intent(in), optional :: key_
integer, intent(inout), allocatable, optional :: lo_bounds_(:)
integer, intent(inout), allocatable, optional :: hi_bounds_(:)
real(8) :: key
integer :: i, jmin, jmax
integer, allocatable :: lo_bounds(:), hi_bounds(:)
if (present(key_)) then
key = key_
else
key = ZERO
end if
allocate(lo_bounds(size(dense, dim=2)))
allocate(hi_bounds(size(dense, dim=2)))
if (allocated(jagged)) deallocate(jagged)
allocate(jagged(size(dense, dim=2)))
do i = 1, size(dense, dim=2)
! Find the min and max j values
do jmin = 1, size(dense, dim=1)
if (dense(jmin, i) /= key) exit
end do
do jmax = size(dense, dim=1), 1, -1
if (dense(jmax, i) /= key) exit
end do
! Treat the case of all values matching the key
if (jmin > jmax) then
jmin = i
jmax = i
end if
! Now store the jagged row
allocate(jagged(i) % data(jmin:jmax))
jagged(i) % data(jmin:jmax) = dense(jmin:jmax, i)
lo_bounds(i) = jmin
hi_bounds(i) = jmax
end do
if (present(lo_bounds_)) then
if (allocated(lo_bounds_)) deallocate(lo_bounds_)
allocate(lo_bounds_(size(dense, dim=2)))
lo_bounds_ = lo_bounds
end if
if (present(hi_bounds_)) then
if (allocated(hi_bounds_)) deallocate(hi_bounds_)
allocate(hi_bounds_(size(dense, dim=2)))
hi_bounds_ = hi_bounds
end if
end subroutine jagged_from_dense_1D
subroutine jagged_from_dense_2D(dense, jagged, lo_bounds_, hi_bounds_, key_)
real(8), intent(in) :: dense(:, :, :)
type(Jagged2D), allocatable, intent(inout) :: jagged(:)
real(8), intent(in), optional :: key_
integer, intent(inout), allocatable, optional :: lo_bounds_(:)
integer, intent(inout), allocatable, optional :: hi_bounds_(:)
real(8) :: key
integer :: i, jmin, jmax
integer, allocatable :: lo_bounds(:), hi_bounds(:)
if (present(key_)) then
key = key_
else
key = ZERO
end if
allocate(lo_bounds(size(dense, dim=3)))
allocate(hi_bounds(size(dense, dim=3)))
if (allocated(jagged)) deallocate(jagged)
allocate(jagged(size(dense, dim=3)))
do i = 1, size(dense, dim=3)
! Find the min and max j values
do jmin = 1, size(dense, dim=2)
if (any(dense(:, jmin, i) /= key)) exit
end do
do jmax = size(dense, dim=2), 1, -1
if (any(dense(:, jmax, i) /= key)) exit
end do
! Treat the case of all values matching the key
if (jmin > jmax) then
jmin = i
jmax = i
end if
! Now store the jagged row
allocate(jagged(i) % data(size(dense, dim=1), jmin:jmax))
jagged(i) % data(:, jmin:jmax) = dense(:, jmin:jmax, i)
lo_bounds(i) = jmin
hi_bounds(i) = jmax
end do
if (present(lo_bounds_)) then
if (allocated(lo_bounds_)) deallocate(lo_bounds_)
allocate(lo_bounds_(size(dense, dim=3)))
lo_bounds_ = lo_bounds
end if
if (present(hi_bounds_)) then
if (allocated(hi_bounds_)) deallocate(hi_bounds_)
allocate(hi_bounds_(size(dense, dim=3)))
hi_bounds_ = hi_bounds
end if
end subroutine jagged_from_dense_2D
end module scattdata_header

14
src/settings.F90
View file

@ -18,9 +18,9 @@ module settings
logical :: urr_ptables_on = .true.
! Default temperature and method for choosing temperatures
integer :: temperature_method = TEMPERATURE_NEAREST
integer(C_INT) :: temperature_method = TEMPERATURE_NEAREST
logical :: temperature_multipole = .false.
real(8) :: temperature_tolerance = 10.0_8
real(C_DOUBLE) :: temperature_tolerance = 10.0_8
real(8) :: temperature_default = 293.6_8
real(8) :: temperature_range(2) = [ZERO, ZERO]
@ -30,13 +30,16 @@ module settings
! MULTI-GROUP CROSS SECTION RELATED VARIABLES
! Maximum Data Order
integer :: max_order
integer(C_INT) :: max_order
! Whether or not to convert Legendres to tabulars
logical :: legendre_to_tabular = .true.
! Number of points to use in the Legendre to tabular conversion
integer :: legendre_to_tabular_points = 33
integer(C_INT) :: legendre_to_tabular_points = C_NONE
! ============================================================================
! SIMULATION VARIABLES
! Assume all tallies are spatially distinct
logical :: assume_separate = .false.
@ -44,9 +47,6 @@ module settings
! Use confidence intervals for results instead of standard deviations
logical :: confidence_intervals = .false.
! ============================================================================
! SIMULATION VARIABLES
integer(C_INT64_T), bind(C) :: n_particles = 0 ! # of particles per generation
integer(C_INT32_T), bind(C) :: n_batches ! # of batches
integer(C_INT32_T), bind(C) :: n_inactive ! # of inactive batches

2
src/simulation.F90
View file

@ -20,7 +20,7 @@ module simulation
use geometry_header, only: n_cells
use material_header, only: n_materials, materials
use message_passing
use mgxs_header, only: energy_bins, energy_bin_avg
use mgxs_interface, only: energy_bins, energy_bin_avg
use nuclide_header, only: micro_xs, n_nuclides
use output, only: header, print_columns, &
print_batch_keff, print_generation, print_runtime, &

2
src/source.F90
View file

@ -14,7 +14,7 @@ module source
use hdf5_interface
use math
use message_passing, only: rank
use mgxs_header, only: rev_energy_bins, num_energy_groups
use mgxs_interface, only: rev_energy_bins, num_energy_groups
use output, only: write_message
use particle_header, only: Particle
use random_lcg, only: prn, set_particle_seed, prn_set_stream

11
src/state_point.F90
View file

@ -22,7 +22,7 @@ module state_point
use hdf5_interface
use mesh_header, only: RegularMesh, meshes, n_meshes
use message_passing
use mgxs_header, only: nuclides_MG
use mgxs_interface
use nuclide_header, only: nuclides
use output, only: time_stamp
use random_lcg, only: openmc_get_seed, openmc_set_seed
@ -73,6 +73,7 @@ contains
character(MAX_WORD_LEN), allocatable :: str_array(:)
character(C_CHAR), pointer :: string(:)
character(len=:, kind=C_CHAR), allocatable :: filename_
character(MAX_WORD_LEN, kind=C_CHAR) :: temp_name
err = 0
if (present(filename)) then
@ -307,11 +308,13 @@ contains
str_array(j) = nuclides(tally % nuclide_bins(j)) % name
end if
else
i_xs = index(nuclides_MG(tally % nuclide_bins(j)) % obj % name, '.')
call get_name_c(tally % nuclide_bins(j), len(temp_name), &
temp_name)
i_xs = index(temp_name, '.')
if (i_xs > 0) then
str_array(j) = nuclides_MG(tally % nuclide_bins(j)) % obj % name(1 : i_xs-1)
str_array(j) = trim(temp_name(1 : i_xs-1))
else
str_array(j) = nuclides_MG(tally % nuclide_bins(j)) % obj % name
str_array(j) = trim(temp_name)
end if
end if
else

30
src/string_functions.cpp Normal file
View file

@ -0,0 +1,30 @@
#include "string_functions.h"
namespace openmc {
std::string& strtrim(std::string& s)
{
const char* t = " \t\n\r\f\v";
s.erase(s.find_last_not_of(t) + 1);
s.erase(0, s.find_first_not_of(t));
return s;
}
char* strtrim(char* c_str)
{
std::string std_str;
std_str.assign(c_str);
strtrim(std_str);
int length = std_str.copy(c_str, std_str.size());
c_str[length] = '\0';
return c_str;
}
void to_lower(std::string& str)
{
for (int i = 0; i < str.size(); i++) str[i] = std::tolower(str[i]);
}
} // namespace openmc

18
src/string_functions.h Normal file
View file

@ -0,0 +1,18 @@
//! \file string_functions.h
//! A collection of helper routines for C-strings and STL strings
#ifndef STRING_FUNCTIONS_H
#define STRING_FUNCTIONS_H
#include <string>
namespace openmc {
std::string& strtrim(std::string& s);
char* strtrim(char* c_str);
void to_lower(std::string& str);
} // namespace openmc
#endif // STRING_FUNCTIONS_H

24
src/summary.F90
View file

@ -8,7 +8,7 @@ module summary
use material_header, only: Material, n_materials
use mesh_header, only: RegularMesh
use message_passing
use mgxs_header, only: nuclides_MG
use mgxs_interface
use nuclide_header
use output, only: time_stamp
use settings, only: run_CE
@ -93,7 +93,7 @@ contains
num_nuclides = 0
num_macros = 0
do i = 1, n_nuclides
if (nuclides_MG(i) % obj % awr /= MACROSCOPIC_AWR) then
if (get_awr_c(i) /= MACROSCOPIC_AWR) then
num_nuclides = num_nuclides + 1
else
num_macros = num_macros + 1
@ -118,12 +118,14 @@ contains
nuc_names(i) = nuclides(i) % name
awrs(i) = nuclides(i) % awr
else
if (nuclides_MG(i) % obj % awr /= MACROSCOPIC_AWR) then
nuc_names(j) = nuclides_MG(i) % obj % name
awrs(j) = nuclides_MG(i) % obj % awr
if (get_awr_c(i) /= MACROSCOPIC_AWR) then
call get_name_c(i, len(nuc_names(j)), nuc_names(j))
nuc_names(j) = trim(nuc_names(j))
awrs(j) = get_awr_c(i)
j = j + 1
else
macro_names(k) = nuclides_MG(i) % obj % name
call get_name_c(i, len(macro_names(k)), macro_names(k))
macro_names(k) = trim(macro_names(k))
k = k + 1
end if
end if
@ -353,7 +355,7 @@ contains
num_nuclides = 0
num_macros = 0
do j = 1, m % n_nuclides
if (nuclides_MG(m % nuclide(j)) % obj % awr /= MACROSCOPIC_AWR) then
if (get_awr_c(m % nuclide(j)) /= MACROSCOPIC_AWR) then
num_nuclides = num_nuclides + 1
else
num_macros = num_macros + 1
@ -379,12 +381,14 @@ contains
k = 1
n = 1
do j = 1, m % n_nuclides
if (nuclides_MG(m % nuclide(j)) % obj % awr /= MACROSCOPIC_AWR) then
nuc_names(k) = nuclides_MG(m % nuclide(j)) % obj % name
if (get_awr_c(m % nuclide(j)) /= MACROSCOPIC_AWR) then
call get_name_c(m % nuclide(j), len(nuc_names(k)), nuc_names(k))
nuc_names(k) = trim(nuc_names(k))
nuc_densities(k) = m % atom_density(j)
k = k + 1
else
macro_names(n) = nuclides_MG(m % nuclide(j)) % obj % name
call get_name_c(m % nuclide(j), len(macro_names(n)), macro_names(n))
macro_names(n) = trim(macro_names(n))
n = n + 1
end if
end do

2
src/surface.h
View file

@ -6,7 +6,7 @@
#include <string>
#include "hdf5.h"
#include "pugixml/pugixml.hpp"
#include "pugixml.hpp"
#include "constants.h"

243
src/tallies/tally.F90
View file

@ -10,7 +10,7 @@ module tally
use math, only: t_percentile
use mesh_header, only: RegularMesh, meshes
use message_passing
use mgxs_header
use mgxs_interface
use nuclide_header
use output, only: header
use particle_header, only: LocalCoord, Particle
@ -1229,8 +1229,6 @@ contains
real(8) :: p_uvw(3) ! Particle's current uvw
integer :: p_g ! Particle group to use for getting info
! to tally with.
class(Mgxs), pointer :: matxs
class(Mgxs), pointer :: nucxs
! Set the direction and group to use with get_xs
if (t % estimator == ESTIMATOR_ANALOG .or. &
@ -1268,13 +1266,13 @@ contains
! To significantly reduce de-referencing, point matxs to the
! macroscopic Mgxs for the material of interest
matxs => macro_xs(p % material) % obj
call set_macro_angle_index_c(p % material, p_uvw)
! Do same for nucxs, point it to the microscopic nuclide data of interest
if (i_nuclide > 0) then
nucxs => nuclides_MG(i_nuclide) % obj
! And since we haven't calculated this temperature index yet, do so now
call nucxs % find_temperature(p % sqrtkT)
call set_nuclide_temperature_index_c(i_nuclide, p % sqrtkT)
call set_nuclide_angle_index_c(i_nuclide, p_uvw)
end if
i = 0
@ -1328,14 +1326,14 @@ contains
end if
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('total', p_g, UVW=p_uvw) / &
matxs % get_xs('total', p_g, UVW=p_uvw) * flux
score = score * flux * atom_density * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_TOTAL, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_TOTAL, p_g)
end if
else
if (i_nuclide > 0) then
score = nucxs % get_xs('total', p_g, UVW=p_uvw) * &
score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_TOTAL, p_g) * &
atom_density * flux
else
score = material_xs % total * flux
@ -1358,19 +1356,23 @@ contains
end if
if (i_nuclide > 0) then
score = score * nucxs % get_xs('inverse-velocity', p_g, UVW=p_uvw) &
/ matxs % get_xs('absorption', p_g, UVW=p_uvw) * flux
score = score * flux * get_nuclide_xs_c(i_nuclide, &
MG_GET_XS_INVERSE_VELOCITY, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
else
score = score * matxs % get_xs('inverse-velocity', p_g, UVW=p_uvw) &
/ matxs % get_xs('absorption', p_g, UVW=p_uvw) * flux
score = score * flux * get_macro_xs_c(p % material, &
MG_GET_XS_INVERSE_VELOCITY, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
else
if (i_nuclide > 0) then
score = flux * nucxs % get_xs('inverse-velocity', p_g, UVW=p_uvw)
score = flux * get_nuclide_xs_c(i_nuclide, &
MG_GET_XS_INVERSE_VELOCITY, p_g)
else
score = flux * matxs % get_xs('inverse-velocity', p_g, UVW=p_uvw)
score = flux * get_macro_xs_c(p % material, &
MG_GET_XS_INVERSE_VELOCITY, p_g)
end if
end if
@ -1392,21 +1394,23 @@ contains
! adjust the score by the actual probability for that nuclide.
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('scatter*f_mu/mult', p % last_g, p % g, &
UVW=p_uvw, MU=p % mu) / &
matxs % get_xs('scatter*f_mu/mult', p % last_g, p % g, &
UVW=p_uvw, MU=p % mu)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_SCATTER_FMU_MULT, &
p % last_g, p % g, MU=p % mu) / &
get_macro_xs_c(p % material, MG_GET_XS_SCATTER_FMU_MULT, &
p % last_g, p % g, MU=p % mu)
end if
else
if (i_nuclide > 0) then
score = atom_density * flux * &
nucxs % get_xs('scatter/mult', p_g, UVW=p_uvw)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_SCATTER_MULT, &
p_g, MU=p % mu)
else
! Get the scattering x/s and take away
! the multiplication baked in to sigS
score = flux * &
matxs % get_xs('scatter/mult', p_g, UVW=p_uvw)
get_macro_xs_c(p % material, MG_GET_XS_SCATTER_MULT, &
p_g, MU=p % mu)
end if
end if
@ -1428,19 +1432,20 @@ contains
! adjust the score by the actual probability for that nuclide.
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('scatter*f_mu', p % last_g, p % g, &
UVW=p_uvw, MU=p % mu) / &
matxs % get_xs('scatter*f_mu', p % last_g, p % g, &
UVW=p_uvw, MU=p % mu)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_SCATTER_FMU, &
p % last_g, p % g, MU=p % mu) / &
get_macro_xs_c(p % material, MG_GET_XS_SCATTER_FMU, &
p % last_g, p % g, MU=p % mu)
end if
else
if (i_nuclide > 0) then
score = nucxs % get_xs('scatter', p_g, UVW=p_uvw) * &
atom_density * flux
score = atom_density * flux * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_SCATTER, p_g)
else
! Get the scattering x/s, which includes multiplication
score = matxs % get_xs('scatter', p_g, UVW=p_uvw) * flux
score = flux * &
get_macro_xs_c(p % material, MG_GET_XS_SCATTER, p_g)
end if
end if
@ -1460,13 +1465,13 @@ contains
end if
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('absorption', p_g, UVW=p_uvw) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_ABSORPTION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
else
if (i_nuclide > 0) then
score = nucxs % get_xs('absorption', p_g, UVW=p_uvw) * &
atom_density * flux
score = atom_density * flux * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_ABSORPTION, p_g)
else
score = material_xs % absorption * flux
end if
@ -1491,19 +1496,19 @@ contains
end if
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('fission', p_g, UVW=p_uvw) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
else
score = score * &
matxs % get_xs('fission', p_g, UVW=p_uvw) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
else
if (i_nuclide > 0) then
score = nucxs % get_xs('fission', p_g, UVW=p_uvw) * &
score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) * &
atom_density * flux
else
score = matxs % get_xs('fission', p_g, UVW=p_uvw) * flux
score = get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g) * flux
end if
end if
@ -1529,12 +1534,12 @@ contains
score = p % absorb_wgt * flux
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('nu-fission', p_g, UVW=p_uvw) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_NU_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
else
score = score * &
matxs % get_xs('nu-fission', p_g, UVW=p_uvw) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
get_macro_xs_c(p % material, MG_GET_XS_NU_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
else
! Skip any non-fission events
@ -1547,17 +1552,17 @@ contains
score = keff * p % wgt_bank * flux
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('fission', p_g, UVW=p_uvw) / &
matxs % get_xs('fission', p_g, UVW=p_uvw)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g)
end if
end if
else
if (i_nuclide > 0) then
score = nucxs % get_xs('nu-fission', p_g, UVW=p_uvw) * &
score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_NU_FISSION, p_g) * &
atom_density * flux
else
score = matxs % get_xs('nu-fission', p_g, UVW=p_uvw) * flux
score = get_macro_xs_c(p % material, MG_GET_XS_NU_FISSION, p_g) * flux
end if
end if
@ -1583,12 +1588,12 @@ contains
score = p % absorb_wgt * flux
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('prompt-nu-fission', p_g, UVW=p_uvw) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_PROMPT_NU_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
else
score = score * &
matxs % get_xs('prompt-nu-fission', p_g, UVW=p_uvw) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
get_macro_xs_c(p % material, MG_GET_XS_PROMPT_NU_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
else
! Skip any non-fission events
@ -1602,17 +1607,17 @@ contains
/ real(p % n_bank, 8)) * flux
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('fission', p_g, UVW=p_uvw) / &
matxs % get_xs('fission', p_g, UVW=p_uvw)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g)
end if
end if
else
if (i_nuclide > 0) then
score = nucxs % get_xs('prompt-nu-fission', p_g, UVW=p_uvw) * &
score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_PROMPT_NU_FISSION, p_g) * &
atom_density * flux
else
score = matxs % get_xs('prompt-nu-fission', p_g, UVW=p_uvw) * flux
score = get_macro_xs_c(p % material, MG_GET_XS_PROMPT_NU_FISSION, p_g) * flux
end if
end if
@ -1638,7 +1643,7 @@ contains
! No fission events occur if survival biasing is on -- need to
! calculate fraction of absorptions that would have resulted in
! nu-fission
if (matxs % get_xs('absorption', p_g, UVW=p_uvw) > ZERO) then
if (get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) > ZERO) then
if (dg_filter > 0) then
select type(filt => filters(t % filter(dg_filter)) % obj)
@ -1653,13 +1658,13 @@ contains
score = p % absorb_wgt * flux
if (i_nuclide > 0) then
score = score * nucxs % get_xs('delayed-nu-fission', &
p_g, UVW=p_uvw, dg=d) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
score = score * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
else
score = score * matxs % get_xs('delayed-nu-fission', &
p_g, UVW=p_uvw, dg=d) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
score = score * &
get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
call score_fission_delayed_dg(t, d_bin, score, score_index)
@ -1669,11 +1674,13 @@ contains
else
score = p % absorb_wgt * flux
if (i_nuclide > 0) then
score = score * nucxs % get_xs('delayed-nu-fission', p_g, &
UVW=p_uvw) / matxs % get_xs('absorption', p_g, UVW=p_uvw)
score = score * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
else
score = score * matxs % get_xs('delayed-nu-fission', p_g, &
UVW=p_uvw) / matxs % get_xs('absorption', p_g, UVW=p_uvw)
score = score * &
get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
end if
end if
@ -1703,8 +1710,8 @@ contains
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('fission', p_g, UVW=p_uvw) / &
matxs % get_xs('fission', p_g, UVW=p_uvw)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g)
end if
call score_fission_delayed_dg(t, d_bin, score, score_index)
@ -1715,8 +1722,8 @@ contains
score = keff * p % wgt_bank / p % n_bank * sum(p % n_delayed_bank) * flux
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('fission', p_g, UVW=p_uvw) / &
matxs % get_xs('fission', p_g, UVW=p_uvw)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g)
end if
end if
end if
@ -1735,11 +1742,11 @@ contains
d = filt % groups(d_bin)
if (i_nuclide > 0) then
score = nucxs % get_xs('delayed-nu-fission', p_g, &
UVW=p_uvw, dg=d) * atom_density * flux
score = atom_density * flux * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
else
score = matxs % get_xs('delayed-nu-fission', p_g, &
UVW=p_uvw, dg=d) * flux
score = flux * &
get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
end if
call score_fission_delayed_dg(t, d_bin, score, score_index)
@ -1748,11 +1755,12 @@ contains
end select
else
if (i_nuclide > 0) then
score = nucxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw) &
* atom_density * flux
score = atom_density * flux * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g)
else
score = matxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw) &
* flux
score = flux * &
get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g)
end if
end if
end if
@ -1767,7 +1775,7 @@ contains
! No fission events occur if survival biasing is on -- need to
! calculate fraction of absorptions that would have resulted in
! nu-fission
if (matxs % get_xs('absorption', p_g, UVW=p_uvw) > ZERO) then
if (get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) > ZERO) then
if (dg_filter > 0) then
select type(filt => filters(t % filter(dg_filter)) % obj)
@ -1782,17 +1790,15 @@ contains
score = p % absorb_wgt * flux
if (i_nuclide > 0) then
score = score * nucxs % get_xs('decay rate', p_g, &
UVW=p_uvw, dg=d) * &
nucxs % get_xs('delayed-nu-fission', p_g, &
UVW=p_uvw, dg=d) / matxs % get_xs('absorption', &
p_g, UVW=p_uvw)
score = score * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
else
score = score * matxs % get_xs('decay rate', p_g, &
UVW=p_uvw, dg=d) * &
matxs % get_xs('delayed-nu-fission', p_g, &
UVW=p_uvw, dg=d) / matxs % get_xs('absorption', &
p_g, UVW=p_uvw)
score = score * &
get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
call score_fission_delayed_dg(t, d_bin, score, score_index)
@ -1809,15 +1815,15 @@ contains
! for all delayed groups.
do d = 1, num_delayed_groups
if (i_nuclide > 0) then
score = score + p % absorb_wgt * &
nucxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * &
nucxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, &
dg=d) / matxs % get_xs('absorption', p_g, UVW=p_uvw) * flux
score = score + p % absorb_wgt * flux * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
else
score = score + p % absorb_wgt * &
matxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * &
matxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, &
dg=d) / matxs % get_xs('absorption', p_g, UVW=p_uvw) * flux
score = score + p % absorb_wgt * flux * &
get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
end do
end if
@ -1846,13 +1852,13 @@ contains
if (i_nuclide > 0) then
score = score + keff * atom_density * &
fission_bank(n_bank - p % n_bank + k) % wgt * &
nucxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=g) * &
nucxs % get_xs('fission', p_g, UVW=p_uvw) / &
matxs % get_xs('fission', p_g, UVW=p_uvw) * flux
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g) * flux
else
score = score + keff * &
fission_bank(n_bank - p % n_bank + k) % wgt * &
matxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=g) * flux
get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * flux
end if
! if the delayed group filter is present, tally to corresponding
@ -1904,13 +1910,13 @@ contains
d = filt % groups(d_bin)
if (i_nuclide > 0) then
score = nucxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * &
nucxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, &
dg=d) * atom_density * flux
score = atom_density * flux * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
else
score = matxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * &
matxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, &
dg=d) * flux
score = flux * &
get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
end if
call score_fission_delayed_dg(t, d_bin, score, score_index)
@ -1927,12 +1933,12 @@ contains
do d = 1, num_delayed_groups
if (i_nuclide > 0) then
score = score + atom_density * flux * &
nucxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * &
nucxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, dg=d)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
else
score = score + flux * &
matxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * &
matxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, dg=d)
get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
end if
end do
end if
@ -1957,19 +1963,20 @@ contains
end if
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('kappa-fission', p_g, UVW=p_uvw) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_KAPPA_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
else
score = score * &
matxs % get_xs('kappa-fission', p_g, UVW=p_uvw) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
get_macro_xs_c(p % material, MG_GET_XS_KAPPA_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
else
if (i_nuclide > 0) then
score = nucxs % get_xs('kappa-fission', p_g, UVW=p_uvw) * &
score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_KAPPA_FISSION, p_g) * &
atom_density * flux
else
score = matxs % get_xs('kappa-fission', p_g, UVW=p_uvw) * flux
score = flux * &
get_macro_xs_c(p % material, MG_GET_XS_KAPPA_FISSION, p_g)
end if
end if
@ -1988,8 +1995,6 @@ contains
t % results(RESULT_VALUE, score_index, filter_index) + score
end do SCORE_LOOP
nullify(matxs, nucxs)
end subroutine score_general_mg
!===============================================================================

2
src/tallies/tally_filter_energy.F90
View file

@ -6,7 +6,7 @@ module tally_filter_energy
use constants
use error
use hdf5_interface
use mgxs_header, only: num_energy_groups, rev_energy_bins
use mgxs_interface, only: num_energy_groups, rev_energy_bins
use particle_header, only: Particle
use settings, only: run_CE
use string, only: to_str

10
src/tracking.F90
View file

@ -1,5 +1,7 @@
module tracking
use, intrinsic :: ISO_C_BINDING
use constants
use error, only: warning, write_message
use geometry_header, only: cells
@ -7,7 +9,7 @@ module tracking
check_cell_overlap
use material_header, only: materials, Material
use message_passing
use mgxs_header
use mgxs_interface
use nuclide_header
use particle_header, only: LocalCoord, Particle
use physics, only: collision
@ -112,8 +114,10 @@ contains
end if
else
! Get the MG data
call macro_xs(p % material) % obj % calculate_xs(p % g, p % sqrtkT, &
p % coord(p % n_coord) % uvw, material_xs)
call calculate_xs_c(p % material, p % g, p % sqrtkT, &
p % coord(p % n_coord) % uvw, material_xs % total, &
material_xs % absorption, material_xs % nu_fission)
! Finally, update the particle group while we have already checked
! for if multi-group

2
src/xml_interface.h
View file

@ -4,7 +4,7 @@
#include <string>
#include <vector>
#include "pugixml/pugixml.hpp"
#include "pugixml.hpp"
namespace openmc {

713
src/xsdata.cpp Normal file
View file

@ -0,0 +1,713 @@
#include <cmath>
#include <cstdlib>
#include <algorithm>
#include <numeric>
#include "constants.h"
#include "error.h"
#include "math_functions.h"
#include "random_lcg.h"
#include "xsdata.h"
namespace openmc {
//==============================================================================
// XsData class methods
//==============================================================================
XsData::XsData(int energy_groups, int num_delayed_groups, bool fissionable,
int scatter_format, int n_pol, int n_azi)
{
int n_ang = n_pol * n_azi;
// check to make sure scatter format is OK before we allocate
if (scatter_format != ANGLE_HISTOGRAM && scatter_format != ANGLE_TABULAR &&
scatter_format != ANGLE_LEGENDRE) {
fatal_error("Invalid scatter_format!");
}
// allocate all [temperature][phi][theta][in group] quantities
total = double_2dvec(n_ang, double_1dvec(energy_groups, 0.));
absorption = double_2dvec(n_ang, double_1dvec(energy_groups, 0.));
inverse_velocity = double_2dvec(n_ang, double_1dvec(energy_groups, 0.));
if (fissionable) {
fission = double_2dvec(n_ang, double_1dvec(energy_groups, 0.));
nu_fission = double_2dvec(n_ang, double_1dvec(energy_groups, 0.));
prompt_nu_fission = double_2dvec(n_ang, double_1dvec(energy_groups, 0.));
kappa_fission = double_2dvec(n_ang, double_1dvec(energy_groups, 0.));
}
// allocate decay_rate; [temperature][phi][theta][delayed group]
decay_rate = double_2dvec(n_ang, double_1dvec(num_delayed_groups, 0.));
if (fissionable) {
// allocate delayed_nu_fission; [temperature][phi][theta][in group][delay group]
delayed_nu_fission = double_3dvec(n_ang, double_2dvec(energy_groups,
double_1dvec(num_delayed_groups, 0.)));
// chi_prompt; [temperature][phi][theta][in group][delayed group]
chi_prompt = double_3dvec(n_ang, double_2dvec(energy_groups,
double_1dvec(energy_groups, 0.)));
// chi_delayed; [temperature][phi][theta][in group][out group][delay group]
chi_delayed = double_4dvec(n_ang, double_3dvec(energy_groups,
double_2dvec(energy_groups, double_1dvec(num_delayed_groups, 0.))));
}
for (int a = 0; a < n_ang; a++) {
if (scatter_format == ANGLE_HISTOGRAM) {
// scatter[a] = std::make_unique(ScattDataHistogram);
scatter.emplace_back(new ScattDataHistogram);
} else if (scatter_format == ANGLE_TABULAR) {
// scatter[a] = std::make_unique(ScattDataTabular);
scatter.emplace_back(new ScattDataTabular);
} else if (scatter_format == ANGLE_LEGENDRE) {
// scatter[a] = std::make_unique(ScattDataLegendre);
scatter.emplace_back(new ScattDataLegendre);
}
}
}
//==============================================================================
void
XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format,
int final_scatter_format, int order_data, int max_order,
int legendre_to_tabular_points, bool is_isotropic, int n_pol, int n_azi)
{
// Reconstruct the dimension information so it doesn't need to be passed
int n_ang = n_pol * n_azi;
int energy_groups = total[0].size();
int delayed_groups = decay_rate[0].size();
// Set the fissionable-specific data
if (fissionable) {
fission_from_hdf5(xsdata_grp, n_pol, n_azi, energy_groups, delayed_groups,
is_isotropic);
}
// Get the non-fission-specific data
read_nd_vector(xsdata_grp, "decay_rate", decay_rate);
read_nd_vector(xsdata_grp, "absorption", absorption, true);
read_nd_vector(xsdata_grp, "inverse-velocity", inverse_velocity);
// Get scattering data
scatter_from_hdf5(xsdata_grp, n_pol, n_azi, energy_groups, scatter_format,
final_scatter_format, order_data, max_order, legendre_to_tabular_points);
// Check absorption to ensure it is not 0 since it is often the
// denominator in tally methods
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
if (absorption[a][gin] == 0.) absorption[a][gin] = 1.e-10;
}
}
// Get or calculate the total x/s
if (object_exists(xsdata_grp, "total")) {
read_nd_vector(xsdata_grp, "total", total);
} else {
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
total[a][gin] = absorption[a][gin] + scatter[a]->scattxs[gin];
}
}
}
// Fix if total is 0, since it is in the denominator when tallying
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
if (total[a][gin] == 0.) total[a][gin] = 1.e-10;
}
}
}
//==============================================================================
void
XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
int energy_groups, int delayed_groups, bool is_isotropic)
{
int n_ang = n_pol * n_azi;
// Get the fission and kappa_fission data xs; these are optional
read_nd_vector(xsdata_grp, "fission", fission);
read_nd_vector(xsdata_grp, "kappa-fission", kappa_fission);
// Set/get beta
double_3dvec temp_beta =double_3dvec(n_ang, double_2dvec(energy_groups,
double_1dvec(delayed_groups, 0.)));
if (object_exists(xsdata_grp, "beta")) {
hid_t xsdata = open_dataset(xsdata_grp, "beta");
int ndims = dataset_ndims(xsdata);
// raise ndims to make the isotropic ndims the same as angular
if (is_isotropic) ndims += 2;
if (ndims == 3) {
// Beta is input as [delayed group]
double_1dvec temp_arr(n_pol * n_azi * delayed_groups);
read_nd_vector(xsdata_grp, "beta", temp_arr);
// Broadcast to all incoming groups
int temp_idx = 0;
for (int a = 0; a < n_ang; a++) {
for (int dg = 0; dg < delayed_groups; dg++) {
// Set the first group index and copy the rest
temp_beta[a][0][dg] = temp_arr[temp_idx++];
for (int gin = 1; gin < energy_groups; gin++) {
temp_beta[a][gin] = temp_beta[a][0];
}
}
}
} else if (ndims == 4) {
// Beta is input as [in group][delayed group]
read_nd_vector(xsdata_grp, "beta", temp_beta);
} else {
fatal_error("beta must be provided as a 3D or 4D array!");
}
}
// If chi is provided, set chi-prompt and chi-delayed
if (object_exists(xsdata_grp, "chi")) {
double_2dvec temp_arr(n_ang, double_1dvec(energy_groups));
read_nd_vector(xsdata_grp, "chi", temp_arr);
for (int a = 0; a < n_ang; a++) {
// First set the first group
for (int gout = 0; gout < energy_groups; gout++) {
chi_prompt[a][0][gout] = temp_arr[a][gout];
}
// Now normalize this data
double chi_sum = std::accumulate(chi_prompt[a][0].begin(),
chi_prompt[a][0].end(),
0.);
if (chi_sum <= 0.) {
fatal_error("Encountered chi for a group that is <= 0!");
}
for (int gout = 0; gout < energy_groups; gout++) {
chi_prompt[a][0][gout] /= chi_sum;
}
// And extend to the remaining incoming groups
for (int gin = 1; gin < energy_groups; gin++) {
chi_prompt[a][gin] = chi_prompt[a][0];
}
// Finally set chi-delayed equal to chi-prompt
// Set chi-delayed to chi-prompt
for(int gin = 0; gin < energy_groups; gin++) {
for (int gout = 0; gout < energy_groups; gout++) {
for (int dg = 0; dg < delayed_groups; dg++) {
chi_delayed[a][gin][gout][dg] =
chi_prompt[a][gin][gout];
}
}
}
}
}
// If nu-fission is provided, set prompt- and delayed-nu-fission;
// if nu-fission is a matrix, set chi-prompt and chi-delayed.
if (object_exists(xsdata_grp, "nu-fission")) {
hid_t xsdata = open_dataset(xsdata_grp, "nu-fission");
int ndims = dataset_ndims(xsdata);
// raise ndims to make the isotropic ndims the same as angular
if (is_isotropic) ndims += 2;
if (ndims == 3) {
// nu-fission is a 3-d array
read_nd_vector(xsdata_grp, "nu-fission", prompt_nu_fission);
// set delayed-nu-fission and correct prompt-nu-fission with beta
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
for (int dg = 0; dg < delayed_groups; dg++) {
delayed_nu_fission[a][gin][dg] =
temp_beta[a][gin][dg] * prompt_nu_fission[a][gin];
}
// Correct the prompt-nu-fission using the delayed neutron fraction
if (delayed_groups > 0) {
double beta_sum = std::accumulate(temp_beta[a][gin].begin(),
temp_beta[a][gin].end(), 0.);
prompt_nu_fission[a][gin] *= (1. - beta_sum);
}
}
}
} else if (ndims == 4) {
// nu-fission is a matrix
read_nd_vector(xsdata_grp, "nu-fission", chi_prompt);
// Normalize the chi info so the CDF is 1.
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
double chi_sum = std::accumulate(chi_prompt[a][gin].begin(),
chi_prompt[a][gin].end(), 0.);
// Set the vector nu-fission from the matrix nu-fission
prompt_nu_fission[a][gin] = chi_sum;
if (chi_sum >= 0.) {
for (int gout = 0; gout < energy_groups; gout++) {
chi_prompt[a][gin][gout] /= chi_sum;
}
} else {
fatal_error("Encountered chi for a group that is <= 0!");
}
}
// set chi-delayed to chi-prompt
for (int gin = 0; gin < energy_groups; gin++) {
for (int gout = 0; gout < energy_groups; gout++) {
for (int dg = 0; dg < delayed_groups; dg++) {
chi_delayed[a][gin][gout][dg] =
chi_prompt[a][gin][gout];
}
}
}
// Set the delayed-nu-fission and correct prompt-nu-fission with beta
for (int gin = 0; gin < energy_groups; gin++) {
for (int dg = 0; dg < delayed_groups; dg++) {
delayed_nu_fission[a][gin][dg] =
temp_beta[a][gin][dg] *
prompt_nu_fission[a][gin];
}
// Correct prompt-nu-fission using the delayed neutron fraction
if (delayed_groups > 0) {
double beta_sum = std::accumulate(temp_beta[a][gin].begin(),
temp_beta[a][gin].end(), 0.);
prompt_nu_fission[a][gin] *= (1. - beta_sum);
}
}
}
} else {
fatal_error("nu-fission must be provided as a 3D or 4D array!");
}
close_dataset(xsdata);
}
// If chi-prompt is provided, set chi-prompt
if (object_exists(xsdata_grp, "chi-prompt")) {
double_2dvec temp_arr(n_ang, double_1dvec(energy_groups));
read_nd_vector(xsdata_grp, "chi-prompt", temp_arr);
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
for (int gout = 0; gout < energy_groups; gout++) {
chi_prompt[a][gin][gout] = temp_arr[a][gout];
}
// Normalize chi so its CDF goes to 1
double chi_sum = std::accumulate(chi_prompt[a][gin].begin(),
chi_prompt[a][gin].end(), 0.);
if (chi_sum >= 0.) {
for (int gout = 0; gout < energy_groups; gout++) {
chi_prompt[a][gin][gout] /= chi_sum;
}
} else {
fatal_error("Encountered chi-prompt for a group that is <= 0.!");
}
}
}
}
// If chi-delayed is provided, set chi-delayed
if (object_exists(xsdata_grp, "chi-delayed")) {
hid_t xsdata = open_dataset(xsdata_grp, "chi-delayed");
int ndims = dataset_ndims(xsdata);
// raise ndims to make the isotropic ndims the same as angular
if (is_isotropic) ndims += 2;
close_dataset(xsdata);
if (ndims == 3) {
// chi-delayed is a [in group] vector
double_2dvec temp_arr(n_ang, double_1dvec(energy_groups));
read_nd_vector(xsdata_grp, "chi-delayed", temp_arr);
for (int a = 0; a < n_ang; a++) {
// normalize the chi CDF to 1
double chi_sum = std::accumulate(temp_arr[a].begin(),
temp_arr[a].end(), 0.);
if (chi_sum <= 0.) {
fatal_error("Encountered chi-delayed for a group that is <= 0!");
}
// set chi-delayed
for (int gin = 0; gin < energy_groups; gin++) {
for (int gout = 0; gout < energy_groups; gout++) {
for (int dg = 0; dg < delayed_groups; dg++) {
chi_delayed[a][gin][gout][dg] = temp_arr[a][gout] / chi_sum;
}
}
}
}
} else if (ndims == 4) {
// chi_delayed is a matrix
read_nd_vector(xsdata_grp, "chi-delayed", chi_delayed);
// Normalize the chi info so the CDF is 1.
for (int a = 0; a < n_ang; a++) {
for (int dg = 0; dg < delayed_groups; dg++) {
for (int gin = 0; gin < energy_groups; gin++) {
double chi_sum = 0.;
for (int gout = 0; gout < energy_groups; gout++) {
chi_sum += chi_delayed[a][gin][gout][dg];
}
if (chi_sum > 0.) {
for (int gout = 0; gout < energy_groups; gout++) {
chi_delayed[a][gin][gout][dg] /= chi_sum;
}
} else {
fatal_error("Encountered chi-delayed for a group that is <= 0!");
}
}
}
}
} else {
fatal_error("chi-delayed must be provided as a 3D or 4D array!");
}
}
// Get prompt-nu-fission, if present
if (object_exists(xsdata_grp, "prompt-nu-fission")) {
hid_t xsdata = open_dataset(xsdata_grp, "prompt-nu-fission");
int ndims = dataset_ndims(xsdata);
// raise ndims to make the isotropic ndims the same as angular
if (is_isotropic) ndims += 2;
close_dataset(xsdata);
if (ndims == 3) {
// prompt-nu-fission is a [in group] vector
read_nd_vector(xsdata_grp, "prompt-nu-fission",
prompt_nu_fission);
} else if (ndims == 4) {
// prompt nu fission is a matrix,
// so set prompt_nu_fiss & chi_prompt
double_3dvec temp_arr(n_ang, double_2dvec(energy_groups,
double_1dvec(energy_groups)));
read_nd_vector(xsdata_grp, "prompt-nu-fission", temp_arr);
// The prompt_nu_fission vector from the matrix form
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
double prompt_sum = std::accumulate(temp_arr[a][gin].begin(),
temp_arr[a][gin].end(), 0.);
prompt_nu_fission[a][gin] = prompt_sum;
}
// The chi_prompt data is just the normalized fission matrix
for (int gin= 0; gin < energy_groups; gin++) {
if (prompt_nu_fission[a][gin] > 0.) {
for (int gout = 0; gout < energy_groups; gout++) {
chi_prompt[a][gin][gout] =
temp_arr[a][gin][gout] / prompt_nu_fission[a][gin];
}
} else {
fatal_error("Encountered chi-prompt for a group that is <= 0!");
}
}
}
} else {
fatal_error("prompt-nu-fission must be provided as a 3D or 4D array!");
}
}
// Get delayed-nu-fission, if present
if (object_exists(xsdata_grp, "delayed-nu-fission")) {
hid_t xsdata = open_dataset(xsdata_grp, "delayed-nu-fission");
int ndims = dataset_ndims(xsdata);
close_dataset(xsdata);
// raise ndims to make the isotropic ndims the same as angular
if (is_isotropic) ndims += 2;
if (ndims == 3) {
// delayed-nu-fission is an [in group] vector
if (temp_beta[0][0][0] == 0.) {
fatal_error("cannot set delayed-nu-fission with a 1D array if "
"beta is not provided");
}
double_2dvec temp_arr(n_ang, double_1dvec(energy_groups));
read_nd_vector(xsdata_grp, "delayed-nu-fission", temp_arr);
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
for (int dg = 0; dg < delayed_groups; dg++) {
// Set delayed-nu-fission using beta
delayed_nu_fission[a][gin][dg] =
temp_beta[a][gin][dg] * temp_arr[a][gin];
}
}
}
} else if (ndims == 4) {
read_nd_vector(xsdata_grp, "delayed-nu-fission",
delayed_nu_fission);
} else if (ndims == 5) {
// This will contain delayed-nu-fision and chi-delayed data
double_4dvec temp_arr(n_ang, double_3dvec(energy_groups,
double_2dvec(energy_groups, double_1dvec(delayed_groups))));
read_nd_vector(xsdata_grp, "delayed-nu-fission", temp_arr);
// Set the 3D delayed-nu-fission matrix and 4D chi-delayed matrix
// from the 4D delayed-nu-fission matrix
for (int a = 0; a < n_ang; a++) {
for (int dg = 0; dg < delayed_groups; dg++) {
for (int gin = 0; gin < energy_groups; gin++) {
double gout_sum = 0.;
for (int gout = 0; gout < energy_groups; gout++) {
gout_sum += temp_arr[a][gin][gout][dg];
chi_delayed[a][gin][gout][dg] = temp_arr[a][gin][gout][dg];
}
delayed_nu_fission[a][gin][dg] = gout_sum;
// Normalize chi-delayed
if (gout_sum > 0.) {
for (int gout = 0; gout < energy_groups; gout++) {
chi_delayed[a][gin][gout][dg] /= gout_sum;
}
} else {
fatal_error("Encountered chi-delayed for a group that is <= 0!");
}
}
}
}
} else {
fatal_error("prompt-nu-fission must be provided as a 3D, 4D, or 5D "
"array!");
}
}
// Combine prompt_nu_fission and delayed_nu_fission into nu_fission
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
nu_fission[a][gin] =
std::accumulate(delayed_nu_fission[a][gin].begin(),
delayed_nu_fission[a][gin].end(),
prompt_nu_fission[a][gin]);
}
}
}
//==============================================================================
void
XsData::scatter_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
int energy_groups, int scatter_format, int final_scatter_format,
int order_data, int max_order, int legendre_to_tabular_points)
{
int n_ang = n_pol * n_azi;
if (!object_exists(xsdata_grp, "scatter_data")) {
fatal_error("Must provide scatter_data group!");
}
hid_t scatt_grp = open_group(xsdata_grp, "scatter_data");
// Get the outgoing group boundary indices
int_2dvec gmin(n_ang, int_1dvec(energy_groups));
read_nd_vector(scatt_grp, "g_min", gmin, true);
int_2dvec gmax(n_ang, int_1dvec(energy_groups));
read_nd_vector(scatt_grp, "g_max", gmax, true);
// Make gmin and gmax start from 0 vice 1 as they do in the library
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
gmin[a][gin] -= 1;
gmax[a][gin] -= 1;
}
}
// Now use this info to find the length of a vector to hold the flattened
// data.
int length = 0;
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
length += order_data * (gmax[a][gin] - gmin[a][gin] + 1);
}
}
double_1dvec temp_arr(length);
read_nd_vector(scatt_grp, "scatter_matrix", temp_arr, true);
// Compare the number of orders given with the max order of the problem;
// strip off the superfluous orders if needed
int order_dim;
if (scatter_format == ANGLE_LEGENDRE) {
order_dim = std::min(order_data - 1, max_order) + 1;
} else {
order_dim = order_data;
}
// convert the flattened temp_arr to a jagged array for passing to
// scatt data
double_4dvec input_scatt(n_ang, double_3dvec(energy_groups));
int temp_idx = 0;
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
input_scatt[a][gin].resize(gmax[a][gin] - gmin[a][gin] + 1);
for (int i_gout = 0; i_gout < input_scatt[a][gin].size(); i_gout++) {
input_scatt[a][gin][i_gout].resize(order_dim);
for (int l = 0; l < order_dim; l++) {
input_scatt[a][gin][i_gout][l] = temp_arr[temp_idx++];
}
// Adjust index for the orders we didnt take
temp_idx += (order_data - order_dim);
}
}
}
temp_arr.clear();
// Get multiplication matrix
double_3dvec temp_mult(n_ang, double_2dvec(energy_groups));
if (object_exists(scatt_grp, "multiplicity_matrix")) {
temp_arr.resize(length / order_data);
read_nd_vector(scatt_grp, "multiplicity_matrix", temp_arr);
// convert the flat temp_arr to a jagged array for passing to scatt data
int temp_idx = 0;
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
temp_mult[a][gin].resize(gmax[a][gin] - gmin[a][gin] + 1);
for (int i_gout = 0; i_gout < temp_mult[a][gin].size(); i_gout++) {
temp_mult[a][gin][i_gout] = temp_arr[temp_idx++];
}
}
}
} else {
// Use a default: multiplicities are 1.0.
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
temp_mult[a][gin].resize(gmax[a][gin] - gmin[a][gin] + 1);
for (int i_gout = 0; i_gout < temp_mult[a][gin].size(); i_gout++) {
temp_mult[a][gin][i_gout] = 1.;
}
}
}
}
temp_arr.clear();
close_group(scatt_grp);
// Finally, convert the Legendre data to tabular, if needed
if (scatter_format == ANGLE_LEGENDRE &&
final_scatter_format == ANGLE_TABULAR) {
for (int a = 0; a < n_ang; a++) {
ScattDataLegendre legendre_scatt;
legendre_scatt.init(gmin[a], gmax[a], temp_mult[a], input_scatt[a]);
// Now create a tabular version of legendre_scatt
convert_legendre_to_tabular(legendre_scatt,
*static_cast<ScattDataTabular*>(scatter[a].get()),
legendre_to_tabular_points);
scatter_format = final_scatter_format;
}
} else {
// We are sticking with the current representation
// Initialize the ScattData object with this data
for (int a = 0; a < n_ang; a++) {
scatter[a]->init(gmin[a], gmax[a], temp_mult[a], input_scatt[a]);
}
}
}
//==============================================================================
void
XsData::combine(const std::vector<XsData*>& those_xs,
const double_1dvec& scalars)
{
// Combine the non-scattering data
for (int i = 0; i < those_xs.size(); i++) {
XsData* that = those_xs[i];
if (!equiv(*that)) fatal_error("Cannot combine the XsData objects!");
double scalar = scalars[i];
for (int a = 0; a < total.size(); a++) {
for (int gin = 0; gin < total[a].size(); gin++) {
total[a][gin] += scalar * that->total[a][gin];
absorption[a][gin] += scalar * that->absorption[a][gin];
if (i == 0) {
inverse_velocity[a][gin] = that->inverse_velocity[a][gin];
}
if (that->prompt_nu_fission.size() > 0) {
nu_fission[a][gin] += scalar * that->nu_fission[a][gin];
prompt_nu_fission[a][gin] +=
scalar * that->prompt_nu_fission[a][gin];
kappa_fission[a][gin] += scalar * that->kappa_fission[a][gin];
fission[a][gin] += scalar * that->fission[a][gin];
for (int dg = 0; dg < delayed_nu_fission[a][gin].size(); dg++) {
delayed_nu_fission[a][gin][dg] +=
scalar * that->delayed_nu_fission[a][gin][dg];
}
for (int gout = 0; gout < chi_prompt[a][gin].size(); gout++) {
chi_prompt[a][gin][gout] +=
scalar * that->chi_prompt[a][gin][gout];
for (int dg = 0; dg < chi_delayed[a][gin][gout].size(); dg++) {
chi_delayed[a][gin][gout][dg] +=
scalar * that->chi_delayed[a][gin][gout][dg];
}
}
}
}
for (int dg = 0; dg < decay_rate[a].size(); dg++) {
decay_rate[a][dg] += scalar * that->decay_rate[a][dg];
}
// Normalize chi
if (chi_prompt.size() > 0) {
for (int gin = 0; gin < chi_prompt[a].size(); gin++) {
double norm = std::accumulate(chi_prompt[a][gin].begin(),
chi_prompt[a][gin].end(), 0.);
if (norm > 0.) {
for (int gout = 0; gout < chi_prompt[a][gin].size(); gout++) {
chi_prompt[a][gin][gout] /= norm;
}
}
for (int dg = 0; dg < chi_delayed[a][gin][0].size(); dg++) {
norm = 0.;
for (int gout = 0; gout < chi_delayed[a][gin].size(); gout++) {
norm += chi_delayed[a][gin][gout][dg];
}
if (norm > 0.) {
for (int gout = 0; gout < chi_delayed[a][gin].size(); gout++) {
chi_delayed[a][gin][gout][dg] /= norm;
}
}
}
}
}
}
}
// Allow the ScattData object to combine itself
for (int a = 0; a < total.size(); a++) {
// Build vector of the scattering objects to incorporate
std::vector<ScattData*> those_scatts(those_xs.size());
for (int i = 0; i < those_xs.size(); i++) {
those_scatts[i] = those_xs[i]->scatter[a].get();
}
// Now combine these guys
scatter[a]->combine(those_scatts, scalars);
}
}
//==============================================================================
bool
XsData::equiv(const XsData& that)
{
return ((absorption.size() == that.absorption.size()) &&
(absorption[0].size() == that.absorption[0].size()));
}
} //namespace openmc

118
src/xsdata.h Normal file
View file

@ -0,0 +1,118 @@
//! \file xsdata.h
//! A collection of classes for containing the Multi-Group Cross Section data
#ifndef XSDATA_H
#define XSDATA_H
#include <memory>
#include <vector>
#include "hdf5_interface.h"
#include "scattdata.h"
namespace openmc {
//==============================================================================
// XSDATA contains the temperature-independent cross section data for an MGXS
//==============================================================================
class XsData {
private:
//! \brief Reads scattering data from the HDF5 file
void
scatter_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, int energy_groups,
int scatter_format, int final_scatter_format, int order_data,
int max_order, int legendre_to_tabular_points);
//! \brief Reads fission data from the HDF5 file
void
fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, int energy_groups,
int delayed_groups, bool is_isotropic);
public:
// The following quantities have the following dimensions:
// [angle][incoming group]
double_2dvec total;
double_2dvec absorption;
double_2dvec nu_fission;
double_2dvec prompt_nu_fission;
double_2dvec kappa_fission;
double_2dvec fission;
double_2dvec inverse_velocity;
// decay_rate has the following dimensions:
// [angle][delayed group]
double_2dvec decay_rate;
// delayed_nu_fission has the following dimensions:
// [angle][incoming group][delayed group]
double_3dvec delayed_nu_fission;
// chi_prompt has the following dimensions:
// [angle][incoming group][outgoing group]
double_3dvec chi_prompt;
// chi_delayed has the following dimensions:
// [angle][incoming group][outgoing group][delayed group]
double_4dvec chi_delayed;
// scatter has the following dimensions: [angle]
std::vector<std::shared_ptr<ScattData> > scatter;
XsData() = default;
//! \brief Constructs the XsData object metadata.
//!
//! @param num_groups Number of energy groups.
//! @param num_delayed_groups Number of delayed groups.
//! @param fissionable Is this a fissionable data set or not.
//! @param scatter_format The scattering representation of the file.
//! @param n_pol Number of polar angles.
//! @param n_azi Number of azimuthal angles.
XsData(int num_groups, int num_delayed_groups, bool fissionable,
int scatter_format, int n_pol, int n_azi);
//! \brief Loads the XsData object from the HDF5 file
//!
//! @param xs_id HDF5 group id for the cross section data.
//! @param fissionable Is this a fissionable data set or not.
//! @param scatter_format The scattering representation of the file.
//! @param final_scatter_format The scattering representation after reading;
//! this is different from scatter_format if converting a Legendre to
//! a tabular representation.
//! @param order_data The dimensionality of the scattering data in the file.
//! @param max_order Maximum order requested by the user;
//! this is only used for Legendre scattering.
//! @param legendre_to_tabular Flag to denote if any Legendre provided
//! should be converted to a Tabular representation.
//! @param legendre_to_tabular_points If a conversion is requested, this
//! provides the number of points to use in the tabular representation.
//! @param is_isotropic Is this an isotropic or angular with respect to
//! the incoming particle.
//! @param n_pol Number of polar angles.
//! @param n_azi Number of azimuthal angles.
void
from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format,
int final_scatter_format, int order_data, int max_order,
int legendre_to_tabular_points, bool is_isotropic, int n_pol,
int n_azi);
//! \brief Combines the microscopic data to a macroscopic object.
//!
//! @param micros Microscopic objects to combine.
//! @param scalars Scalars to multiply the microscopic data by.
void
combine(const std::vector<XsData*>& those_xs, const double_1dvec& scalars);
//! \brief Checks to see if this and that are able to be combined
//!
//! This comparison is used when building macroscopic cross sections
//! from microscopic cross sections.
//! @param that The other XsData to compare to this one.
//! @return True if they can be combined.
bool
equiv(const XsData& that);
};
} //namespace openmc
#endif // XSDATA_H

129
tests/unit_tests/test_data_neutron.py
View file

@ -37,9 +37,10 @@ def sm150():
@pytest.fixture(scope='module')
def gd154():
"""Gd154 ENDF data (contains Reich Moore resonance range)"""
"""Gd154 ENDF data (contains Reich Moore resonance range and reosnance
covariance with LCOMP=1)."""
filename = os.path.join(_ENDF_DATA, 'neutrons', 'n-064_Gd_154.endf')
return openmc.data.IncidentNeutron.from_endf(filename)
return openmc.data.IncidentNeutron.from_endf(filename, covariance=True)
@pytest.fixture(scope='module')
@ -77,6 +78,13 @@ def na22():
return openmc.data.IncidentNeutron.from_endf(filename)
@pytest.fixture(scope='module')
def na23():
"""Na23 ENDF data (contains MLBW resonance covariance with LCOMP=0)."""
filename = os.path.join(_ENDF_DATA, 'neutrons', 'n-011_Na_023.endf')
return openmc.data.IncidentNeutron.from_endf(filename, covariance=True)
@pytest.fixture(scope='module')
def be9():
"""Be9 ENDF data (contains laboratory angle-energy distribution)."""
@ -97,6 +105,28 @@ def am244():
return openmc.data.IncidentNeutron.from_njoy(endf_file)
@pytest.fixture(scope='module')
def ti50():
"""Ti50 ENDF data (contains Multi-level Breit-Wigner resonance range and
resonance covariance with LCOMP=1)."""
filename = os.path.join(_ENDF_DATA, 'neutrons', 'n-022_Ti_050.endf')
return openmc.data.IncidentNeutron.from_endf(filename, covariance=True)
@pytest.fixture(scope='module')
def cf252():
"""Cf252 ENDF data (contains RM resonance covariance with LCOMP=0)."""
filename = os.path.join(_ENDF_DATA, 'neutrons', 'n-098_Cf_252.endf')
return openmc.data.IncidentNeutron.from_endf(filename, covariance=True)
@pytest.fixture(scope='module')
def th232():
"""Th232 ENDF data (contains RM resonance covariance with LCOMP=2)."""
filename = os.path.join(_ENDF_DATA, 'neutrons', 'n-090_Th_232.endf')
return openmc.data.IncidentNeutron.from_endf(filename, covariance=True)
def test_attributes(pu239):
assert pu239.name == 'Pu239'
assert pu239.mass_number == 239
@ -277,6 +307,101 @@ def test_rml(cl35):
assert isinstance(group, openmc.data.SpinGroup)
def test_mlbw_cov_lcomp0(cf252):
# Testing on first range only
cov = cf252.resonance_covariance.ranges[0]
res = cf252.resonances.ranges[0]
assert cov.parameters['energy'][0] == pytest.approx(-3.5)
assert res.parameters['energy'][0] == cov.parameters['energy'][0]
assert isinstance(cov, openmc.data.resonance_covariance.MultiLevelBreitWignerCovariance)
assert cov.energy_min == pytest.approx(1e-5)
assert cov.energy_max == pytest.approx(1000.)
assert cov.covariance[0,0] == pytest.approx(1.225e-05)
subset = cov.subset('energy', [0, 100])
assert not subset.parameters.empty
assert (subset.file2res.parameters['energy'] < 100).all()
samples = cov.sample(1)
xs = samples[0].reconstruct([10., 100., 1000.])
assert sorted(xs.keys()) == [2, 18, 102]
def test_mlbw_cov_lcomp1(ti50):
# Testing on first range only
cov = ti50.resonance_covariance.ranges[0]
res = ti50.resonances.ranges[0]
assert cov.parameters['energy'][0] == pytest.approx(-21020.)
assert res.parameters['energy'][0] == cov.parameters['energy'][0]
assert isinstance(cov, openmc.data.resonance_covariance.MultiLevelBreitWignerCovariance)
assert cov.energy_min == pytest.approx(1e-5)
assert cov.energy_max == pytest.approx(587000.)
assert cov.covariance[0,0] == pytest.approx(1.410177e5)
subset = cov.subset('L', [1, 1])
assert not subset.parameters.empty
assert (subset.file2res.parameters['L'] == 1).all()
samples = cov.sample(1)
xs = samples[0].reconstruct([10., 100., 1000.])
assert sorted(xs.keys()) == [2, 18, 102]
def test_mlbw_cov_lcomp2(na23):
# Testing on first range only
cov = na23.resonance_covariance.ranges[0]
res = na23.resonances.ranges[0]
assert cov.parameters['energy'][0] == pytest.approx(2810.)
assert res.parameters['energy'][0] == cov.parameters['energy'][0]
assert isinstance(cov, openmc.data.resonance_covariance.MultiLevelBreitWignerCovariance)
assert cov.energy_min == pytest.approx(600)
assert cov.energy_max == pytest.approx(500000.)
assert cov.covariance[0,0] == pytest.approx(16.1064163584)
subset = cov.subset('L', [1, 1])
assert not subset.parameters.empty
assert (subset.file2res.parameters['L'] == 1).all()
samples = cov.sample(1)
xs = samples[0].reconstruct([10., 100., 1000.])
assert sorted(xs.keys()) == [2, 18, 102]
def test_rmcov_lcomp1(gd154):
# Testing on first range only
cov = gd154.resonance_covariance.ranges[0]
res = gd154.resonances.ranges[0]
assert cov.parameters['energy'][0] == pytest.approx(-2.200001)
assert res.parameters['energy'][0] == cov.parameters['energy'][0]
assert isinstance(cov, openmc.data.resonance_covariance.ReichMooreCovariance)
assert cov.energy_min == pytest.approx(1e-5)
assert cov.energy_max == pytest.approx(2760.)
assert cov.covariance[0,0] == pytest.approx(0.8895997)
subset = cov.subset('energy', [0, 100])
assert not subset.parameters.empty
assert (subset.file2res.parameters['energy'] < 100).all()
samples = cov.sample(1)
xs = samples[0].reconstruct([10., 100., 1000.])
assert sorted(xs.keys()) == [2, 18, 102]
def test_rmcov_lcomp2(th232):
# Testing on first range only
cov = th232.resonance_covariance.ranges[0]
res = th232.resonances.ranges[0]
assert cov.parameters['energy'][0] == pytest.approx(-2000)
assert res.parameters['energy'][0] == cov.parameters['energy'][0]
assert isinstance(cov, openmc.data.resonance_covariance.ReichMooreCovariance)
assert cov.energy_min == pytest.approx(1e-5)
assert cov.energy_max == pytest.approx(4000.)
assert cov.covariance[0,0] == pytest.approx(246.6043092496)
subset = cov.subset('energy', [0, 100])
assert not subset.parameters.empty
assert (subset.file2res.parameters['energy'] < 100).all()
samples = cov.sample(1)
xs = samples[0].reconstruct([10., 100., 1000.])
assert sorted(xs.keys()) == [2, 18, 102]
def test_madland_nix(am241):
fission = am241.reactions[18]
prompt_neutron = fission.products[0]

3
tools/ci/travis-install.sh
View file

@ -4,6 +4,9 @@ set -ex
# Install NJOY 2016
./tools/ci/travis-install-njoy.sh
# Upgrade pip before doing anything else
pip install --upgrade pip
# Running OpenMC's setup.py requires numpy/cython already. NumPy float
# formatting changed in version 1.14, so stick with a lower version until we can
# handle it in our test suite

0
src/faddeeva/Faddeeva.c → vendor/faddeeva/Faddeeva.c vendored
View file

0
src/faddeeva/Faddeeva.cc → vendor/faddeeva/Faddeeva.cc vendored
View file

0
src/faddeeva/Faddeeva.h → vendor/faddeeva/Faddeeva.h vendored
View file

0
src/pugixml/pugiconfig.hpp → vendor/pugixml/pugiconfig.hpp vendored
View file

0
src/pugixml/pugixml.cpp → vendor/pugixml/pugixml.cpp vendored
View file

0
src/pugixml/pugixml.hpp → vendor/pugixml/pugixml.hpp vendored
View file