crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-29 06:35:48 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge branch 'develop' into cmfd-helper2

Browse source
This commit is contained in:
Shikhar Kumar 2018-10-07 14:58:17 -04:00
commit d3b994b318
55 changed files with 1656 additions and 192 deletions
Download patch file Download diff file Expand all files Collapse all files

5
.travis.yml
View file

@ -13,6 +13,8 @@ addons:
- libmpich-dev
- libhdf5-serial-dev
- libhdf5-mpich-dev
- libblas-dev
- liblapack-dev
cache:
directories:
- $HOME/nndc_hdf5
@ -27,6 +29,7 @@ env:
- OPENMC_CROSS_SECTIONS=$HOME/nndc_hdf5/cross_sections.xml
- OPENMC_ENDF_DATA=$HOME/endf-b-vii.1
- OPENMC_MULTIPOLE_LIBRARY=$HOME/WMP_Library
- LD_LIBRARY_PATH=$HOME/MOAB/lib:$HOME/DAGMC/lib
- PATH=$PATH:$HOME/NJOY2016/build
- DISPLAY=:99.0
- COVERALLS_PARALLEL=true
@ -35,6 +38,8 @@ env:
- OMP=y MPI=n PHDF5=n
- OMP=n MPI=y PHDF5=n
- OMP=n MPI=y PHDF5=y
- OMP=n MPI=y PHDF5=y DAGMC=y
- OMP=y MPI=y PHDF5=y DAGMC=y
notifications:
webhooks: https://coveralls.io/webhook?repo_token=$COVERALLS_REPO_TOKEN
install:

26
CMakeLists.txt
View file

@ -19,6 +19,9 @@ 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)
option(dagmc "Enable support for DAGMC (CAD) geometry" OFF)
# Maximum number of nested coordinates levels
set(maxcoord 10 CACHE STRING "Maximum number of nested coordinate levels")
#===============================================================================
@ -36,12 +39,23 @@ if(MPI_ENABLED AND mpif08)
message("-- Using Fortran 2008 MPI bindings")
endif()
#===============================================================================
# DAGMC Geometry Support - need DAGMC/MOAB
#===============================================================================
if(dagmc)
find_package(DAGMC REQUIRED)
if(NOT DAGMC_FOUND)
message(FATAL_ERROR "Could not find DAGMC installation")
endif()
link_directories(${DAGMC_LIBRARY_DIRS})
endif()
#===============================================================================
# HDF5 for binary output
#===============================================================================
# Allow user to specify HDF5_ROOT
if (CMAKE_VERSION VERSION_GREATER_EQUAL 3.12)
if (NOT (CMAKE_VERSION VERSION_LESS 3.12))
cmake_policy(SET CMP0074 NEW)
endif()
@ -295,6 +309,7 @@ add_library(libopenmc SHARED
src/algorithm.F90
src/bank_header.F90
src/api.F90
src/dagmc_header.F90
src/cmfd_data.F90
src/cmfd_execute.F90
src/cmfd_header.F90
@ -383,6 +398,7 @@ add_library(libopenmc SHARED
src/tallies/tally_header.F90
src/tallies/trigger.F90
src/tallies/trigger_header.F90
src/dagmc.cpp
src/cell.cpp
src/cmfd_execute.cpp
src/distribution.cpp
@ -428,6 +444,7 @@ add_library(libopenmc SHARED
src/thermal.cpp
src/xml_interface.cpp
src/xsdata.cpp)
set_target_properties(libopenmc PROPERTIES
OUTPUT_NAME openmc
LINKER_LANGUAGE Fortran)
@ -473,10 +490,15 @@ endif()
target_link_libraries(libopenmc ${ldflags} ${HDF5_LIBRARIES} pugixml
faddeeva xtensor)
if(dagmc)
target_compile_definitions(libopenmc PRIVATE DAGMC)
target_link_libraries(libopenmc ${DAGMC_LIBRARIES})
target_include_directories(libopenmc PRIVATE ${DAGMC_INCLUDE_DIRS})
endif()
#===============================================================================
# openmc executable
#===============================================================================
add_executable(openmc src/main.cpp)
target_compile_options(openmc PRIVATE ${cxxflags})
target_link_libraries(openmc libopenmc)

39
Dockerfile Normal file
View file

@ -0,0 +1,39 @@
FROM ubuntu:latest
# Setup environment variables for Docker image
ENV FC=/usr/bin/mpif90 CC=/usr/bin/mpicc CXX=/usr/bin/mpicxx \
PATH=/opt/openmc/bin:/opt/NJOY2016/build:$PATH \
LD_LIBRARY_PATH=/opt/openmc/lib:$LD_LIBRARY_PATH \
OPENMC_CROSS_SECTIONS=/root/nndc_hdf5/cross_sections.xml \
OPENMC_MULTIPOLE_LIBRARY=/root/WMP_Library \
OPENMC_ENDF_DATA=/root/endf-b-vii.1
# Install dependencies from Debian package manager
RUN apt-get update -y && \
apt-get upgrade -y && \
apt-get install -y python3-pip && \
apt-get install -y wget git emacs && \
apt-get install -y gfortran g++ cmake && \
apt-get install -y mpich libmpich-dev && \
apt-get install -y libhdf5-serial-dev libhdf5-mpich-dev && \
apt-get install -y imagemagick && \
apt-get autoremove
# Update system-provided pip
RUN pip3 install --upgrade pip
# Clone and install NJOY2016
RUN git clone https://github.com/njoy/NJOY2016 /opt/NJOY2016 && \
cd /opt/NJOY2016 && \
mkdir build && cd build && \
cmake -Dstatic=on .. && make 2>/dev/null && make install
# Clone and install OpenMC
RUN git clone https://github.com/openmc-dev/openmc.git /opt/openmc && \
cd /opt/openmc && mkdir -p build && cd build && \
cmake -Doptimize=on -DHDF5_PREFER_PARALLEL=on .. && \
make && make install && \
cd .. && pip install -e .[test]
# Download cross sections (NNDC and WMP) and ENDF data needed by test suite
RUN ./opt/openmc/tools/ci/download-xs.sh

18
cmake/Modules/FindDAGMC.cmake Normal file
View file

@ -0,0 +1,18 @@
# Try to find DAGMC
#
# Once done this will define
#
# DAGMC_FOUND - system has DAGMC
# DAGMC_INCLUDE_DIRS - the DAGMC include directory
# DAGMC_LIBRARIES - Link these to use DAGMC
# DAGMC_DEFINITIONS - Compiler switches required for using DAGMC
find_path(DAGMC_CMAKE_CONFIG NAMES DAGMCConfig.cmake
HINTS ${DAGMC_ROOT} $ENV{DAGMC_ROOT}
PATHS ENV LD_LIBRARY_PATH
PATH_SUFFIXES lib Lib cmake lib/cmake
NO_DEFAULT_PATH)
message(STATUS "Found DAGMC in ${DAGMC_CMAKE_CONFIG}")
include(${DAGMC_CMAKE_CONFIG}/DAGMCConfig.cmake)

57
docs/source/devguide/docker.rst Normal file
View file

@ -0,0 +1,57 @@
.. _devguide_docker:
======================
Deployment with Docker
======================
OpenMC can be easily deployed using `Docker <https://www.docker.com/>`_ on any
Windows, Mac or Linux system. With Docker running, execute the following
command in the shell to build a `Docker image`_ called ``debian/openmc:latest``:
.. code-block:: sh
docker build -t debian/openmc:latest https://github.com/openmc-dev/openmc.git#develop
.. note:: This may take 5 -- 10 minutes to run to completion.
This command will execute the instructions in OpenMC's ``Dockerfile`` to
build a Docker image with OpenMC installed. The image includes OpenMC with
MPICH and parallel HDF5 in the ``/opt/openmc`` directory, and
`Miniconda3 <https://conda.io/miniconda.html>`_ with all of the Python
pre-requisites (NumPy, SciPy, Pandas, etc.) installed. The
`NJOY2016 <http://www.njoy21.io/NJOY2016/>`_ codebase is installed in
``/opt/NJOY2016`` to support full functionality and testing of the
``openmc.data`` Python module. The publicly available nuclear data libraries
necessary to run OpenMC's test suite -- including NNDC and WMP cross sections
and ENDF data -- are in the ``/opt/openmc/data directory``, and the
corresponding :envvar:`OPENMC_CROSS_SECTIONS`,
:envvar:`OPENMC_MULTIPOLE_LIBRARY`, and :envvar:`OPENMC_ENDF_DATA`
environment variables are initialized.
After building the Docker image, you can run the following to see the names of
all images on your machine, including ``debian/openmc:latest``:
.. code-block:: sh
docker image ls
Now you can run the following to create a `Docker container`_ called
``my_openmc`` based on the ``debian/openmc:latest`` image:
.. code-block:: sh
docker run -it --name=my_openmc debian/openmc:latest
This command will open an interactive shell running from within the
Docker container where you have access to use OpenMC.
.. note:: The ``docker run`` command supports many
`options <https://docs.docker.com/engine/reference/commandline/run/>`_
for spawning containers -- including `mounting volumes`_ from the
host filesystem -- which many users will find useful.
.. _Docker image: https://docs.docker.com/engine/reference/commandline/images/
.. _Docker container: https://www.docker.com/resources/what-container
.. _options: https://docs.docker.com/engine/reference/commandline/run/
.. _mounting volumes: https://docs.docker.com/storage/volumes/

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

@ -18,3 +18,4 @@ other related topics.
tests
user-input
docbuild
docker

53
docs/source/io_formats/fission_energy.rst
View file

@ -1,53 +0,0 @@
.. _usersguide_fission_energy:
==================================
Fission Energy Release File Format
==================================
This file is a compact HDF5 representation of the ENDF MT=1, MF=458 data (see
ENDF-102_ for details). It gives the information needed to compute the energy
carried away from fission reactions by each reaction product (e.g. fragment
nuclei, neutrons) which depends on the incident neutron energy. OpenMC is
distributed with one of these files under
data/fission_Q_data_endfb71.h5. More files of this format can be created from
ENDF files with the
``openmc.data.write_compact_458_library`` function. They can be read with the
``openmc.data.FissionEnergyRelease.from_compact_hdf5`` class method.
:Attributes: - **comment** (*char[]*) -- An optional text comment
- **component order** (*char[][]*) -- An array of strings
specifying the order each reaction product occurs in the data
arrays. The components use the 2-3 letter abbreviations
specified in ENDF-102 e.g. EFR for fission fragments and ENP for
prompt neutrons.
**/<nuclide name>/**
Nuclides are named by concatenating their atomic symbol and mass number. For
example, 'U235' or 'Pu239'. Metastable nuclides are appended with an
'_m' and their metastable number. For example, 'Am242_m1'
:Datasets:
- **data** (*double[][][]*) -- The energy release coefficients. The
first axis indexes the component type. The second axis specifies
values or uncertainties. The third axis indexes the polynomial
order. If the data uses the Sher-Beck format, then the last axis
will have a length of one and ENDF-102 should be consulted for
energy dependence. Otherwise, the data uses the Madland format
which is a polynomial of incident energy.
For example, if 'EFR' is given first in the **component order**
attribute and the data uses the Madland format, then the energy
released in the form of fission fragments at an incident energy
:math:`E` is given by
.. math::
\text{data}[0, 0, 0] + \text{data}[0, 0, 1] \cdot E
+ \text{data}[0, 0, 2] \cdot E^2 + \ldots
And its uncertainty is
.. math::
\text{data}[0, 1, 0] + \text{data}[0, 1, 1] \cdot E
+ \text{data}[0, 1, 2] \cdot E^2 + \ldots
.. _ENDF-102: http://www.nndc.bnl.gov/endfdocs/ENDF-102-2012.pdf

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

@ -34,7 +34,6 @@ Data Files
nuclear_data
mgxs_library
data_wmp
fission_energy
------------
Output Files

2
docs/source/io_formats/summary.rst
View file

@ -24,6 +24,8 @@ The current version of the summary file format is 6.0.
- **n_universes** (*int*) -- Number of unique universes in the
problem.
- **n_lattices** (*int*) -- Number of lattices in the problem.
- **dagmc** (*int*) -- Indicates that a DAGMC geometry was used
if present.
**/geometry/cells/cell <uid>/**

693
docs/source/methods/photon_physics.rst
View file

@ -85,26 +85,26 @@ momentum transfer is traditionally expressed as
.. math::
:label: momentum-transfer
x = \kappa \alpha \sqrt{1 - \mu}
x = a k \sqrt{1 - \mu}
where :math:`\alpha` is the ratio of the photon energy to the electron rest
mass, and the coefficient :math:`\kappa` can be shown to be
where :math:`k` is the ratio of the photon energy to the electron rest
mass, and the coefficient :math:`a` can be shown to be
.. math::
:label: kappa
:label: omega
\kappa = \frac{m_e c^2}{\sqrt{2}hc} \approx 29.14329,
a = \frac{m_e c^2}{\sqrt{2}hc} \approx 29.14329~\unicode{x212B},
where :math:`m_e` is the mass of the electron, :math:`c` is the speed of light
in a vacuum, and :math:`h` is Planck's constant. Using :eq:`momentum-transfer`,
we have :math:`\mu = 1 - [x/(\kappa\alpha)]^2` and :math:`d\mu/dx^2 =
-1/(\kappa\alpha)^2`. The probability density in :math:`x^2` is
we have :math:`\mu = 1 - [x/(ak)]^2` and :math:`d\mu/dx^2 =
-1/(ak)^2`. The probability density in :math:`x^2` is
.. math::
:label: coherent-pdf-x2
p(x^2) dx^2 = p(\mu) \left | \frac{d\mu}{dx^2} \right | dx^2 = \frac{2\pi
r_e^2 A(\bar{x}^2,Z)}{(\kappa\alpha)^2 \sigma(E)} \left (
r_e^2 A(\bar{x}^2,Z)}{(ak)^2 \sigma(E)} \left (
\frac{1 + \mu^2}{2} \right ) \left ( \frac{F(x, Z)^2}{A(\bar{x}^2, Z)} \right ) dx^2
where :math:`\bar{x}` is the maximum value of :math:`x` that occurs for
@ -113,7 +113,7 @@ where :math:`\bar{x}` is the maximum value of :math:`x` that occurs for
.. math::
:label: xmax
\bar{x} = \kappa \alpha \sqrt{2} = \frac{m_e c^2}{hc} \alpha,
\bar{x} = a k \sqrt{2} = \frac{m_e c^2}{hc} k,
and :math:`A(x^2, Z)` is the integral of the square of the form factor:
@ -154,6 +154,8 @@ section. The complete algorithm is as follows:
6. If :math:`\xi_2 < (1 + \mu^2)/2`, accept :math:`\mu`. Otherwise, repeat the
sampling at step 3.
.. _incoherent-sampling:
Incoherent (Compton) Scattering
-------------------------------
@ -166,12 +168,11 @@ the two authors who discovered it:
.. math::
:label: klein-nishina
\frac{d\sigma_{KN}}{d\mu} = \pi r_e^2 \left ( \frac{\alpha'}{\alpha} \right
)^2 \left [ \frac{\alpha'}{\alpha} + \frac{\alpha}{\alpha'} + \mu^2 - 1
\right ]
\frac{d\sigma_{KN}}{d\mu} = \pi r_e^2 \left ( \frac{k'}{k} \right)^2 \left
[ \frac{k'}{k} + \frac{k}{k'} + \mu^2 - 1 \right ]
where :math:`\alpha` and :math:`\alpha'` are the ratios of the incoming and
exiting photon energies to the electron rest mass energy equivalent (0.511 MeV),
where :math:`k` and :math:`k'` are the ratios of the incoming and exiting
photon energies to the electron rest mass energy equivalent (0.511 MeV),
respectively. Although it appears that the outgoing energy and angle are
separate, there is actually a one-to-one relationship between them such that
only one needs to be sampled:
@ -179,32 +180,31 @@ only one needs to be sampled:
.. math::
:label: compton-energy-angle
\alpha' = \frac{\alpha}{1 + \alpha(1 - \mu)}.
k' = \frac{k}{1 + k(1 - \mu)}.
Note that when :math:`\alpha'/\alpha` goes to one, i.e., scattering is elastic,
the Klein-Nishina cross section becomes identical to the Thomson cross
section. In general though, the scattering is inelastic and is known as Compton
scattering. When a photon interacts with a bound electron in an atom, the
Klein-Nishina formula must be modified to account for the binding effects. As in
the case of coherent scattering, this is done by means of a form factor. The
differential cross section for incoherent scattering is given by
Note that when :math:`k'/k` goes to one, i.e., scattering is elastic, the
Klein-Nishina cross section becomes identical to the Thomson cross section. In
general though, the scattering is inelastic and is known as Compton scattering.
When a photon interacts with a bound electron in an atom, the Klein-Nishina
formula must be modified to account for the binding effects. As in the case of
coherent scattering, this is done by means of a form factor. The differential
cross section for incoherent scattering is given by
.. math::
:label: incoherent-xs
\frac{d\sigma}{d\mu} = \frac{d\sigma_{KN}}{d\mu} S(x,Z) = \pi r_e^2 \left (
\frac{\alpha'}{\alpha} \right )^2 \left [ \frac{\alpha'}{\alpha} +
\frac{\alpha}{\alpha'} + \mu^2 - 1 \right ] S(x,Z)
\frac{k'}{k} \right )^2 \left [ \frac{k'}{k} + \frac{k}{k'} + \mu^2 - 1
\right ] S(x,Z)
where :math:`S(x,Z)` is the form factor. The approach in OpenMC is to first
sample the Klein-Nishina cross section and then perform rejection sampling on
the form factor. As in other codes, `Kahn's rejection method`_ is used for
:math:`\alpha < 3` and a direct method by Koblinger_ is used for :math:`\alpha
\ge 3`. The complete algorithm is as follows:
:math:`k < 3` and a direct method by Koblinger_ is used for :math:`k \ge 3`.
The complete algorithm is as follows:
1. If :math:`\alpha < 3`, sample :math:`\mu` from the Klein-Nishina cross
section using Kahn's rejection method. Otherwise, use Koblinger's direct
method.
1. If :math:`k < 3`, sample :math:`\mu` from the Klein-Nishina cross section
using Kahn's rejection method. Otherwise, use Koblinger's direct method.
2. Calculate :math:`x` and :math:`\bar{x}` using :eq:`momentum-transfer` and
:eq:`xmax`, respectively.
@ -215,19 +215,370 @@ the form factor. As in other codes, `Kahn's rejection method`_ is used for
Doppler Energy Broadening
+++++++++++++++++++++++++
LA-UR-04-0487_ and LA-UR-04-0488_
Bound electrons are not at rest but have a momentum distribution that will
cause the energy of the scattered photon to be Doppler broadened. More tightly
bound electrons have a wider momentum distribution, so the energy spectrum of
photons scattering off inner shell electrons will be broadened the most.
In addition, scattering from bound electrons places a limit on the maximum
scattered photon energy:
.. math::
:label: max-energy-out
E'_{\text{max}} = E - E_{b,i},
where :math:`E_{b,i}` is the binding energy of the :math:`i`-th subshell.
Compton profiles :math:`J_i(p_z)` are used to account for the binding effects.
The quantity :math:`p_z = {\bf p} \cdot {\bf q}/q` is the projection of the
initial electron momentum on :math:`{\bf q}`, where the scattering vector
:math:`{\bf q} = {\bf p} - {\bf p'}` is the momentum gained by the photon,
:math:`{\bf p}` is the initial momentum of the electron, and :math:`{\bf p'}`
is the momentum of the scattered electron. Applying the conservation of energy
and momentum, :math:`p_z` can be written in terms of the photon energy and
scattering angle:
.. math::
:label: pz
p_z = \frac{E - E' - EE'(1 - \mu)/(m_e c^2)}{-\alpha \sqrt{E^2 + E'^2 -
2EE'\mu}},
where :math:`\alpha` is the fine structure constant. The maximum momentum
transferred, :math:`p_{z,\text{max}}`, can be calculated from :eq:`pz` using
:math:`E' = E'_{\text{max}}`. The Compton profile of the :math:`i`-th electron
subshell is defined as
.. math::
:label: compton-profile
J_i(p_z) = \int \int \rho_i({\bf p}) dp_x dp_y,
where :math:`\rho_i({\bf p})` is the initial electron momentum distribution.
:math:`J_i(p_z)` can be interpreted as the probability density function of
:math:`p_z`.
The Doppler broadened energy of the Compton-scattered photon can be sampled by
selecting an electron shell, sampling a value of :math:`p_z` using the Compton
profile, and calculating the scattered photon energy. The theory and methods
used to do this are described in detail in LA-UR-04-0487_ and LA-UR-04-0488_.
The sampling algorithm is summarized below:
1. Sample :math:`\mu` from :eq:`incoherent-xs` using the algorithm described in
:ref:`incoherent-sampling`.
2. Sample the electron subshell :math:`i` using the number of electrons per
shell as the probability mass function.
3. Sample :math:`p_z` using :math:`J_i(p_z)` as the PDF.
4. Calculate :math:`E'` by solving :eq:`pz` for :math:`E'` using the sampled
value of :math:`p_z`.
5. If :math:`p_z < p_{z,\text{max}}` for shell :math:`i`, accept :math:`E'`.
Otherwise repeat from step 2.
Compton Electrons
+++++++++++++++++
Because the Compton-scattered photons can transfer a large fraction of their
energy to the kinetic energy of the recoil electron, which may in turn go on to
lose its energy as bremsstrahlung radiation, it is necessary to accurately
model the angular and energy distributions of Compton electrons. The energy of
the recoil electron ejected from the :math:`i`-th subshell is given by
.. math::
:label: compton-electron-energy
E_{-} = E - E' - E_{b,i}.
The direction of the electron is assumed to be in the direction of the momentum
transfer, with the cosine of the polar angle given by
.. math::
:label: compton-electron-mu
\mu_{-} = \frac{E - E'\mu}{\sqrt{E^2 +E'^2 - 2EE'\mu}}
and the azimuthal angle :math:`\phi_{-} = \phi + \pi`, where :math:`\phi` is
the azimuthal angle of the photon. The vacancy left by the ejected electron is
filled through atomic relaxation.
Photoelectric Effect
--------------------
In the photoelectric effect, the incident photon is absorbed by an atomic
electron, which is then emitted from the :math:`i`-th shell with kinetic energy
.. math::
:label: photoelectron-kinetic-energy
E_{-} = E - E_{b,i}.
Photoelectric emission is only possible when the photon energy exceeds the
binding energy of the shell. These binding energies are often referred to as
edge energies because the otherwise continuously decreasing cross section has
discontinuities at these points, creating the characteristic sawtooth shape.
The photoelectric effect dominates at low energies and is more important for
heavier elements.
When simulating the photoelectric effect, the first step is to sample the
electron shell. The shell :math:`i` where the ionization occurs can be
considered a discrete random variable with probability mass function
.. math::
:label: photoelectron-shell-pdf
p_i = \frac{\sigma_{\text{pe},i}}{\sigma_{\text{pe}}},
where :math:`\sigma_{\text{pe},i}` is the cross section of the :math:`i`-th
shell, and the total photoelectric cross section of the atom,
:math:`\sigma_{\text{pe}}`, is the sum over the shell cross sections. Once the
shell has been sampled, the energy of the photoelectron is calculated using
:eq:`photoelectron-kinetic-energy`.
To determine the direction of the photoelectron, we implement the method
described in Kaltiaisenaho_, which models the angular distribution of the
photoelectrons using the K-shell cross section derived by Sauter (K-shell
electrons are the most tightly bound, and they contribute the most to
:math:`\sigma_{\text{pe}}`). The non-relativistic Sauter distribution for
unpolarized photons can be approximated as
.. math::
:label: sauter
\frac{d\sigma_{\text{pe}}}{d\mu_{-}} \propto
\frac{1 - \mu_{-}^2}{(1 - \beta_{-} \mu_{-})^4},
where :math:`\beta_{-}` is the ratio of the velocity of the electron to the
speed of light,
.. math::
:label: beta-2
\beta_{-} = \frac{\sqrt{(E_{-}(E_{-} + 2m_e c^2)}}{E_{-} + m_e c^2}.
To sample :math:`\mu_{-}` from the Sauter distribution, we first express
:eq:`sauter` in the form:
.. math::
:label: photoelectron-mu-pdf
f(\mu_{-}) = \frac{3}{2} \psi(\mu_{-}) g(\mu_{-}),
where
.. math::
:label: mu-pdf-factors
\psi(\mu_{-}) &= \frac{(1 - \beta_{-}^2)(1 - \mu_{-}^2)}{(1 -
\beta_{-}\mu_{-})^2}, \\
g(\mu_{-}) &= \frac{1 - \beta_{-}^2}{2 (1 - \beta_{-}\mu_{-})^2}.
In the interval :math:`[-1, 1]`, :math:`g(\mu_{-})` is a normalized PDF and
:math:`\psi(\mu_{-})` satisfies the condition :math:`0 < \psi(\mu_{-}) < 1`.
The following algorithm can now be used to sample :math:`\mu_{-}`:
1. Using the inverse transform method, sample :math:`\mu_{-}` from
:math:`g(\mu_{-})` using the sampling formula
.. math::
\mu_{-} = \frac{2\xi_1 + \beta_{-} - 1}{2\beta_{-}\xi_1 - \beta_{-} + 1}.
2. If :math:`\xi_2 \le \psi(\mu_{-})`, accept :math:`\mu_{-}`. Otherwise,
repeat the sampling from step 1.
The azimuthal angle is sampled uniformly on :math:`[0, 2\pi)`.
The atom is left in an excited state with a vacancy in the :math:`i`-th shell
and decays to its ground state through a cascade of transitions that produce
fluorescent photons and Auger electrons.
Pair Production
---------------
In electron-positron pair production, a photon is absorbed in the vicinity of
an atomic nucleus or an electron and an electron and positron are created. Pair
production is the dominant interaction with matter at high photon energies and
is more important for high-Z elements. When it takes place in the field of a
nucleus, energy is essentially conserved among the incident photon and the
resulting charged particles. Therefore, in order for pair production to occur,
the photon energy must be greater than the sum of the rest mass energies of the
electron and positron, i.e., :math:`E_{\text{threshold,pp}} = 2 m_e c^2 =
1.022` MeV.
The photon can also interact in the field of an atomic electron. This process
is referred to as "triplet production" because the target electron is ejected
from the atom and three charged particles emerge from the interaction. In this
case, the recoiling electron also absorbs some energy, so the energy threshold
for triplet production is greater than that of pair production from atomic
nuclei, with :math:`E_{\text{threshold,tp}} = 4 m_e c^2 = 2.044` MeV. The ratio
of the triplet production cross section to the pair production cross section is
approximately 1/Z, so triplet production becomes increasingly unimportant for
high-Z elements. Though it can be significant in lighter elements, the momentum
of the recoil electron becomes negligible in the energy regime where pair
production dominates. For our purposes, it is a good approximation to treat
triplet production as pair production and only simulate the electron-positron
pair.
Accurately modeling the creation of electron-positron pair is important because
the charged particles can go on to lose much of their energy as bremsstrahlung
radiation, and the subsequent annihilation of the positron with an electron
produces two additional photons. We sample the energy and direction of the
charged particles using a semiempirical model described in Salvat_. The
Bethe-Heitler differential cross section, given by
.. math::
:label: bethe-heitler
\frac{d\sigma_{\text{pp}}}{d\epsilon} = \alpha r_e^2 Z^2
\left[ (\epsilon^2 + (1-\epsilon)^2) (\Phi_1 - 4f_C) +
\frac{2}{3}\epsilon(1-\epsilon)(\Phi_2 - 4f_C) \right],
is used as a starting point, where :math:`\alpha` is the fine structure
constant, :math:`f_C` is the Coulomb correction function, :math:`\Phi_1` and
:math:`\Phi_2` are screening functions, and :math:`\epsilon = (E_{-} + m_e
c^2)/E` is the electron reduced energy (i.e., the fraction of the photon energy
given to the electron). :math:`\epsilon` can take values between
:math:`\epsilon_{\text{min}} = k^{-1}` (when the kinetic energy of the electron
is zero) and :math:`\epsilon_{\text{max}} = 1 - k^{-1}` (when the kinetic
energy of the positron is zero).
The Coulomb correction, given by
.. math::
:label: coulomb-correction
f_C = \alpha^{2}Z^{2} \big[&(1 + \alpha^{2}Z^{2})^{-1} + 0.202059
- 0.03693\alpha^{2}Z^{2} + 0.00835\alpha^{4}Z^{4} \\
&- 0.00201\alpha^{6}Z^{6} + 0.00049\alpha^{8}Z^{8}
- 0.00012\alpha^{10}Z^{10} + 0.00003\alpha^{12}Z^{12}\big]
is introduced to correct for the fact that the Bethe-Heitler differential cross
section was derived using the Born approximation, which treats the Coulomb
interaction as a small perturbation.
The screening functions :math:`\Phi_1` and :math:`\Phi_2` account for the
screening of the Coulomb field of the atomic nucleus by outer electrons. Since
they are given by integrals which include the atomic form factor, they must be
computed numerically for a realistic form factor. However, by assuming
exponential screening and using a simplified form factor, analytical
approximations of the screening functions can be derived:
.. math::
:label: screening-functions
\Phi_1 &= 2 - 2\ln(1 + b^2) - 4b\arctan(b^{-1}) + 4\ln(Rm_{e}c/\hbar) \\
\Phi_2 &= \frac{4}{3} - 2\ln(1 + b^2) + 2b^2 \left[ 4 - 4b\arctan(b^{-1})
- 3\ln(1 + b^{-2}) \right] + 4\ln(Rm_{e}c/\hbar)
where
.. math::
:label: b
b = \frac{Rm_{e}c}{2k\epsilon(1 - \epsilon)\hbar}.
and :math:`R` is the screening radius.
The differential cross section in :eq:`bethe-heitler` with the approximations
described above will not be accurate at low energies: the lower boundary of
:math:`\epsilon` will be shifted above :math:`\epsilon_{\text{min}}` and the
upper boundary of :math:`\epsilon` will be shifted below
:math:`\epsilon_{\text{max}}`. To offset this behavior, a correcting factor
:math:`F_0(k, Z)` is used:
.. math::
:label: correcting-factor
F_0(k, Z) =~& (0.1774 + 12.10\alpha Z - 11.18\alpha^{2}Z^{2})(2/k)^{1/2} \\
&+ (8.523 + 73.26\alpha Z - 44.41\alpha^{2}Z^{2})(2/k) \\
&- (13.52 + 121.1\alpha Z - 96.41\alpha^{2}Z^{2})(2/k)^{3/2} \\
&+ (8.946 + 62.05\alpha Z - 63.41\alpha^{2}Z^{2})(2/k)^{2}.
To aid sampling, the differential cross section used to sample :math:`\epsilon`
(minus the normalization constant) can now be expressed in the form
.. math::
:label: pp-pdf
\frac{d\sigma_{\text{pp}}}{d\epsilon} =
u_1 \frac{\phi_1(\epsilon)}{\phi_1(1/2)} \pi_1(\epsilon)
+ u_2 \frac{\phi_2(\epsilon)}{\phi_2(1/2)} \pi_2(\epsilon)
where
.. math::
:label: u
u_1 &= \frac{2}{3} \left(\frac{1}{2} - \frac{1}{k}\right)^2 \phi_1(1/2), \\
u_2 &= \phi_2(1/2),
.. math::
:label: phi
\phi_1(\epsilon) &= \frac{1}{2}(3\Phi_1 - \Phi_2) - 4f_{C}(Z) + F_0(k, Z), \\
\phi_2(\epsilon) &= \frac{1}{4}(3\Phi_1 + \Phi_2) - 4f_{C}(Z) + F_0(k, Z),
and
.. math::
:label: pi
\pi_1(\epsilon) &= \frac{3}{2} \left(\frac{1}{2} - \frac{1}{k}\right)^{-3}
\left(\frac{1}{2} - \epsilon\right)^2, \\
\pi_2(\epsilon) &= \frac{1}{2} \left(\frac{1}{2} - \frac{1}{k}\right)^{-1}.
The functions in :eq:`phi` are non-negative and maximum at :math:`\epsilon =
1/2`. In the interval :math:`(\epsilon_{\text{min}}, \epsilon_{\text{max}})`,
the functions in :eq:`pi` are normalized PDFs and
:math:`\phi_i(\epsilon)/\phi_i(1/2)` satisfies the condition :math:`0 <
\phi_i(\epsilon)/\phi_i(1/2) < 1`. The following algorithm can now be used to
sample the reduced electron energy :math:`\epsilon`:
1. Sample :math:`i` according to the point probabilities
:math:`p(i=1) = u_1/(u_1 + u_2)` and :math:`p(i=2) = u_2/(u_1 + u_2)`.
2. Using the inverse transform method, sample :math:`\epsilon` from
:math:`\pi_i(\epsilon)` using the sampling formula
.. math::
\epsilon &= \frac{1}{2} + \left(\frac{1}{2} - \frac{1}{k}\right)
(2\xi_1 - 1)^{1/3} ~~~~&\text{if}~~ i = 1 \\
\epsilon &= \frac{1}{k} + \left(\frac{1}{2} -
\frac{1}{k}\right) 2\xi_1 ~~~~&\text{if}~~ i = 2.
3. If :math:`\xi_2 \le \phi_i(\epsilon)/\phi_i(1/2)`, accept
:math:`\epsilon`. Otherwise, repeat the sampling from step 1.
Because charged particles have a much smaller range than the mean free path of
photons and because they immediately undergo multiple scattering events which
randomize their direction, it is sufficient to use a simplified model to sample
the direction of the electron and positron. The cosines of the polar angles are
sampled using the leading order term of the Sauter–Gluckstern–Hull
distribution,
.. math::
:label: sauter–gluckstern–hull
p(\mu_{\pm}) = C(1 - \beta_{\pm}\mu_{\pm})^{-2},
where :math:`C` is a normalization constant and :math:`\beta_{\pm}` is the
ratio of the velocity of the charged particle to the speed of light given in
:eq:`beta-2`.
The inverse transform method is used to sample :math:`\mu_{-}` and
:math:`\mu_{+}` from :eq:`sauter–gluckstern–hull`, using the sampling formula
.. math::
:label: sample-mu
\mu_{\pm} = \frac{2\xi - 1 + \beta_{\pm}}{(2\xi - 1)\beta_{\pm} + 1}.
The azimuthal angles for the electron and positron are sampled independently
and uniformly on :math:`[0, 2\pi)`.
-------------------
Secondary Processes
@ -248,19 +599,297 @@ additional photons.
Atomic Relaxation
-----------------
When an electron is ejected from an atom and a vacancy is left in an inner
shell, an electron from a higher energy level will fill the vacancy. This
results in either a radiative transition, in which a photon with a
characteristic energy (fluorescence photon) is emitted, or non-radiative
transition, in which an electron from a shell that is farther out (Auger
electron) is emitted. If a non-radiative transition occurs, the new vacancy is
filled in the same manner, and as the process repeats a shower of photons and
electrons can be produced.
The energy of a fluorescence photon is the equal to the energy difference
between the transition states, i.e.,
.. math::
:label: fluorescence-photon-energy
E = E_{b,v} - E_{b,i},
where :math:`E_{b,v}` is the binding energy of the vacancy shell and
:math:`E_{b,i}` is the binding energy of the shell from which the electron
transitioned. The energy of an Auger electron is given by
.. math::
:label: auger-electron-energy
E_{-} = E_{b,v} - E_{b,i} - E_{b,a},
where :math:`E_{b,a}` is the binding energy of the shell from which the Auger
electron is emitted. While Auger electrons are low-energy so their range and
bremsstrahlung yield is small, fluorescence photons can travel far before
depositing their energy, so the relaxation process should be modeled in detail.
Transition energies and probabilities are needed for each subshell to simulate
atomic relaxation. Starting with the initial shell vacancy, the following
recursive algorithm is used to fill vacancies and create fluorescence photons
and Auger electrons:
1. If there are no transitions for the vacancy shell, create a fluorescence
photon assuming it is from a captured free electron and terminate.
2. Sample a transition using the transition probabilities for the vacancy
shell as the probability mass function.
3. Create either a fluorescence photon or Auger electron, sampling the
direction of the particle isotropically.
4. If a non-radiative transition occurred, repeat from step 1 for the vacancy
left by the emitted Auger electron.
5. Repeat from step 1 for vacancy left by the transition electron.
Electron-Positron Annihilation
------------------------------
When a positron collides with an electron, both particles are annihilated and
generally two photons with equal energy are created. If the kinetic energy of
the positron is high enough, the two photons can have different energies, and
the higher-energy photon is emitted preferentially in the direction of flight
of the positron. It is also possible to produce a single photon if the
interaction occurs with a bound electron, and in some cases three (or, rarely,
even more) photons can be emitted. However, the annihilation cross section is
largest for low-energy positrons, and as the positron energy decreases, the
angular distribution of the emitted photons becomes isotropic.
In OpenMC, we assume the most likely case in which a low-energy positron (which
has already lost most of its energy to bremsstrahlung radiation) interacts with
an electron which is free and at rest. Two photons with energy equal to the
electron rest mass energy :math:`m_e c^2 = 0.511` MeV are emitted isotropically
in opposite directions.
Bremsstrahlung
--------------
When a charged particle is decelerated in the field of an atom, some of its
kinetic energy is converted into electromagnetic radiation known as
bremsstrahlung, or 'braking radiation'. In each event, an electron or positron
with kinetic energy :math:`T` generates a photon with an energy :math:`E`
between :math:`0` and :math:`T`. Bremsstrahlung is described by a cross section
that is differential in photon energy, in the direction of the emitted photon,
and in the final direction of the charged particle. However, in Monte Carlo
simulations it is typical to integrate over the angular variables to obtain a
single differential cross section with respect to photon energy, which is often
expressed in the form
.. math::
:label: bremsstrahlung-dcs
\frac{d\sigma_{\text{br}}}{dE} = \frac{Z^2}{\beta^2} \frac{1}{E}
\chi(Z, T, \kappa),
where :math:`\kappa = E/T` is the reduced photon energy and :math:`\chi(Z, T,
\kappa)` is the scaled bremsstrahlung cross section, which is experimentally
measured.
Because electrons are attracted to atomic nuclei whereas positrons are
repulsed, the cross section for positrons is smaller, though it approaches that
of electrons in the high energy limit. To obtain the positron cross section, we
multiply :eq:`bremsstrahlung-dcs` by the :math:`\kappa`-independent factor used
in Salvat_,
.. math::
:label: positron-factor
F_{\text{p}}(Z,T) =
& 1 - \text{exp}(-1.2359\times 10^{-1}t + 6.1274\times 10^{-2}t^2 - 3.1516\times 10^{-2}t^3 \\
& + 7.7446\times 10^{-3}t^4 - 1.0595\times 10^{-3}t^5 + 7.0568\times 10^{-5}t^6 \\
& - 1.8080\times 10^{-6}t^7),
where
.. math::
:label: positron-factor-t
t = \ln\left(1 + \frac{10^6}{Z^2}\frac{T}{\text{m}_\text{e}c^2} \right).
:math:`F_{\text{p}}(Z,T)` is the ratio of the radiative stopping powers for
positrons and electrons. Stopping power describes the average energy loss per
unit path length of a charged particle as it passes through matter:
.. math::
:label: stopping-power
-\frac{dT}{ds} = n \int E \frac{d\sigma}{dE} dE \equiv S(T),
where :math:`n` is the number density of the material and :math:`d\sigma/dE` is
the cross section differential in energy loss. The total stopping power
:math:`S(T)` can be separated into two components: the radiative stopping
power :math:`S_{\text{rad}}(T)`, which refers to energy loss due to
bremsstrahlung, and the collision stopping power :math:`S_{\text{col}}(T)`,
which refers to the energy loss due to inelastic collisions with bound
electrons in the material that result in ionization and excitation. To obtain
the radiative stopping power for positrons, the radiative stopping power for
electrons is multiplied by :eq:`positron-factor`. Currently, the collision
stopping power for electrons is also used for positrons.
While the models for photon interactions with matter described above can safely
assume interactions occur with free atoms, sampling the target atom based on
the macroscopic cross sections, molecular effects cannot necessarily be
disregarded for charged particle treatment. For compounds and mixtures, the
bremsstrahlung cross section is calculated using Bragg's additivity rule as
.. math::
:label: material-bremsstrahlung-dcs
\frac{d\sigma_{\text{br}}}{dE} = \frac{1}{\beta^2 E} \sum_i \gamma_i Z^2_i
\chi(Z_i, T, \kappa),
where the sum is over the constituent elements and :math:`\gamma_i` is the
atomic fraction of the :math:`i`-th element. Similarly, the radiative stopping
power is calculated using Bragg's additivity rule as
.. math::
:label: material-radiative-stopping-power
S_{\text{rad}}(T) = \sum_i w_i S_{\text{rad},i}(T),
where :math:`w_i` is the mass fraction of the :math:`i`-th element. The
collision stopping power, however, is a function of certain quantities such as
the mean excitation energy :math:`I` and the density effect correction
:math:`\delta_F` that depend on molecular properties. These quantities cannot
simply be summed over constituent elements in a compound, but should instead be
calculated for the material. Currently, we use Bragg's additivity rule to
calculate the collision stopping power as well, but this is not a good
approximation and should be fixed in the future.
.. _ttb:
Thick-Target Bremsstrahlung Approximation
+++++++++++++++++++++++++++++++++++++++++
Since charged particles lose their energy on a much shorter distance scale than
neutral particles, not much error should be introduced by neglecting to
transport electrons. However, the bremsstrahlung emitted from high energy
electrons and positrons can travel far from the interaction site. Thus, even
without a full electron transport mode it is necessary to model bremsstrahlung.
We use a thick-target bremsstrahlung (TTB) approximation based on the models in
Salvat_ and Kaltiaisenaho_ for generating bremsstrahlung photons, which assumes
the charged particle loses all its energy in a single homogeneous material
region.
To model bremsstrahlung using the TTB approximation, we need to know the number
of photons emitted by the charged particle and the energy distribution of the
photons. These quantities can be calculated using the continuous slowing down
approximation (CSDA). The CSDA assumes charged particles lose energy
continuously along their trajectory with a rate of energy loss equal to the
total stopping power, ignoring fluctuations in the energy loss. The
approximation is useful for expressing average quantities that describe how
charged particles slow down in matter. For example, the CSDA range approximates
the average path length a charged particle travels as it slows to rest:
.. math::
:label: csda-range
R(T) = \int^T_0 \frac{dT'}{S(T')}.
Actual path lengths will fluctuate around :math:`R(T)`. The average number of
photons emitted per unit path length is given by the inverse bremsstrahlung
mean free path:
.. math::
:label: inverse-bremsstrahlung-mfp
\lambda_{\text{br}}^{-1}(T,E_{\text{cut}})
= n\int_{E_{\text{cut}}}^T\frac{d\sigma_{\text{br}}}{dE}dE
= n\frac{Z^2}{\beta^2}\int_{\kappa_{\text{cut}}}^1\frac{1}{\kappa}
\chi(Z,T,\kappa)d\kappa.
The lower limit of the integral in :eq:`inverse-bremsstrahlung-mfp` is non-zero
because the bremsstrahlung differential cross section diverges for small photon
energies but is finite for photon energies above some cutoff energy
:math:`E_{\text{cut}}`. The mean free path
:math:`\lambda_{\text{br}}^{-1}(T,E_{\text{cut}})` is used to calculate the
photon number yield, defined as the average number of photons emitted with
energy greater than :math:`E_{\text{cut}}` as the charged particle slows down
from energy :math:`T` to :math:`E_{\text{cut}}`. The photon number yield is
given by
.. math::
:label: photon-number-yield
Y(T,E_{\text{cut}}) = \int^{R(T)}_{R(E_{\text{cut}})}
\lambda_{\text{br}}^{-1}(T',E_{\text{cut}})ds = \int_{E_{\text{cut}}}^T
\frac{\lambda_{\text{br}}^{-1}(T',E_{\text{cut}})}{S(T')}dT'.
:math:`Y(T,E_{\text{cut}})` can be used to construct the energy spectrum of
bremsstrahlung photons: the number of photons created with energy between
:math:`E_1` and :math:`E_2` by a charged particle with initial kinetic energy
:math:`T` as it comes to rest is given by :math:`Y(T,E_1) - Y(T,E_2)`.
To simulate the emission of bremsstrahlung photons, the total stopping power
and bremsstrahlung differential cross section for positrons and electrons must
be calculated for a given material using :eq:`material-bremsstrahlung-dcs` and
:eq:`material-radiative-stopping-power`. These quantities are used to build the
tabulated bremsstrahlung energy PDF and CDF for that material for each incident
energy :math:`T_k` on the energy grid. The following algorithm is then applied
to sample the photon energies:
1. For an incident charged particle with energy :math:`T`, sample the number of
emitted photons as
.. math::
N = \lfloor Y(T,E_{\text{cut}}) + \xi_1 \rfloor.
2. Rather than interpolate the PDF between indices :math:`k` and :math:`k+1`
for which :math:`T_k < T < T_{k+1}`, which is computationally expensive, use
the composition method and sample from the PDF at either :math:`k` or
:math:`k+1`. Using linear interpolation on a logarithmic scale, the PDF can
be expressed as
.. math::
p_{\text{br}}(T,E) = \pi_k p_{\text{br}}(T_k,E) + \pi_{k+1}
p_{\text{br}}(T_{k+1},E),
where the interpolation weights are
.. math::
\pi_k = \frac{\ln T_{k+1} - \ln T}{\ln T_{k+1} - \ln T_k},~~~
\pi_{k+1} = \frac{\ln T - \ln T_k}{\ln T_{k+1} - \ln T_k}.
Sample either the index :math:`i = k` or :math:`i = k+1` according to the
point probabilities :math:`\pi_{k}` and :math:`\pi_{k+1}`.
3. Determine the maximum value of the CDF :math:`P_{\text{br,max}}`.
3. Sample the photon energies using the inverse transform method with the
tabulated CDF :math:`P_{\text{br}}(T_i, E)` i.e.,
.. math::
E = E_j \left[ (1 + a_j) \frac{\xi_2 P_{\text{br,max}} -
P_{\text{br}}(T_i, E_j)} {E_j p_{\text{br}}(T_i, E_j)} + 1
\right]^{\frac{1}{1 + a_j}}
where the interpolation factor :math:`a_j` is given by
.. math::
a_j = \frac{\ln p_{\text{br}}(T_i,E_{j+1}) - \ln p_{\text{br}}(T_i,E_j)}
{\ln E_{j+1} - \ln E_j}
and :math:`P_{\text{br}}(T_i, E_j) \le \xi_2 P_{\text{br,max}} \le
P_{\text{br}}(T_i, E_{j+1})`.
We ignore the range of the electron or positron, i.e., the bremsstrahlung
photons are produced in the same location that the charged particle was
created. The direction of the photons is assumed to be the same as the
direction of the incident charged particle, which is a reasonable approximation
at higher energies when the bremsstrahlung radiation is emitted at small
angles.
.. _Koblinger: https://doi.org/10.13182/NSE75-A26663
@ -273,3 +902,7 @@ Thick-Target Bremsstrahlung Approximation
.. _LA-UR-04-0487: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-04-0487.pdf
.. _LA-UR-04-0488: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-04-0488.pdf
.. _Kaltiaisenaho: https://aaltodoc.aalto.fi/bitstream/handle/123456789/21004/master_Kaltiaisenaho_Toni_2016.pdf
.. _Salvat: http://www.oecd-nea.org/globalsearch/download.php?doc=77434

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

@ -39,7 +39,6 @@ Core Functions
openmc.data.linearize
openmc.data.thin
openmc.data.water_density
openmc.data.write_compact_458_library
openmc.data.zam
Angle-Energy Distributions

23
docs/source/quickinstall.rst
View file

@ -58,6 +58,29 @@ are no longer supported.
.. _Personal Package Archive: https://launchpad.net/~paulromano/+archive/staging
.. _APT package manager: https://help.ubuntu.com/community/AptGet/Howto
-------------------------------------------
Installing on Linux/Mac/Windows with Docker
-------------------------------------------
OpenMC can be easily deployed using `Docker <https://www.docker.com/>`_ on any
Windows, Mac or Linux system. With Docker running, execute the following
command in the shell to download and run a `Docker image`_ with the most recent release of OpenMC from `DockerHub <https://hub.docker.com/>`_ called ``openmc/openmc:v0.10.0``:
.. code-block:: sh
docker run openmc/openmc:v0.10.0
This will take several minutes to run depending on your internet download speed. The command will place you in an interactive shell running in a `Docker container`_ with OpenMC installed.
.. note:: The ``docker run`` command supports many `options`_ for spawning
containers -- including `mounting volumes`_ from the host
filesystem -- which many users will find useful.
.. _Docker image: https://docs.docker.com/engine/reference/commandline/images/
.. _Docker container: https://www.docker.com/resources/what-container
.. _options: https://docs.docker.com/engine/reference/commandline/run/
.. _mounting volumes: https://docs.docker.com/storage/volumes/
---------------------------------------
Installing from Source on Ubuntu 15.04+
---------------------------------------

48
include/openmc/cell.h
View file

@ -13,6 +13,9 @@
#include "openmc/constants.h"
#include "openmc/position.h"
#ifdef DAGMC
#include "DagMC.hpp"
#endif
namespace openmc {
@ -107,9 +110,8 @@ public:
std::vector<int32_t> offset_; //!< Distribcell offset table
Cell() {};
explicit Cell(pugi::xml_node cell_node);
Cell() {};
//! \brief Determine if a cell contains the particle at a given location.
//!
@ -130,21 +132,57 @@ public:
//! \param on_surface The signed index of a surface that the coordinate is
//! known to be on. This index takes precedence over surface sense
//! calculations.
virtual bool
contains(Position r, Direction u, int32_t on_surface) const = 0;
//! Find the oncoming boundary of this cell.
virtual std::pair<double, int32_t>
distance(Position r, Direction u, int32_t on_surface) const = 0;
//! Write all information needed to reconstruct the cell to an HDF5 group.
//! @param group_id An HDF5 group id.
virtual void to_hdf5(hid_t group_id) const = 0;
virtual ~Cell() {}
};
class CSGCell : public Cell
{
public:
CSGCell();
explicit CSGCell(pugi::xml_node cell_node);
bool
contains(Position r, Direction u, int32_t on_surface) const;
//! Find the oncoming boundary of this cell.
std::pair<double, int32_t>
distance(Position r, Direction u, int32_t on_surface) const;
//! \brief Write cell information to an HDF5 group.
//! \param group_id An HDF5 group id.
void to_hdf5(hid_t group_id) const;
protected:
bool contains_simple(Position r, Direction u, int32_t on_surface) const;
bool contains_complex(Position r, Direction u, int32_t on_surface) const;
};
#ifdef DAGMC
class DAGCell : public Cell
{
public:
moab::DagMC* dagmc_ptr_;
DAGCell();
std::pair<double, int32_t> distance(Position r, Direction u, int32_t on_surface) const;
bool contains(Position r, Direction u, int32_t on_surface) const;
void to_hdf5(hid_t group_id) const;
};
#endif
} // namespace openmc
#endif // OPENMC_CELL_H

28
include/openmc/dagmc.h Normal file
View file

@ -0,0 +1,28 @@
#ifndef OPENMC_DAGMC_H
#define OPENMC_DAGMC_H
#ifdef DAGMC
#include "DagMC.hpp"
#include "openmc/cell.h"
#include "openmc/surface.h"
namespace openmc {
extern moab::DagMC* DAG;
extern "C" void load_dagmc_geometry();
extern "C" void free_memory_dagmc();
}
#endif // DAGMC
#endif // OPENMC_DAGMC_H
#ifdef DAGMC
extern "C" constexpr bool dagmc_enabled = true;
#else
extern "C" constexpr bool dagmc_enabled = false;
#endif

1
include/openmc/position.h
View file

@ -25,6 +25,7 @@ struct Position {
Position& operator-=(double);
Position& operator*=(Position);
Position& operator*=(double);
const double& operator[](int i) const {
switch (i) {
case 0: return x;

2
include/openmc/settings.h
View file

@ -45,6 +45,7 @@ extern "C" bool ufs_on; //!< uniform fission site method on?
extern "C" bool urr_ptables_on; //!< use unresolved resonance prob. tables?
extern "C" bool write_all_tracks; //!< write track files for every particle?
extern "C" bool write_initial_source; //!< write out initial source file?
extern "C" bool dagmc; //!< indicator of DAGMC geometry
// Paths to various files
extern std::string path_cross_sections; //!< path to cross_sections.xml
@ -86,7 +87,6 @@ extern "C" int trigger_batch_interval; //!< Batch interval for triggers
extern "C" int verbosity; //!< How verbose to make output
extern "C" double weight_cutoff; //!< Weight cutoff for Russian roulette
extern "C" double weight_survive; //!< Survival weight after Russian roulette
} // namespace settings
//! Read settings from XML file

50
include/openmc/surface.h
View file

@ -12,6 +12,9 @@
#include "openmc/constants.h"
#include "openmc/position.h"
#ifdef DAGMC
#include "DagMC.hpp"
#endif
namespace openmc {
@ -65,6 +68,7 @@ public:
std::vector<int> neighbor_neg_; //!< List of cells on negative side
explicit Surface(pugi::xml_node surf_node);
Surface();
virtual ~Surface() {}
@ -104,12 +108,40 @@ public:
//! Write all information needed to reconstruct the surface to an HDF5 group.
//! \param group_id An HDF5 group id.
//TODO: this probably needs to include i_periodic for PeriodicSurface
virtual void to_hdf5(hid_t group_id) const = 0;
};
class CSGSurface : public Surface
{
public:
explicit CSGSurface(pugi::xml_node surf_node);
CSGSurface();
void to_hdf5(hid_t group_id) const;
protected:
virtual void to_hdf5_inner(hid_t group_id) const = 0;
};
//==============================================================================
//! A `Surface` representing a DAGMC-based surface in DAGMC.
//==============================================================================
#ifdef DAGMC
class DAGSurface : public Surface
{
public:
moab::DagMC* dagmc_ptr_;
DAGSurface();
double evaluate(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
//! Get the bounding box of this surface.
BoundingBox bounding_box() const;
void to_hdf5(hid_t group_id) const;
};
#endif
//==============================================================================
//! A `Surface` that supports periodic boundary conditions.
//!
@ -118,7 +150,7 @@ protected:
//! `XPlane`-`YPlane` pairs.
//==============================================================================
class PeriodicSurface : public Surface
class PeriodicSurface : public CSGSurface
{
public:
int i_periodic_{C_NONE}; //!< Index of corresponding periodic surface
@ -227,7 +259,7 @@ public:
//! \f$(y - y_0)^2 + (z - z_0)^2 - R^2 = 0\f$
//==============================================================================
class SurfaceXCylinder : public Surface
class SurfaceXCylinder : public CSGSurface
{
double y0_, z0_, radius_;
public:
@ -245,7 +277,7 @@ public:
//! \f$(x - x_0)^2 + (z - z_0)^2 - R^2 = 0\f$
//==============================================================================
class SurfaceYCylinder : public Surface
class SurfaceYCylinder : public CSGSurface
{
double x0_, z0_, radius_;
public:
@ -263,7 +295,7 @@ public:
//! \f$(x - x_0)^2 + (y - y_0)^2 - R^2 = 0\f$
//==============================================================================
class SurfaceZCylinder : public Surface
class SurfaceZCylinder : public CSGSurface
{
double x0_, y0_, radius_;
public:
@ -281,7 +313,7 @@ public:
//! \f$(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 - R^2 = 0\f$
//==============================================================================
class SurfaceSphere : public Surface
class SurfaceSphere : public CSGSurface
{
double x0_, y0_, z0_, radius_;
public:
@ -299,7 +331,7 @@ public:
//! \f$(y - y_0)^2 + (z - z_0)^2 - R^2 (x - x_0)^2 = 0\f$
//==============================================================================
class SurfaceXCone : public Surface
class SurfaceXCone : public CSGSurface
{
double x0_, y0_, z0_, radius_sq_;
public:
@ -317,7 +349,7 @@ public:
//! \f$(x - x_0)^2 + (z - z_0)^2 - R^2 (y - y_0)^2 = 0\f$
//==============================================================================
class SurfaceYCone : public Surface
class SurfaceYCone : public CSGSurface
{
double x0_, y0_, z0_, radius_sq_;
public:
@ -335,7 +367,7 @@ public:
//! \f$(x - x_0)^2 + (y - y_0)^2 - R^2 (z - z_0)^2 = 0\f$
//==============================================================================
class SurfaceZCone : public Surface
class SurfaceZCone : public CSGSurface
{
double x0_, y0_, z0_, radius_sq_;
public:
@ -352,7 +384,7 @@ public:
//! \f$A x^2 + B y^2 + C z^2 + D x y + E y z + F x z + G x + H y + J z + K = 0\f$
//==============================================================================
class SurfaceQuadric : public Surface
class SurfaceQuadric : public CSGSurface
{
// Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0
double A_, B_, C_, D_, E_, F_, G_, H_, J_, K_;

4
openmc/capi/__init__.py
View file

@ -12,7 +12,7 @@ objects in the :mod:`openmc.capi` subpackage, for example:
"""
from ctypes import CDLL
from ctypes import CDLL, c_bool
import os
import sys
@ -38,6 +38,8 @@ else:
from unittest.mock import Mock
_dll = Mock()
dagmc_enabled = bool(c_bool.in_dll(_dll, "dagmc_enabled"))
from .error import *
from .core import *
from .nuclide import *

4
openmc/data/fission_energy.py
View file

@ -43,10 +43,6 @@ class FissionEnergyRelease(EqualityMixin):
dataset. The :meth:`FissionEnergyRelease.from_hdf5` method builds this
class from the usual OpenMC HDF5 data files.
:meth:`FissionEnergyRelease.from_endf` uses ENDF-formatted data.
:meth:`FissionEnergyRelease.from_compact_hdf5` uses a different HDF5 format
that is meant to be compact and store the exact same data as the ENDF
format. Files with this format can be generated with the
:func:`openmc.data.write_compact_458_library` function.
References
----------

31
openmc/filter.py
View file

@ -815,6 +815,7 @@ class MeshSurfaceFilter(MeshFilter):
mesh_key = 'mesh {}'.format(self.mesh.id)
# Find mesh dimensions - use 3D indices for simplicity
n_surfs = 4 * len(self.mesh.dimension)
if len(self.mesh.dimension) == 3:
nx, ny, nz = self.mesh.dimension
elif len(self.mesh.dimension) == 2:
@ -826,31 +827,33 @@ class MeshSurfaceFilter(MeshFilter):
# Generate multi-index sub-column for x-axis
filter_bins = np.arange(1, nx + 1)
repeat_factor = 12 * stride
repeat_factor = n_surfs * stride
filter_bins = np.repeat(filter_bins, repeat_factor)
tile_factor = data_size // len(filter_bins)
filter_bins = np.tile(filter_bins, tile_factor)
filter_dict[(mesh_key, 'x')] = filter_bins
# Generate multi-index sub-column for y-axis
filter_bins = np.arange(1, ny + 1)
repeat_factor = 12 * nx * stride
filter_bins = np.repeat(filter_bins, repeat_factor)
tile_factor = data_size // len(filter_bins)
filter_bins = np.tile(filter_bins, tile_factor)
filter_dict[(mesh_key, 'y')] = filter_bins
if len(self.mesh.dimension) > 1:
filter_bins = np.arange(1, ny + 1)
repeat_factor = n_surfs * nx * stride
filter_bins = np.repeat(filter_bins, repeat_factor)
tile_factor = data_size // len(filter_bins)
filter_bins = np.tile(filter_bins, tile_factor)
filter_dict[(mesh_key, 'y')] = filter_bins
# Generate multi-index sub-column for z-axis
filter_bins = np.arange(1, nz + 1)
repeat_factor = 12 * nx * ny * stride
filter_bins = np.repeat(filter_bins, repeat_factor)
tile_factor = data_size // len(filter_bins)
filter_bins = np.tile(filter_bins, tile_factor)
filter_dict[(mesh_key, 'z')] = filter_bins
if len(self.mesh.dimension) > 2:
filter_bins = np.arange(1, nz + 1)
repeat_factor = n_surfs * nx * ny * stride
filter_bins = np.repeat(filter_bins, repeat_factor)
tile_factor = data_size // len(filter_bins)
filter_bins = np.tile(filter_bins, tile_factor)
filter_dict[(mesh_key, 'z')] = filter_bins
# Generate multi-index sub-column for surface
repeat_factor = stride
filter_bins = np.repeat(_CURRENT_NAMES, repeat_factor)
filter_bins = np.repeat(_CURRENT_NAMES[:n_surfs], repeat_factor)
tile_factor = data_size // len(filter_bins)
filter_bins = np.tile(filter_bins, tile_factor)
filter_dict[(mesh_key, 'surf')] = filter_bins

5
openmc/geometry.py
View file

@ -266,7 +266,10 @@ class Geometry(object):
Dictionary mapping cell IDs to :class:`openmc.Cell` instances
"""
return self.root_universe.get_all_cells()
if self.root_universe is not None:
return self.root_universe.get_all_cells()
else:
return []
def get_all_universes(self):
"""Return all universes in the geometry.

21
openmc/settings.py
View file

@ -37,6 +37,8 @@ class Settings(object):
weight assigned to particles that are not killed after Russian
roulette. Value of energy should be a float indicating energy in eV
below which particle type will be killed.
dagmc : bool
Indicate that a CAD-based DAGMC geometry will be used.
electron_treatment : {'led', 'ttb'}
Whether to deposit all energy from electrons locally ('led') or create
secondary bremsstrahlung photons ('ttb').
@ -224,6 +226,8 @@ class Settings(object):
self._create_fission_neutrons = None
self._log_grid_bins = None
self._dagmc = None
@property
def run_mode(self):
return self._run_mode
@ -368,6 +372,10 @@ class Settings(object):
def log_grid_bins(self):
return self._log_grid_bins
@property
def dagmc(self):
return self._dagmc
@run_mode.setter
def run_mode(self, run_mode):
cv.check_value('run mode', run_mode, _RUN_MODES)
@ -511,6 +519,11 @@ class Settings(object):
cv.check_type('photon transport', photon_transport, bool)
self._photon_transport = photon_transport
@dagmc.setter
def dagmc(self, dagmc):
cv.check_type('dagmc geometry', dagmc, bool)
self._dagmc = dagmc
@ptables.setter
def ptables(self, ptables):
cv.check_type('probability tables', ptables, bool)
@ -945,6 +958,11 @@ class Settings(object):
elem = ET.SubElement(root, "log_grid_bins")
elem.text = str(self._log_grid_bins)
def _create_dagmc_subelement(self, root):
if self._dagmc is not None:
elem = ET.SubElement(root, "dagmc")
element.text = str(self._dagmc).lower()
def export_to_xml(self, path='settings.xml'):
"""Export simulation settings to an XML file.
@ -992,7 +1010,8 @@ class Settings(object):
self._create_volume_calcs_subelement(root_element)
self._create_create_fission_neutrons_subelement(root_element)
self._create_log_grid_bins_subelement(root_element)
self._create_dagmc_subelement(root_element)
# Clean the indentation in the file to be user-readable
clean_indentation(root_element)

3
openmc/summary.py
View file

@ -97,6 +97,9 @@ class Summary(object):
self._macroscopics = name.decode()
def _read_geometry(self):
if "dagmc" in self._f['geometry'].attrs.keys():
return
# Read in and initialize the Materials and Geometry
self._read_materials()
self._read_surfaces()

9
scripts/openmc-convert-mcnp71-data
View file

@ -25,8 +25,6 @@ parser = argparse.ArgumentParser(
)
parser.add_argument('-d', '--destination', default='mcnp_endfb71',
help='Directory to create new library in')
parser.add_argument('-f', '--fission_energy_release',
help='HDF5 file containing fission energy release data')
parser.add_argument('--libver', choices=['earliest', 'latest'],
default='earliest', help="Output HDF5 versioning. Use "
"'earliest' for backwards compatibility or 'latest' for "
@ -64,13 +62,6 @@ for basename, xs_list in sorted(suffixes.items()):
else:
data = openmc.data.IncidentNeutron.from_ace(filename, 'mcnp')
# Add fission energy release data, if available
if args.fission_energy_release is not None:
fer = openmc.data.FissionEnergyRelease.from_compact_hdf5(
args.fission_energy_release, data)
if fer is not None:
data.fission_energy = fer
# For each higher temperature, add cross sections to the existing table
for xs in xs_list[1:]:
filename = '.'.join((basename, xs))

8
src/api.F90
View file

@ -29,6 +29,10 @@ module openmc_api
use timer_header
use volume_calc, only: openmc_calculate_volumes
#ifdef DAGMC
use dagmc_header, only: free_memory_dagmc
#endif
implicit none
private
@ -149,6 +153,7 @@ contains
root_universe = -1
run_CE = .true.
run_mode = -1
dagmc = .false.
satisfy_triggers = .false.
call openmc_set_seed(DEFAULT_SEED)
source_latest = .false.
@ -325,6 +330,9 @@ contains
call free_memory_tally_filter()
call free_memory_tally_derivative()
call free_memory_bank()
#ifdef DAGMC
call free_memory_dagmc()
#endif
! Deallocate CMFD
call deallocate_cmfd(cmfd)

80
src/cell.cpp
View file

@ -15,7 +15,6 @@
#include "openmc/surface.h"
#include "openmc/xml_interface.h"
namespace openmc {
//==============================================================================
@ -209,7 +208,9 @@ Universe::to_hdf5(hid_t universes_group) const
// Cell implementation
//==============================================================================
Cell::Cell(pugi::xml_node cell_node)
CSGCell::CSGCell() {} // empty constructor
CSGCell::CSGCell(pugi::xml_node cell_node)
{
if (check_for_node(cell_node, "id")) {
id_ = std::stoi(get_node_value(cell_node, "id"));
@ -393,7 +394,7 @@ Cell::Cell(pugi::xml_node cell_node)
//==============================================================================
bool
Cell::contains(Position r, Direction u, int32_t on_surface) const
CSGCell::contains(Position r, Direction u, int32_t on_surface) const
{
if (simple_) {
return contains_simple(r, u, on_surface);
@ -405,7 +406,7 @@ Cell::contains(Position r, Direction u, int32_t on_surface) const
//==============================================================================
std::pair<double, int32_t>
Cell::distance(Position r, Direction u, int32_t on_surface) const
CSGCell::distance(Position r, Direction u, int32_t on_surface) const
{
double min_dist {INFTY};
int32_t i_surf {std::numeric_limits<int32_t>::max()};
@ -434,12 +435,12 @@ Cell::distance(Position r, Direction u, int32_t on_surface) const
//==============================================================================
void
Cell::to_hdf5(hid_t cells_group) const
CSGCell::to_hdf5(hid_t cell_group) const
{
// Create a group for this cell.
std::stringstream group_name;
group_name << "cell " << id_;
auto group = create_group(cells_group, group_name);
auto group = create_group(cell_group, group_name);
if (!name_.empty()) {
write_string(group, "name", name_, false);
@ -513,7 +514,7 @@ Cell::to_hdf5(hid_t cells_group) const
//==============================================================================
bool
Cell::contains_simple(Position r, Direction u, int32_t on_surface) const
CSGCell::contains_simple(Position r, Direction u, int32_t on_surface) const
{
for (int32_t token : rpn_) {
if (token < OP_UNION) {
@ -537,7 +538,7 @@ Cell::contains_simple(Position r, Direction u, int32_t on_surface) const
//==============================================================================
bool
Cell::contains_complex(Position r, Direction u, int32_t on_surface) const
CSGCell::contains_complex(Position r, Direction u, int32_t on_surface) const
{
// Make a stack of booleans. We don't know how big it needs to be, but we do
// know that rpn.size() is an upper-bound.
@ -585,6 +586,50 @@ Cell::contains_complex(Position r, Direction u, int32_t on_surface) const
}
}
//==============================================================================
// DAGMC Cell implementation
//==============================================================================
#ifdef DAGMC
DAGCell::DAGCell() : Cell{} {};
std::pair<double, int32_t>
DAGCell::distance(Position r, Direction u, int32_t on_surface) const
{
moab::ErrorCode rval;
moab::EntityHandle vol = dagmc_ptr_->entity_by_id(3, id_);
moab::EntityHandle hit_surf;
double dist;
double pnt[3] = {r.x, r.y, r.z};
double dir[3] = {u.x, u.y, u.z};
rval = dagmc_ptr_->ray_fire(vol, pnt, dir, hit_surf, dist);
MB_CHK_ERR_CONT(rval);
int surf_idx;
if (hit_surf != 0) {
surf_idx = dagmc_ptr_->index_by_handle(hit_surf);
} else { // indicate that particle is lost
surf_idx = -1;
}
return {dist, surf_idx};
}
bool DAGCell::contains(Position r, Direction u, int32_t on_surface) const
{
moab::ErrorCode rval;
moab::EntityHandle vol = dagmc_ptr_->entity_by_id(3, id_);
int result = 0;
double pnt[3] = {r.x, r.y, r.z};
double dir[3] = {u.x, u.y, u.z};
rval = dagmc_ptr_->point_in_volume(vol, pnt, result, dir);
MB_CHK_ERR_CONT(rval);
return result;
}
void DAGCell::to_hdf5(hid_t group_id) const { return; }
#endif
//==============================================================================
// Non-method functions
//==============================================================================
@ -601,7 +646,7 @@ read_cells(pugi::xml_node* node)
// Loop over XML cell elements and populate the array.
cells.reserve(n_cells);
for (pugi::xml_node cell_node: node->children("cell")) {
cells.push_back(new Cell(cell_node));
cells.push_back(new CSGCell(cell_node));
}
// Populate the Universe vector and map.
@ -725,6 +770,21 @@ extern "C" {
int cell_type(Cell* c) {return c->type_;}
#ifdef DAGMC
int32_t next_cell(DAGCell* cur_cell, DAGSurface* surf_xed )
{
moab::EntityHandle surf = surf_xed->dagmc_ptr_->entity_by_id(2,surf_xed->id_);
moab::EntityHandle vol = cur_cell->dagmc_ptr_->entity_by_id(3,cur_cell->id_);
moab::EntityHandle new_vol;
cur_cell->dagmc_ptr_->next_vol(surf, vol, new_vol);
return cur_cell->dagmc_ptr_->index_by_handle(new_vol);
}
#endif
int32_t cell_universe(Cell* c) {return c->universe_;}
int32_t cell_fill(Cell* c) {return c->fill_;}
@ -753,7 +813,7 @@ extern "C" {
{
cells.reserve(cells.size() + n);
for (int32_t i = 0; i < n; i++) {
cells.push_back(new Cell());
cells.push_back(new CSGCell());
}
n_cells = cells.size();
}

144
src/dagmc.cpp Normal file
View file

@ -0,0 +1,144 @@
#include "openmc/dagmc.h"
#include "openmc/error.h"
#include "openmc/string_functions.h"
#include "openmc/settings.h"
#include "openmc/geometry.h"
#include <string>
#include <sstream>
#include <algorithm>
#ifdef DAGMC
namespace openmc {
moab::DagMC* DAG;
void load_dagmc_geometry()
{
if (!DAG) {
DAG = new moab::DagMC();
}
int32_t dagmc_univ_id = 0; // universe is always 0 for DAGMC
moab::ErrorCode rval = DAG->load_file("dagmc.h5m");
MB_CHK_ERR_CONT(rval);
rval = DAG->init_OBBTree();
MB_CHK_ERR_CONT(rval);
std::vector<std::string> prop_keywords;
prop_keywords.push_back("mat");
prop_keywords.push_back("boundary");
std::map<std::string, std::string> ph;
DAG->parse_properties(prop_keywords, ph, ":");
MB_CHK_ERR_CONT(rval);
// initialize cell objects
n_cells = DAG->num_entities(3);
// Allocate the cell overlap count if necessary.
if (settings::check_overlaps) overlap_check_count.resize(n_cells, 0);
for (int i = 0; i < n_cells; i++) {
moab::EntityHandle vol_handle = DAG->entity_by_index(3, i+1);
// set cell ids using global IDs
DAGCell* c = new DAGCell();
c->id_ = DAG->id_by_index(3, i+1);
c->dagmc_ptr_ = DAG;
c->universe_ = dagmc_univ_id; // set to zero for now
c->fill_ = C_NONE; // no fill, single universe
cells.push_back(c);
cell_map[c->id_] = c->id_;
// Populate the Universe vector and dict
auto it = universe_map.find(dagmc_univ_id);
if (it == universe_map.end()) {
universes.push_back(new Universe());
universes.back()-> id_ = dagmc_univ_id;
universes.back()->cells_.push_back(i);
universe_map[dagmc_univ_id] = universes.size() - 1;
} else {
universes[it->second]->cells_.push_back(i);
}
if (DAG->is_implicit_complement(vol_handle)) {
// assuming implicit complement is void for now
c->material_.push_back(MATERIAL_VOID);
continue;
}
if (DAG->has_prop(vol_handle, "mat")){
std::string mat_value;
rval = DAG->prop_value(vol_handle, "mat", mat_value);
MB_CHK_ERR_CONT(rval);
to_lower(mat_value);
if (mat_value == "void" || mat_value == "vacuum") {
c->material_.push_back(MATERIAL_VOID);
} else {
c->material_.push_back(std::stoi(mat_value));
}
} else {
std::stringstream err_msg;
err_msg << "Volume " << c->id_ << " has no material assignment.";
fatal_error(err_msg.str());
}
}
// initialize surface objects
n_surfaces = DAG->num_entities(2);
surfaces.resize(n_surfaces);
for (int i = 0; i < n_surfaces; i++) {
moab::EntityHandle surf_handle = DAG->entity_by_index(2, i+1);
// set cell ids using global IDs
DAGSurface* s = new DAGSurface();
s->id_ = DAG->id_by_index(2, i+1);
s->dagmc_ptr_ = DAG;
if (DAG->has_prop(surf_handle, "boundary")) {
std::string bc_value;
rval = DAG->prop_value(surf_handle, "boundary", bc_value);
MB_CHK_ERR_CONT(rval);
to_lower(bc_value);
if (bc_value == "transmit" || bc_value == "transmission") {
s->bc_ = BC_TRANSMIT;
} else if (bc_value == "vacuum") {
s->bc_ = BC_VACUUM;
} else if (bc_value == "reflective" || bc_value == "reflect" || bc_value == "reflecting") {
s->bc_ = BC_REFLECT;
} else if (bc_value == "periodic") {
fatal_error("Periodic boundary condition not supported in DAGMC.");
} else {
std::stringstream err_msg;
err_msg << "Unknown boundary condition \"" << s->bc_
<< "\" specified on surface " << s->id_;
fatal_error(err_msg);
}
} else { // if no BC property is found, set to transmit
s->bc_ = BC_TRANSMIT;
}
// add to global array and map
surfaces[i] = s;
surface_map[s->id_] = s->id_;
}
return;
}
void free_memory_dagmc()
{
delete DAG;
}
}
#endif

20
src/dagmc_header.F90 Normal file
View file

@ -0,0 +1,20 @@
#ifdef DAGMC
module dagmc_header
use, intrinsic :: ISO_C_BINDING
implicit none
interface
subroutine load_dagmc_geometry() bind(C)
end subroutine load_dagmc_geometry
subroutine free_memory_dagmc() bind(C)
end subroutine free_memory_dagmc
end interface
end module dagmc_header
#endif

24
src/geometry.F90
View file

@ -55,10 +55,34 @@ module geometry
subroutine neighbor_lists() bind(C)
end subroutine neighbor_lists
#ifdef DAGMC
function next_cell_c(current_cell, surface_crossed) &
bind(C, name="next_cell") result(new_cell)
import C_PTR, C_INT32_T
type(C_PTR), intent(in), value :: current_cell
type(C_PTR), intent(in), value :: surface_crossed
integer(C_INT32_T) :: new_cell
end function next_cell_c
#endif
end interface
contains
#ifdef DAGMC
function next_cell(c, s) result(new_cell)
type(Cell), intent(in) :: c
type(Surface), intent(in) :: s
integer :: new_cell
new_cell = next_cell_c(c%ptr, s%ptr)
end function next_cell
#endif
!===============================================================================
! FIND_CELL determines what cell a source particle is in within a particular
! universe. If the base universe is passed, the particle should be found as long

108
src/input_xml.F90
View file

@ -11,6 +11,9 @@ module input_xml
use error, only: fatal_error, warning, write_message, openmc_err_msg
use geometry, only: neighbor_lists
use geometry_header
#ifdef DAGMC
use dagmc_header
#endif
use hdf5_interface
use list_header, only: ListChar, ListInt, ListReal
use material_header
@ -176,7 +179,9 @@ contains
! After reading input and basic geometry setup is complete, build lists of
! neighboring cells for efficient tracking
call neighbor_lists()
if (.not. dagmc) then
call neighbor_lists()
end if
! Assign temperatures to cells that don't have temperatures already assigned
call assign_temperatures()
@ -347,6 +352,63 @@ contains
end subroutine read_settings_xml_f
#ifdef DAGMC
!===============================================================================
! READ_GEOMETRY_DAGMC reads data from a DAGMC .h5m file, checking
! for material properties and surface boundary conditions
! some universe information is spoofed for now
!===============================================================================
subroutine read_geometry_dagmc()
integer :: i, j
integer :: univ_id
integer :: n_cells_in_univ
logical :: file_exists
character(MAX_LINE_LEN) :: filename
type(Cell), pointer :: c
type(VectorInt) :: univ_ids ! List of all universe IDs
type(DictIntInt) :: cells_in_univ_dict ! Used to count how many cells each
! universe contains
! Check if dagmc.h5m exists
filename = trim(path_input) // "dagmc.h5m"
inquire(FILE=filename, EXIST=file_exists)
if (.not. file_exists) then
call fatal_error("Geometry DAGMC file '" // trim(filename) // "' does not &
&exist!")
end if
call write_message("Reading DAGMC geometry...", 5)
call load_dagmc_geometry()
call allocate_surfaces()
call allocate_cells()
! setup universe data structs
do i = 1, n_cells
c => cells(i)
! additional metadata spoofing
univ_id = c % universe()
if (.not. cells_in_univ_dict % has(univ_id)) then
n_universes = n_universes + 1
n_cells_in_univ = 1
call universe_dict % set(univ_id, n_universes)
call univ_ids % push_back(univ_id)
else
n_cells_in_univ = 1 + cells_in_univ_dict % get(univ_id)
end if
call cells_in_univ_dict % set(univ_id, n_cells_in_univ)
end do
root_universe = find_root_universe()
end subroutine read_geometry_dagmc
#endif
!===============================================================================
! READ_GEOMETRY_XML reads data from a geometry.xml file and parses it, checking
! for errors and placing properly-formatted data in the right data structures
@ -374,6 +436,12 @@ contains
type(VectorInt) :: univ_ids ! List of all universe IDs
type(DictIntInt) :: cells_in_univ_dict ! Used to count how many cells each
! universe contains
#ifdef DAGMC
if (dagmc) then
call read_geometry_dagmc()
return
end if
#endif
! Display output message
call write_message("Reading geometry XML file...", 5)
@ -401,7 +469,6 @@ contains
! Allocate surfaces array
allocate(surfaces(n_surfaces))
do i = 1, n_surfaces
surfaces(i) % ptr = surface_pointer(i - 1);
@ -559,6 +626,40 @@ contains
end subroutine read_geometry_xml
subroutine allocate_surfaces()
integer :: i
! Allocate surfaces array
allocate(surfaces(n_surfaces))
do i = 1, n_surfaces
surfaces(i) % ptr = surface_pointer(i - 1);
! Add surface to dictionary
call surface_dict % set(surfaces(i) % id(), i)
end do
end subroutine allocate_surfaces
subroutine allocate_cells()
integer :: i
type(Cell), pointer :: c
! Allocate cells array
allocate(cells(n_cells))
do i = 1, n_cells
c => cells(i)
c % ptr = cell_pointer(i - 1)
! Check to make sure 'id' hasn't been used
if (cell_dict % has(c % id())) then
call fatal_error("Two or more cells use the same unique ID: " &
// to_str(c % id()))
end if
! Add cell to dictionary
call cell_dict % set(c % id(), i)
end do
end subroutine allocate_cells
!===============================================================================
! READ_MATERIAL_XML reads data from a materials.xml file and parses it, checking
! for errors and placing properly-formatted data in the right data structures
@ -2993,6 +3094,9 @@ contains
&OPENMC_MULTIPOLE_LIBRARY environment variable.")
end if
call already_read % clear()
call element_already_read % clear()
end subroutine read_ce_cross_sections
!===============================================================================

13
src/nuclide_header.F90
View file

@ -266,6 +266,19 @@ contains
subroutine nuclide_clear(this)
class(Nuclide), intent(inout) :: this ! The Nuclide object to clear
integer :: i
interface
subroutine reaction_delete(rx) bind(C)
import C_PTR
type(C_PTR), value :: rx
end subroutine reaction_delete
end interface
do i = 1, size(this % reactions)
call reaction_delete(this % reactions(i) % ptr)
end do
deallocate(this % reactions)
if (associated(this % multipole)) deallocate(this % multipole)

2
src/relaxng/settings.rnc
View file

@ -57,6 +57,8 @@ element settings {
element ptables { xsd:boolean }? &
element dagmc { xsd:boolean }? &
element run_cmfd { xsd:boolean }? &
element run_mode { xsd:string }? &

5
src/relaxng/settings.rng
View file

@ -254,6 +254,11 @@
<data type="boolean"/>
</element>
</optional>
<optional>
<element name="dagmc">
<data type="boolean"/>
</element>
</optional>
<optional>
<element name="run_cmfd">
<data type="boolean"/>

6
src/secondary_correlated.cpp
View file

@ -124,6 +124,12 @@ CorrelatedAngleEnergy::CorrelatedAngleEnergy(hid_t group)
m = mu.shape()[1] - offset_mu;
}
// For incoherent inelastic thermal scattering, the angle distributions
// may be given as discrete mu values. In this case, interpolation values
// of zero appear in the HDF5 file. Here we change it to a 1 so that
// int2interp doesn't fail.
if (interp_mu == 0) interp_mu = 1;
auto interp = int2interp(interp_mu);
auto xs = xt::view(mu, 0, xt::range(offset_mu, offset_mu + m));
auto ps = xt::view(mu, 1, xt::range(offset_mu, offset_mu + m));

3
src/settings.F90
View file

@ -80,6 +80,9 @@ module settings
! Mode to run in (fixed source, eigenvalue, plotting, etc)
integer(C_INT), bind(C) :: run_mode
! flag for use of DAGMC geometry
logical(C_BOOL), bind(C) :: dagmc
! Restart run
logical(C_BOOL), bind(C) :: restart_run

21
src/settings.cpp
View file

@ -56,7 +56,8 @@ bool ufs_on {false};
bool urr_ptables_on {true};
bool write_all_tracks {false};
bool write_initial_source {false};
bool dagmc {false};
std::string path_cross_sections;
std::string path_input;
std::string path_multipole;
@ -207,6 +208,17 @@ void read_settings_xml()
verbosity = std::stoi(get_node_value(root, "verbosity"));
}
// DAGMC geometry check
if (check_for_node(root, "dagmc")) {
dagmc = get_node_value_bool(root, "dagmc");
}
#ifndef DAGMC
if (dagmc) {
fatal_error("DAGMC mode unsupported for this build of OpenMC");
}
#endif
// To this point, we haven't displayed any output since we didn't know what
// the verbosity is. Now that we checked for it, show the title if necessary
if (openmc_master) {
@ -396,6 +408,13 @@ void read_settings_xml()
#endif
}
#ifdef _OPENMP
if (dagmc && omp_get_max_threads() > 1) {
warning("Forcing number of threads to 1 for DAGMC simulation.");
omp_set_num_threads(1);
}
#endif
// ==========================================================================
// EXTERNAL SOURCE

10
src/summary.F90
View file

@ -82,9 +82,9 @@ contains
integer(HID_T) :: nuclide_group
integer(HID_T) :: macro_group
integer :: i
character(12), allocatable :: nuc_names(:)
character(12), allocatable :: macro_names(:)
real(8), allocatable :: awrs(:)
character(kind=C_CHAR, len=20), allocatable :: nuc_names(:)
character(kind=C_CHAR, len=20), allocatable :: macro_names(:)
real(C_DOUBLE), allocatable :: awrs(:)
integer :: num_nuclides
integer :: num_macros
integer :: j
@ -172,8 +172,8 @@ contains
integer :: k
integer :: n
integer :: err
character(20), allocatable :: nuc_names(:)
character(20), allocatable :: macro_names(:)
character(kind=C_CHAR, len=20), allocatable :: nuc_names(:)
character(kind=C_CHAR, len=20), allocatable :: macro_names(:)
real(8) :: volume
real(8), allocatable :: nuc_densities(:)
integer :: num_nuclides

11
src/summary.cpp
View file

@ -2,13 +2,22 @@
#include "openmc/hdf5_interface.h"
#include "openmc/lattice.h"
#include "openmc/surface.h"
#include "openmc/settings.h"
namespace openmc {
extern "C" void
write_geometry(hid_t file_id) {
auto geom_group = create_group(file_id, "geometry");
#ifdef DAGMC
if (settings::dagmc) {
write_attribute(geom_group, "dagmc", 1);
return;
}
#endif
write_attribute(geom_group, "n_cells", cells.size());
write_attribute(geom_group, "n_surfaces", surfaces.size());
write_attribute(geom_group, "n_universes", universes.size());

75
src/surface.cpp
View file

@ -138,6 +138,8 @@ void read_coeffs(pugi::xml_node surf_node, int surf_id, double &c1, double &c2,
// Surface implementation
//==============================================================================
Surface::Surface() {} // empty constructor
Surface::Surface(pugi::xml_node surf_node)
{
if (check_for_node(surf_node, "id")) {
@ -207,8 +209,11 @@ Surface::reflect(Position r, Direction u) const
return u -= (2.0 * projection / magnitude) * n;
}
CSGSurface::CSGSurface() : Surface{} {};
CSGSurface::CSGSurface(pugi::xml_node surf_node) : Surface{surf_node} {};
void
Surface::to_hdf5(hid_t group_id) const
CSGSurface::to_hdf5(hid_t group_id) const
{
std::string group_name {"surface "};
group_name += std::to_string(id_);
@ -239,12 +244,62 @@ Surface::to_hdf5(hid_t group_id) const
close_group(surf_group);
}
//==============================================================================
// DAGSurface implementation
//==============================================================================
#ifdef DAGMC
DAGSurface::DAGSurface() : Surface{} {} // empty constructor
double DAGSurface::evaluate(Position r) const
{
return 0.0;
}
double
DAGSurface::distance(Position r, Direction u, bool coincident) const
{
moab::ErrorCode rval;
moab::EntityHandle surf = dagmc_ptr_->entity_by_id(2, id_);
moab::EntityHandle hit_surf;
double dist;
double pnt[3] = {r.x, r.y, r.z};
double dir[3] = {u.x, u.y, u.z};
rval = dagmc_ptr_->ray_fire(surf, pnt, dir, hit_surf, dist, NULL, 0, 0);
MB_CHK_ERR_CONT(rval);
if (dist < 0.0) dist = INFTY;
return dist;
}
Direction DAGSurface::normal(Position r) const
{
moab::ErrorCode rval;
Direction u;
moab::EntityHandle surf = dagmc_ptr_->entity_by_id(2, id_);
double pnt[3] = {r.x, r.y, r.z};
double dir[3] = {u.x, u.y, u.z};
rval = dagmc_ptr_->get_angle(surf, pnt, dir);
MB_CHK_ERR_CONT(rval);
return u;
}
BoundingBox DAGSurface::bounding_box() const
{
moab::ErrorCode rval;
moab::EntityHandle surf = dagmc_ptr_->entity_by_id(2, id_);
double min[3], max[3];
rval = dagmc_ptr_->getobb(surf, min, max);
MB_CHK_ERR_CONT(rval);
return {min[0], max[0], min[1], max[1], min[2], max[2]};
}
void DAGSurface::to_hdf5(hid_t group_id) const {}
#endif
//==============================================================================
// PeriodicSurface implementation
//==============================================================================
PeriodicSurface::PeriodicSurface(pugi::xml_node surf_node)
: Surface {surf_node}
: CSGSurface {surf_node}
{
if (check_for_node(surf_node, "periodic_surface_id")) {
i_periodic_ = std::stoi(get_node_value(surf_node, "periodic_surface_id"));
@ -580,7 +635,7 @@ axis_aligned_cylinder_normal(Position r, double offset1, double offset2)
//==============================================================================
SurfaceXCylinder::SurfaceXCylinder(pugi::xml_node surf_node)
: Surface(surf_node)
: CSGSurface(surf_node)
{
read_coeffs(surf_node, id_, y0_, z0_, radius_);
}
@ -614,7 +669,7 @@ void SurfaceXCylinder::to_hdf5_inner(hid_t group_id) const
//==============================================================================
SurfaceYCylinder::SurfaceYCylinder(pugi::xml_node surf_node)
: Surface(surf_node)
: CSGSurface(surf_node)
{
read_coeffs(surf_node, id_, x0_, z0_, radius_);
}
@ -647,7 +702,7 @@ void SurfaceYCylinder::to_hdf5_inner(hid_t group_id) const
//==============================================================================
SurfaceZCylinder::SurfaceZCylinder(pugi::xml_node surf_node)
: Surface(surf_node)
: CSGSurface(surf_node)
{
read_coeffs(surf_node, id_, x0_, y0_, radius_);
}
@ -680,7 +735,7 @@ void SurfaceZCylinder::to_hdf5_inner(hid_t group_id) const
//==============================================================================
SurfaceSphere::SurfaceSphere(pugi::xml_node surf_node)
: Surface(surf_node)
: CSGSurface(surf_node)
{
read_coeffs(surf_node, id_, x0_, y0_, z0_, radius_);
}
@ -833,7 +888,7 @@ axis_aligned_cone_normal(Position r, double offset1, double offset2,
//==============================================================================
SurfaceXCone::SurfaceXCone(pugi::xml_node surf_node)
: Surface(surf_node)
: CSGSurface(surf_node)
{
read_coeffs(surf_node, id_, x0_, y0_, z0_, radius_sq_);
}
@ -866,7 +921,7 @@ void SurfaceXCone::to_hdf5_inner(hid_t group_id) const
//==============================================================================
SurfaceYCone::SurfaceYCone(pugi::xml_node surf_node)
: Surface(surf_node)
: CSGSurface(surf_node)
{
read_coeffs(surf_node, id_, x0_, y0_, z0_, radius_sq_);
}
@ -899,7 +954,7 @@ void SurfaceYCone::to_hdf5_inner(hid_t group_id) const
//==============================================================================
SurfaceZCone::SurfaceZCone(pugi::xml_node surf_node)
: Surface(surf_node)
: CSGSurface(surf_node)
{
read_coeffs(surf_node, id_, x0_, y0_, z0_, radius_sq_);
}
@ -932,7 +987,7 @@ void SurfaceZCone::to_hdf5_inner(hid_t group_id) const
//==============================================================================
SurfaceQuadric::SurfaceQuadric(pugi::xml_node surf_node)
: Surface(surf_node)
: CSGSurface(surf_node)
{
read_coeffs(surf_node, id_, A_, B_, C_, D_, E_, F_, G_, H_, J_, K_);
}

2
src/tallies/tally_filter_universe.F90
View file

@ -78,7 +78,7 @@ contains
allocate(universe_ids(size(this % universes)))
do i = 1, size(this % universes)
universe_ids(i) = universe_id(this % universes(i)-1)
universe_ids(i) = universe_id(this % universes(i))
end do
call write_dataset(filter_group, "bins", universe_ids)
end subroutine to_statepoint_universe

16
src/thermal.cpp
View file

@ -487,12 +487,12 @@ ThermalData::sample(const NuclideMicroXS* micro_xs, double E,
// Determine endpoints on grid i
auto n = inelastic_data_[i].e_out.size();
double E_i_1 = inelastic_data_[i].e_out(0);
double E_i_J = inelastic_data_[i].e_out(n);
double E_i_J = inelastic_data_[i].e_out(n - 1);
// Determine endpoints on grid i + 1
n = inelastic_data_[i].e_out.size();
n = inelastic_data_[i + 1].e_out.size();
double E_i1_1 = inelastic_data_[i + 1].e_out(0);
double E_i1_J = inelastic_data_[i + 1].e_out(n);
double E_i1_J = inelastic_data_[i + 1].e_out(n - 1);
double E_1 = E_i_1 + f * (E_i1_1 - E_i_1);
double E_J = E_i_J + f * (E_i1_J - E_i_J);
@ -504,7 +504,7 @@ ThermalData::sample(const NuclideMicroXS* micro_xs, double E,
double c_j = inelastic_data_[l].e_out_cdf[0];
double c_j1;
std::size_t j;
for (j = 0; j < n - 2; ++j) {
for (j = 0; j < n - 1; ++j) {
c_j1 = inelastic_data_[l].e_out_cdf[j + 1];
if (r1 < c_j1) break;
c_j = c_j1;
@ -532,9 +532,9 @@ ThermalData::sample(const NuclideMicroXS* micro_xs, double E,
// Now interpolate between incident energy bins i and i + 1
if (l == i) {
*E_out = E_1 + (E - E_i_1) * (E_J - E_1) / (E_i_J - E_i_1);
*E_out = E_1 + (*E_out - E_i_1) * (E_J - E_1) / (E_i_J - E_i_1);
} else {
*E_out = E_1 + (E - E_i1_1) * (E_J - E_1) / (E_i1_J - E_i1_1);
*E_out = E_1 + (*E_out - E_i1_1) * (E_J - E_1) / (E_i1_J - E_i1_1);
}
// Sample outgoing cosine bin
@ -560,9 +560,9 @@ ThermalData::sample(const NuclideMicroXS* micro_xs, double E,
// Determine (k+1)th mu value
double mu_right;
if (k == n_inelastic_mu_ - 1) {
mu_right = 1.0 - *mu;
mu_right = 1.0;
} else {
mu_right = mu_l(j, k+1) + f*(mu_l(j+1, k+1) - mu_l(j, k+1)) - *mu;
mu_right = mu_l(j, k+1) + f*(mu_l(j+1, k+1) - mu_l(j, k+1));
}
// Smear angle

20
src/tracking.F90
View file

@ -7,6 +7,10 @@ module tracking
use geometry_header, only: cells
use geometry, only: find_cell, distance_to_boundary, cross_lattice,&
check_cell_overlap
#ifdef DAGMC
use geometry, only: next_cell
#endif
use material_header, only: materials, Material
use message_passing
use mgxs_interface
@ -309,6 +313,7 @@ contains
real(8) :: norm ! "norm" of surface normal
real(8) :: xyz(3) ! Saved global coordinate
integer :: i_surface ! index in surfaces
integer :: i_cell ! index of new cell
logical :: rotational ! if rotational periodic BC applied
logical :: found ! particle found in universe?
class(Surface), pointer :: surf
@ -467,6 +472,21 @@ contains
! ==========================================================================
! SEARCH NEIGHBOR LISTS FOR NEXT CELL
#ifdef DAGMC
if (dagmc) then
i_cell = next_cell(cells(p % last_cell(1) + 1), surfaces(abs(p % surface)))
! save material and temp
p % last_material = p % material
p % last_sqrtkT = p % sqrtKT
! set new cell value
p % coord(1) % cell = i_cell-1 ! decrement for C++ indexing
p % cell_instance = 1
p % material = cells(i_cell) % material(1)
p % sqrtKT = cells(i_cell) % sqrtKT(1)
return
end if
#endif
call find_cell(p, found, p % surface)
if (found) return

1
src/xml_interface.F90
View file

@ -12,6 +12,7 @@ module xml_interface
public :: check_for_node
public :: get_node_list
public :: get_node_value
public :: get_node_value_bool
public :: get_node_array
public :: node_value_string
public :: node_word_count

0
tests/regression_tests/dagmc/__init__.py Normal file
View file

BIN
tests/regression_tests/dagmc/dagmc.h5m Normal file
View file

Binary file not shown.

16
tests/regression_tests/dagmc/materials.xml Normal file
View file

@ -0,0 +1,16 @@
<?xml version="1.0"?>
<materials>
<material id="40">
<density value="11" units="g/cc" />
<nuclide name="U235" ao="1.0" />
</material>
<material id="41">
<density value="1.0" units="g/cc" />
<nuclide name="H1" ao="2.0" />
<nuclide name="O16" ao="1.0" />
<sab name="c_H_in_H2O"/>
</material>
</materials>

5
tests/regression_tests/dagmc/results_true.dat Normal file
View file

@ -0,0 +1,5 @@
k-combined:
1.115067E+00 5.423808E-02
tally 1:
8.543144E+00
1.530584E+01

16
tests/regression_tests/dagmc/settings.xml Normal file
View file

@ -0,0 +1,16 @@
<?xml version="1.0"?>
<settings>
<run_mode>eigenvalue</run_mode>
<dagmc>true</dagmc>
<batches>5</batches>
<inactive>0</inactive>
<particles>100</particles>
<!-- Starting source -->
<source>
<space type="box">
<parameters>-4 -4 -4 4 4 4</parameters>
</space>
</source>
</settings>

13
tests/regression_tests/dagmc/tallies.xml Normal file
View file

@ -0,0 +1,13 @@
<?xml version="1.0"?>
<tallies>
<filter id="1" type="cell">
<bins>1</bins>
</filter>
<tally id="1">
<filters>1</filters>
<scores>total </scores>
</tally>
</tallies>

12
tests/regression_tests/dagmc/test.py Normal file
View file

@ -0,0 +1,12 @@
from tests.testing_harness import TestHarness
import os
import pytest
import openmc
pytestmark = pytest.mark.skipif(
not openmc.capi.dagmc_enabled,
reason="DAGMC CAD geometry is not enabled.")
def test_dagmc():
harness = TestHarness('statepoint.5.h5')
harness.main()

19
tools/ci/download-xs.sh Executable file
View file

@ -0,0 +1,19 @@
#!/bin/bash
set -ex
# Download NNDC HDF5 data
if [[ ! -e $HOME/nndc_hdf5/cross_sections.xml ]]; then
wget -q -O - https://anl.box.com/shared/static/a0eflty17atnpd0pp7460exagr3nuhm7.xz | tar -C $HOME -xJ
fi
# Download ENDF/B-VII.1 distribution
ENDF=$HOME/endf-b-vii.1/
if [[ ! -d $ENDF/neutrons || ! -d $ENDF/photoat || ! -d $ENDF/atomic_relax ]]; then
wget -q -O - https://anl.box.com/shared/static/4kd2gxnf4gtk4w1c8eua5fsua22kvgjb.xz | tar -C $HOME -xJ
fi
# Download multipole library
if [[ ! -e $HOME/WMP_Library/092235.h5 ]]; then
wget -q https://github.com/mit-crpg/WMP_Library/releases/download/v1.0/WMP_Library_v1.0.tar.gz
tar -C $HOME -xzf WMP_Library_v1.0.tar.gz
fi

18
tools/ci/travis-before-script.sh
View file

@ -5,19 +5,5 @@ set -ex
# https://docs.travis-ci.com/user/gui-and-headless-browsers/#Using-xvfb-to-Run-Tests-That-Require-a-GUI
sh -e /etc/init.d/xvfb start
# Download NNDC HDF5 data
if [[ ! -e $HOME/nndc_hdf5/cross_sections.xml ]]; then
wget https://anl.box.com/shared/static/a0eflty17atnpd0pp7460exagr3nuhm7.xz -O - | tar -C $HOME -xvJ
fi
# Download ENDF/B-VII.1 distribution
ENDF=$HOME/endf-b-vii.1/
if [[ ! -d $ENDF/neutrons || ! -d $ENDF/photoat || ! -d $ENDF/atomic_relax ]]; then
wget https://anl.box.com/shared/static/4kd2gxnf4gtk4w1c8eua5fsua22kvgjb.xz -O - | tar -C $HOME -xvJ
fi
# Download multipole library
if [[ ! -e $HOME/WMP_Library/092235.h5 ]]; then
wget https://github.com/mit-crpg/WMP_Library/releases/download/v1.0/WMP_Library_v1.0.tar.gz
tar -C $HOME -xzvf WMP_Library_v1.0.tar.gz
fi
# Download NNDC HDF5 data, ENDF/B-VII.1 distribution, multipole library
source tools/ci/download-xs.sh

36
tools/ci/travis-install-dagmc.sh Executable file
View file

@ -0,0 +1,36 @@
#!/bin/bash
set -ex
# MOAB Variables
MOAB_BRANCH='Version5.0'
MOAB_REPO='https://bitbucket.org/fathomteam/moab/'
MOAB_INSTALL_DIR=$HOME/MOAB/
# DAGMC Variables
DAGMC_BRANCH='develop'
DAGMC_REPO='https://github.com/svalinn/dagmc'
DAGMC_INSTALL_DIR=$HOME/DAGMC/
CURRENT_DIR=$(pwd)
# MOAB Install
cd $HOME
mkdir MOAB && cd MOAB
git clone -b $MOAB_BRANCH $MOAB_REPO
mkdir build && cd build
cmake ../moab -DENABLE_HDF5=ON -DCMAKE_INSTALL_PREFIX=$MOAB_INSTALL_DIR
make -j && make -j test install
rm -rf $HOME/MOAB/moab
export LD_LIBRARY_PATH=$MOAB_INSTALL_DIR/lib:$LD_LIBRARY_PATH
# DAGMC Install
mkdir DAGMC && cd DAGMC
git clone -b $DAGMC_BRANCH $DAGMC_REPO
mkdir build && cd build
cmake ../dagmc -DBUILD_TALLY=ON -DCMAKE_INSTALL_PREFIX=$DAGMC_INSTALL_DIR
make -j install
rm -rf $HOME/DAGMC/dagmc
export LD_LIBRARY_PATH=$DAGMC_INSTALL_DIR/lib:$LD_LIBRARY_PATH
cd $CURRENT_DIR

14
tools/ci/travis-install.py
View file

@ -2,7 +2,6 @@ import os
import shutil
import subprocess
def which(program):
def is_exe(fpath):
return os.path.isfile(fpath) and os.access(fpath, os.X_OK)
@ -20,7 +19,7 @@ def which(program):
return None
def install(omp=False, mpi=False, phdf5=False):
def install(omp=False, mpi=False, phdf5=False, dagmc=False):
# Create build directory and change to it
shutil.rmtree('build', ignore_errors=True)
os.mkdir('build')
@ -48,23 +47,26 @@ def install(omp=False, mpi=False, phdf5=False):
else:
cmake_cmd.append('-DHDF5_PREFER_PARALLEL=OFF')
if dagmc:
cmake_cmd.append('-Ddagmc=ON')
# Build and install
cmake_cmd.append('..')
print(' '.join(cmake_cmd))
subprocess.check_call(cmake_cmd)
subprocess.check_call(['make', '-j'])
subprocess.check_call(['make', '-j4'])
subprocess.check_call(['sudo', 'make', 'install'])
def main():
# Convert Travis matrix environment variables into arguments for install()
omp = (os.environ.get('OMP') == 'y')
mpi = (os.environ.get('MPI') == 'y')
phdf5 = (os.environ.get('PHDF5') == 'y')
# Build and install
install(omp, mpi, phdf5)
dagmc = (os.environ.get('DAGMC') == 'y')
# Build and install
install(omp, mpi, phdf5, dagmc)
if __name__ == '__main__':
main()

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

@ -4,6 +4,11 @@ set -ex
# Install NJOY 2016
./tools/ci/travis-install-njoy.sh
# Install DAGMC if needed
if [[ $DAGMC = 'y' ]]; then
./tools/ci/travis-install-dagmc.sh
fi
# Upgrade pip before doing anything else
pip install --upgrade pip