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

Merge branch 'develop' into fix_coupled_tr_with_esr

Browse source
This commit is contained in:
Lorenzo Chierici 2026-07-21 09:53:47 +02:00 • committed by GitHub
commit 9fbf7aca06
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
239 changed files with 10120 additions and 1878 deletions
Download patch file Download diff file Expand all files Collapse all files

2
.github/workflows/ci.yml vendored
View file

@ -177,7 +177,7 @@ jobs:
- name: Setup tmate debug session
continue-on-error: true
if: ${{ contains(env.COMMIT_MESSAGE, '[gha-debug]') }}
if: ${{ failure() && contains(env.COMMIT_MESSAGE, '[gha-debug]') }}
uses: mxschmitt/action-tmate@v3
timeout-minutes: 10

1
.gitignore vendored
View file

@ -25,6 +25,7 @@ examples/**/*.xml
# Documentation builds
docs/build
docs/doxygen/xml
docs/source/_images/*.pdf
docs/source/_images/*.aux
docs/source/pythonapi/generated/

2
.readthedocs.yaml
View file

@ -7,6 +7,8 @@ build:
jobs:
post_checkout:
- git fetch --unshallow || true
- cd docs/doxygen && doxygen && cd -
sphinx:
configuration: docs/source/conf.py

26
CMakeLists.txt
View file

@ -532,7 +532,7 @@ else()
endif()
if(OPENMC_USE_DAGMC)
target_compile_definitions(libopenmc PRIVATE OPENMC_DAGMC_ENABLED)
target_compile_definitions(libopenmc PUBLIC OPENMC_DAGMC_ENABLED)
target_link_libraries(libopenmc dagmc-shared)
if(OPENMC_USE_UWUW)
@ -611,9 +611,7 @@ add_custom_command(TARGET libopenmc POST_BUILD
#===============================================================================
# Install executable, scripts, manpage, license
#===============================================================================
configure_file(cmake/OpenMCConfig.cmake.in "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/OpenMCConfig.cmake" @ONLY)
configure_file(cmake/OpenMCConfigVersion.cmake.in "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/OpenMCConfigVersion.cmake" @ONLY)
include(CMakePackageConfigHelpers)
set(INSTALL_CONFIGDIR ${CMAKE_INSTALL_LIBDIR}/cmake/OpenMC)
install(TARGETS openmc libopenmc
@ -627,10 +625,24 @@ install(EXPORT openmc-targets
NAMESPACE OpenMC::
DESTINATION ${INSTALL_CONFIGDIR})
configure_package_config_file(
"cmake/OpenMCConfig.cmake.in"
"${CMAKE_BINARY_DIR}/${CMAKE_FILES_DIRECTORY}/OpenMCConfig.cmake"
INSTALL_DESTINATION ${INSTALL_CONFIGDIR}
)
write_basic_package_version_file(
"${CMAKE_BINARY_DIR}/${CMAKE_FILES_DIRECTORY}/OpenMCConfigVersion.cmake"
VERSION ${OPENMC_VERSION}
COMPATIBILITY AnyNewerVersion
)
install(FILES
"${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/OpenMCConfig.cmake"
"${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/OpenMCConfigVersion.cmake"
DESTINATION ${INSTALL_CONFIGDIR})
"${CMAKE_BINARY_DIR}/${CMAKE_FILES_DIRECTORY}/OpenMCConfig.cmake"
"${CMAKE_BINARY_DIR}/${CMAKE_FILES_DIRECTORY}/OpenMCConfigVersion.cmake"
DESTINATION "${INSTALL_CONFIGDIR}"
)
install(FILES man/man1/openmc.1 DESTINATION ${CMAKE_INSTALL_MANDIR}/man1)
install(FILES LICENSE DESTINATION "${CMAKE_INSTALL_DOCDIR}" RENAME copyright)
install(DIRECTORY include/ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})

34
cmake/OpenMCConfig.cmake.in
View file

@ -1,12 +1,18 @@
get_filename_component(OpenMC_CMAKE_DIR "${CMAKE_CURRENT_LIST_FILE}" DIRECTORY)
@PACKAGE_INIT@
# Compute the install prefix from this file's location
get_filename_component(_OPENMC_PREFIX "${OpenMC_CMAKE_DIR}/../../.." ABSOLUTE)
include("${CMAKE_CURRENT_LIST_DIR}/OpenMCConfigVersion.cmake")
include(CMakeFindDependencyMacro)
# Explicitly calculate prefix if it was not generated above
if(NOT DEFINED PACKAGE_PREFIX_DIR)
get_filename_component(PACKAGE_PREFIX_DIR "${CMAKE_CURRENT_LIST_DIR}/../../.." ABSOLUTE)
endif()
find_dependency(fmt CONFIG REQUIRED HINTS ${PACKAGE_PREFIX_DIR})
find_dependency(pugixml CONFIG REQUIRED HINTS ${PACKAGE_PREFIX_DIR})
find_package(fmt CONFIG REQUIRED HINTS ${_OPENMC_PREFIX})
find_package(pugixml CONFIG REQUIRED HINTS ${_OPENMC_PREFIX})
if(@OPENMC_USE_DAGMC@)
find_package(DAGMC REQUIRED HINTS @DAGMC_DIR@)
find_dependency(DAGMC REQUIRED HINTS @DAGMC_DIR@)
endif()
if(@OPENMC_USE_LIBMESH@)
@ -16,20 +22,24 @@ if(@OPENMC_USE_LIBMESH@)
pkg_check_modules(LIBMESH REQUIRED @LIBMESH_PC_FILE@>=1.7.0 IMPORTED_TARGET)
endif()
find_package(PNG)
if(NOT TARGET OpenMC::libopenmc)
include("${OpenMC_CMAKE_DIR}/OpenMCTargets.cmake")
if("@PNG_FOUND@")
find_dependency(PNG)
endif()
if(@OPENMC_USE_MPI@)
find_package(MPI REQUIRED)
find_dependency(MPI REQUIRED)
endif()
if(@OPENMC_USE_OPENMP@)
find_package(OpenMP REQUIRED)
find_dependency(OpenMP REQUIRED)
endif()
if(@OPENMC_USE_UWUW@ AND NOT ${DAGMC_BUILD_UWUW})
message(FATAL_ERROR "UWUW is enabled in OpenMC but the DAGMC installation discovered was not configured with UWUW.")
endif()
include("${CMAKE_CURRENT_LIST_DIR}/OpenMCTargets.cmake")
if(NOT OpenMC_FIND_QUIETLY)
message(STATUS "Found OpenMC: ${PACKAGE_VERSION} (found in ${PACKAGE_PREFIX_DIR})")
endif()

11
cmake/OpenMCConfigVersion.cmake.in
View file

@ -1,11 +0,0 @@
set(PACKAGE_VERSION "@OPENMC_VERSION@")
# Check whether the requested PACKAGE_FIND_VERSION is compatible
if("${PACKAGE_VERSION}" VERSION_LESS "${PACKAGE_FIND_VERSION}")
set(PACKAGE_VERSION_COMPATIBLE FALSE)
else()
set(PACKAGE_VERSION_COMPATIBLE TRUE)
if ("${PACKAGE_VERSION}" VERSION_EQUAL "${PACKAGE_FIND_VERSION}")
set(PACKAGE_VERSION_EXACT TRUE)
endif()
endif()

1
docs/Makefile
View file

@ -45,6 +45,7 @@ help:
clean:
-rm -rf $(BUILDDIR)/*
-rm -rf source/pythonapi/generated/
-rm -rf doxygen/xml
html:
$(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html

13
docs/doxygen/Doxyfile Normal file
View file

@ -0,0 +1,13 @@
# Doxyfile 1.9.1
# This file describes the settings to be used by the documentation system
# doxygen (www.doxygen.org) for a project.
# Difference with default Doxyfile 1.9.1
PROJECT_NAME = OpenMC
QUIET = YES
WARN_IF_UNDOCUMENTED = NO
INPUT = ../../include/openmc/capi.h
GENERATE_HTML = NO
GENERATE_LATEX = NO
GENERATE_XML = YES

43
docs/source/capi/index.rst
View file

@ -46,43 +46,20 @@ Type Definitions
Functions
---------
.. c:function:: int openmc_calculate_volumes()
..
Once documentation is complete in capi.h, use:
.. doxygenfile:: capi.h
to populate this documentation without using
.. doxygenfunction::
for every function.
Run a stochastic volume calculation
.. doxygenfunction:: openmc_calculate_volumes
:return: Return status (negative if an error occurred)
:rtype: int
.. doxygenfunction:: openmc_cell_get_fill
.. c:function:: int openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n)
.. doxygenfunction:: openmc_cell_get_id
Get the fill for a cell
:param int32_t index: Index in the cells array
:param int* type: Type of the fill
:param int32_t** indices: Array of material indices for cell
:param int32_t* n: Length of indices array
:return: Return status (negative if an error occurred)
:rtype: int
.. c:function:: int openmc_cell_get_id(int32_t index, int32_t* id)
Get the ID of a cell
:param int32_t index: Index in the cells array
:param int32_t* id: ID of the cell
:return: Return status (negative if an error occurred)
:rtype: int
.. c:function:: int openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T)
Get the temperature of a cell
:param int32_t index: Index in the cells array
:param int32_t* instance: Which instance of the cell. If a null pointer is passed, the temperature
of the first instance is returned.
:param double* T: temperature of the cell
:return: Return status (negative if an error occurred)
:rtype: int
.. doxygenfunction:: openmc_cell_get_temperature
.. c:function:: int openmc_cell_get_density(int32_t index, const int32_t* instance, double* density)

13
docs/source/conf.py
View file

@ -11,7 +11,10 @@
# All configuration values have a default; values that are commented out
# serve to show the default.
import sys, os
import os
from pathlib import Path
import subprocess
import sys
# Determine if we're on Read the Docs server
on_rtd = os.environ.get('READTHEDOCS', None) == 'True'
@ -37,6 +40,7 @@ sys.path.insert(0, os.path.abspath('../..'))
# Add any Sphinx extension module names here, as strings. They can be extensions
# coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
extensions = [
'breathe',
'sphinx.ext.autodoc',
'sphinx.ext.napoleon',
'sphinx.ext.autosummary',
@ -47,6 +51,8 @@ extensions = [
]
if not on_rtd:
extensions.append('sphinxcontrib.rsvgconverter')
doxygen_dir = Path(__file__).parents[1] / 'doxygen'
subprocess.run(['doxygen'], cwd=doxygen_dir, check=True)
# Add any paths that contain templates here, relative to this directory.
templates_path = ['_templates']
@ -117,6 +123,11 @@ pygments_style = 'tango'
# A list of ignored prefixes for module index sorting.
#modindex_common_prefix = []
# -- Options breathe + doxygen -------------------------------------------------
breathe_projects = {"OpenMC": "../doxygen/xml"}
breathe_default_project = "OpenMC"
breathe_domain_by_file_pattern = {"*capi.h": "c"}
# -- Options for HTML output ---------------------------------------------------

3
docs/source/devguide/contributing.rst
View file

@ -111,7 +111,8 @@ The TC consists of the following individuals:
- `Paul Romano <https://github.com/paulromano>`_
- `Patrick Shriwise <https://github.com/pshriwise>`_
- `Adam Nelson <https://github.com/nelsonag>`_
- `Benoit Forget <https://github.com/bforget>`_
- `Jonathan Shimwell <https://github.com/shimwell>`_
- `John Tramm <https://github.com/jtramm>`_
The Project Lead is Paul Romano.

5
docs/source/devguide/docbuild.rst
View file

@ -14,6 +14,11 @@ Python API. That is, from the root directory of the OpenMC repository:
python -m pip install ".[docs]"
The OpenMC documentation also uses Doxygen to automatically generate its
C/C++ API documentation directly from the docstrings available in the source
code. You will need to have a working installation of Doxygen to generate the
documentation locally.
-----------------------------------
Building Documentation as a Webpage
-----------------------------------

4
docs/source/io_formats/collision_track.rst
View file

@ -10,7 +10,7 @@ may also be written after each batch when multiple files are requested
(``collision_track.N.h5``) or when the run is performed in parallel. The file
contains the information needed to reconstruct each recorded collision.
The current revision of the collision track file format is 1.1.
The current revision of the collision track file format is 1.2.
**/**
@ -33,7 +33,7 @@ The current revision of the collision track file format is 1.1.
- ``event_mt`` (*int*) -- ENDF MT number identifying the reaction.
- ``delayed_group`` (*int*) -- Delayed neutron group index (non-zero for delayed events).
- ``cell_id`` (*int*) -- ID of the cell in which the collision occurred.
- ``nuclide_id`` (*int*) -- ZA identifier of the nuclide (ZZZAAAM format).
- ``nuclide_id`` (*int*) -- PDG number of the nuclide (100ZZZAAAM).
- ``material_id`` (*int*) -- ID of the material containing the collision site.
- ``universe_id`` (*int*) -- ID of the universe containing the collision site.
- ``n_collision`` (*int*) -- Collision counter for the particle history.

61
docs/source/io_formats/geometry.rst
View file

@ -406,24 +406,55 @@ Each ``<dagmc_universe>`` element can have the following attributes or sub-eleme
*Default*: None
:material_overrides:
This element contains information on material overrides to be applied to the
DAGMC universe. It has the following attributes and sub-elements:
:cell:
Zero or more ``<cell>`` sub-elements may appear to override properties of
individual DAGMC volumes. Each ``<cell>`` element supports the following
attributes and sub-elements:
:cell:
Material override information for a single cell. It contains the following
attributes and sub-elements:
:id:
The integer cell ID in the DAGMC geometry to override. Required.
:id:
The cell ID in the DAGMC geometry for which the material override will
apply.
:name:
An optional string label for the cell.
:materials:
A list of material IDs that will apply to instances of the cell. If the
list contains only one ID, it will replace the original material
assignment of all instances of the DAGMC cell. If the list contains more
than one material, each material ID of the list will be assigned to the
various instances of the DAGMC cell.
*Default*: None
:material:
The material ID to assign to this cell. Use ``void`` for vacuum. Multiple
space-separated IDs may be given to specify a distribmat (distributed
material) assignment. Required.
:temperature:
Temperature(s) in [K] to assign to the cell. Must be greater than or equal
to 0. Multiple space-separated values may be given.
*Default*: None
:density:
Density in [g/cm³] to assign to the cell. Must be greater than 0. Requires a non-void
material fill. Multiple space-separated values may be given.
*Default*: None
:volume:
Volume of the cell in [cm³].
.. note:: DAGMC can compute cell volumes exactly from the triangulated
mesh surfaces. Specifying a manual volume risks inconsistency
with that capability.
*Default*: None
The following standard ``<cell>`` attributes are **not** supported inside
``<dagmc_universe>`` and will raise an error if present: ``region``,
``fill``, ``universe``, ``translation``, ``rotation``.
.. deprecated::
The ``<material_overrides>`` sub-element (containing ``<cell_override>``
children with ``<material_ids>``) is deprecated. A deprecation warning is
emitted and the overrides are converted to the ``<cell>`` format at parse
time. It is an error to specify both ``<material_overrides>`` and
``<cell>`` sub-elements on the same ``<dagmc_universe>``.
*Default*: None

209
docs/source/io_formats/settings.rst
View file

@ -98,6 +98,11 @@ sub-elements:
A list of strings representing the nuclide, to define specific
define specific target nuclide collisions to be banked.
.. note::
Electron and positron collision-track events are not associated with
a specific nuclide. If a ``nuclides`` entry is specified, these events
are omitted.
*Default*: None
:reactions:
@ -606,30 +611,30 @@ found in the :ref:`random ray user guide <random_ray>`.
*Default*: None
:adjoint_source:
Specifies an adjoint fixed source for adjoint transport simulations, and
follows the format for :ref:`source_element`. The distributions which make
up the adjoint source are subject to the same restrictions as forward
Specifies an adjoint fixed source for adjoint transport simulations, and
follows the format for :ref:`source_element`. The distributions which make
up the adjoint source are subject to the same restrictions as forward
fixed sources in Random Ray mode.
*Default*: None
:adjoint:
Specifies whether to perform adjoint transport. The default is 'False',
Specifies whether to perform adjoint transport. The default is 'False',
corresponding to forward transport.
*Default*: None
:volume_estimator:
Specifies choice of volume estimator for the random ray solver. Options
Specifies choice of volume estimator for the random ray solver. Options
are 'naive', 'simulation_averaged', or 'hybrid'. The default is 'hybrid'.
*Default*: None
:volume_normalized_flux_tallies:
Specifies whether to normalize flux tallies by volume (bool). The
default is 'False'. When enabled, flux tallies will be reported in units
of cm/cm^3. When disabled, flux tallies will be reported in units of cm
(i.e., total distance traveled by neutrons in the spatial tally
Specifies whether to normalize flux tallies by volume (bool). The
default is 'False'. When enabled, flux tallies will be reported in units
of cm/cm^3. When disabled, flux tallies will be reported in units of cm
(i.e., total distance traveled by neutrons in the spatial tally
region).
*Default*: None
@ -741,14 +746,16 @@ pseudo-random number generator.
*Default*: 1
--------------------
``<stride>`` Element
--------------------
-----------------------------------
``<shared_secondary_bank>`` Element
-----------------------------------
The ``stride`` element is used to specify how many random numbers are allocated
for each source particle history.
*Default*: 152,917
The ``shared_secondary_bank`` element indicates whether to use a shared
secondary particle bank. When enabled, secondary particles are collected into
a global bank, sorted for reproducibility, and load-balanced across MPI ranks
between generations. If not specified, the shared secondary bank is enabled
automatically for fixed-source simulations with weight windows active, and
disabled otherwise.
.. _source_element:
@ -770,9 +777,9 @@ attributes/sub-elements:
*Default*: 1.0
:type:
Indicator of source type. One of ``independent``, ``file``, ``compiled``, or
``mesh``. The type of the source will be determined by this attribute if it
is present.
Indicator of source type. One of ``independent``, ``file``, ``compiled``,
``mesh``, or ``tokamak``. The type of the source will be determined by this
attribute if it is present.
:particle:
The source particle type, specified as a PDG number or a string alias (e.g.,
@ -1008,6 +1015,80 @@ attributes/sub-elements:
mesh element and follows the format for :ref:`source_element`. The number of
``<source>`` sub-elements should correspond to the number of mesh elements.
For a source with ``type="tokamak"``, the spatial distribution is described by
a Miller-style flux-surface parameterization and the following sub-elements
are used instead of the ``space`` element:
:major_radius:
The major radius :math:`R_0` of the plasma in [cm].
:minor_radius:
The minor radius :math:`a` of the plasma in [cm]. Must be smaller than
``major_radius``.
:elongation:
The plasma elongation :math:`\kappa` (must be > 0).
:triangularity:
The plasma triangularity :math:`\delta` (must be in [-1, 1]). Negative
values describe negative-triangularity plasmas.
:shafranov_shift:
The Shafranov shift :math:`\Delta` in [cm] (must be >= 0 and less than
``minor_radius``/2).
:r_over_a:
A list of normalized minor-radius grid points :math:`r/a`. Must be strictly
increasing, start at 0, and end at 1.
:emission_density:
A list of neutron emission densities :math:`S(r)` evaluated at each
``r_over_a`` grid point (arbitrary units, must be non-negative). Only the
shape matters, since the profile is normalized internally. Values are
interpolated linearly between grid points and the profile is refined on an
internal grid for radial sampling. Must have the same length as
``r_over_a`` and contain at least one positive value.
:phi_start:
The starting toroidal angle in [rad].
*Default*: 0.0
:phi_extent:
The toroidal angle extent in [rad]. The source is sampled uniformly in
:math:`[\phi_\text{start},\ \phi_\text{start} + \phi_\text{extent}]`.
*Default*: :math:`2\pi`
:n_alpha:
The number of poloidal-angle grid points used to build the sampling CDFs
(must be > 2). Larger values reduce discretization bias; values below 51
produce a warning.
*Default*: 101
:vertical_shift:
A vertical shift of the plasma center in [cm].
*Default*: 0.0
:energy:
For a tokamak source, one or more ``energy`` sub-elements specify the
neutron energy distribution(s). Either a single distribution is given (used
at all radii) or exactly one distribution per ``r_over_a`` grid point is
given, in which case the energy is sampled from one of the two
distributions bracketing the sampled radius, selected stochastically with
probability proportional to the proximity of the radius to each grid point
(stochastic interpolation). Each follows the format of a univariate
probability distribution (see :ref:`univariate`).
:time:
An optional ``time`` sub-element specifying the time distribution of source
particles, following the format of a univariate probability distribution
(see :ref:`univariate`).
*Default*: particles are born at :math:`t=0`
.. note:: Biased sampling can be applied to the spatial and energy distributions
of a source by using the ``<bias>`` sub-element (see
:ref:`univariate` for details on how to specify bias distributions).
@ -1058,17 +1139,19 @@ variable and whose sub-elements/attributes are as follows:
:type:
The type of the distribution. Valid options are "uniform", "discrete",
"tabular", "maxwell", "watt", and "mixture". The "uniform" option produces
variates sampled from a uniform distribution over a finite interval. The
"discrete" option produces random variates that can assume a finite number
of values (i.e., a distribution characterized by a probability mass function).
The "tabular" option produces random variates sampled from a tabulated
distribution where the density function is either a histogram or
"tabular", "maxwell", "watt", "mixture", and "decay_spectrum". The "uniform"
option produces variates sampled from a uniform distribution over a finite
interval. The "discrete" option produces random variates that can assume a
finite number of values (i.e., a distribution characterized by a probability
mass function). The "tabular" option produces random variates sampled from a
tabulated distribution where the density function is either a histogram or
linearly-interpolated between tabulated points. The "watt" option produces
random variates is sampled from a Watt fission spectrum (only used for
energies). The "maxwell" option produce variates sampled from a Maxwell
fission spectrum (only used for energies). The "mixture" option produces samples
from univariate sub-distributions with given probabilities.
fission spectrum (only used for energies). The "mixture" option produces
samples from univariate sub-distributions with given probabilities. The
"decay_spectrum" option produces photon energies sampled from decay photon
spectra in a depletion chain (only used for energies).
*Default*: None
@ -1086,6 +1169,10 @@ variable and whose sub-elements/attributes are as follows:
:math:`(x,p)` pairs defining the discrete/tabular distribution. All :math:`x`
points are given first followed by corresponding :math:`p` points.
For a "decay_spectrum" distribution, ``parameters`` gives the atom densities
in [atom/b-cm] for the nuclides listed in the ``nuclides`` element, in the
same order.
For a "watt" distribution, ``parameters`` should be given as two real numbers
:math:`a` and :math:`b` that parameterize the distribution :math:`p(x) dx = c
e^{-x/a} \sinh \sqrt{b \, x} dx`.
@ -1115,6 +1202,21 @@ variable and whose sub-elements/attributes are as follows:
This sub-element of a ``pair`` element provides information on the
corresponding univariate distribution.
:volume:
For a "decay_spectrum" distribution, this attribute specifies the source
region volume in cm\ :sup:`3`. It is used together with atom densities to
determine the absolute photon emission rate. When a source uses a
"decay_spectrum" energy distribution, the source strength is set from this
emission rate.
:nuclides:
For a "decay_spectrum" distribution, this element specifies a
whitespace-separated list of nuclide names contributing to the decay photon
source. The atom densities for these nuclides are given by the ``parameters``
element in the same order. Nuclides are resolved against the depletion chain,
and nuclides without decay photon spectra do not contribute to the
distribution.
:bias:
This optional element specifies a biased distribution for importance sampling.
For continuous distributions, the ``bias`` element should contain another
@ -1135,23 +1237,6 @@ based on constraints.
*Default*: 0.05
-------------------------
``<state_point>`` Element
-------------------------
The ``<state_point>`` element indicates at what batches a state point file
should be written. A state point file can be used to restart a run or to get
tally results at any batch. The default behavior when using this tag is to
write out the source bank in the state_point file. This behavior can be
customized by using the ``<source_point>`` element. This element has the
following attributes/sub-elements:
:batches:
A list of integers separated by spaces indicating at what batches a state
point file should be written.
*Default*: Last batch only
--------------------------
``<source_point>`` Element
--------------------------
@ -1201,6 +1286,32 @@ attributes/sub-elements:
*Default*: false
-------------------------
``<state_point>`` Element
-------------------------
The ``<state_point>`` element indicates at what batches a state point file
should be written. A state point file can be used to restart a run or to get
tally results at any batch. The default behavior when using this tag is to
write out the source bank in the state_point file. This behavior can be
customized by using the ``<source_point>`` element. This element has the
following attributes/sub-elements:
:batches:
A list of integers separated by spaces indicating at what batches a state
point file should be written.
*Default*: Last batch only
--------------------
``<stride>`` Element
--------------------
The ``stride`` element is used to specify how many random numbers are allocated
for each source particle history.
*Default*: 152,917
------------------------------
``<surf_source_read>`` Element
------------------------------
@ -1300,7 +1411,7 @@ The ``<surface_grazing_cutoff>`` element specifies the surface flux cosine cutof
``<surface_grazing_ratio>`` Element
-----------------------------------
The ``<surface_grazing_ratio>`` element specifies the surface flux cosine
The ``<surface_grazing_ratio>`` element specifies the surface flux cosine
substitution ratio.
*Default*: 0.5
@ -1725,11 +1836,11 @@ mesh-based weight windows.
The ratio of the lower to upper weight window bounds.
*Default*: 5.0
For FW-CADIS:
:targets:
A sequence of IDs corresponding to the tallies which cover phase
A sequence of IDs corresponding to the tallies which cover phase
space regions of interest for local variance reduction.
*Default*: None

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

@ -26,6 +26,7 @@ Simulation Settings
openmc.FileSource
openmc.CompiledSource
openmc.MeshSource
openmc.TokamakSource
openmc.SourceParticle
openmc.VolumeCalculation
openmc.Settings

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

@ -22,6 +22,7 @@ Univariate Probability Distributions
openmc.stats.Legendre
openmc.stats.Mixture
openmc.stats.Normal
openmc.stats.DecaySpectrum
.. autosummary::
:toctree: generated

24
docs/source/usersguide/decay_sources.rst
View file

@ -132,8 +132,9 @@ can be run::
r2s.run(timesteps, source_rates, bounding_boxes=bounding_boxes)
If not specified otherwise, a photon transport calculation is run at each time
in the depletion schedule. That means in the case above, we would see three
photon transport calculations. To specify specific times at which photon
in the depletion schedule for which a decay photon source exists. Times without
a decay photon source, such as the initial state of a model containing only
stable nuclides, are omitted. To specify particular times at which photon
transport calculations should be run, pass the ``photon_time_indices`` argument.
For example, if we wanted to run a photon transport calculation only on the last
time (after the 5 hour decay), we would run::
@ -141,6 +142,19 @@ time (after the 5 hour decay), we would run::
r2s.run(timesteps, source_rates, bounding_boxes=bounding_boxes,
photon_time_indices=[2])
To attribute photon tally results to their parent radionuclides, set
``by_parent_nuclide=True``. This automatically adds a
:class:`openmc.ParentNuclideFilter` to every photon tally that does not already
have one. The filter bins are the union of radionuclides contributing to the
prepared decay photon sources. The resulting bins can be used directly when
inspecting the tally results::
r2s.run(timesteps, source_rates, bounding_boxes=bounding_boxes,
photon_time_indices=[2], by_parent_nuclide=True)
photon_tally = r2s.results['photon_tallies'][2][0]
tally_by_parent = photon_tally.get_pandas_dataframe()
After an R2S calculation has been run, the :class:`~openmc.deplete.R2SManager`
instance will have a ``results`` dictionary that allows you to directly access
results from each of the steps. It will also write out all the output files into
@ -148,12 +162,13 @@ a directory that is named "r2s_<timestamp>/". The ``output_dir`` argument to the
:meth:`~openmc.deplete.R2SManager.run` method enables you to override the
default output directory name if desired.
The :meth:`~openmc.deplete.R2SManager.run` method actually runs three
The :meth:`~openmc.deplete.R2SManager.run` method actually runs four
lower-level methods under the hood::
r2s.step1_neutron_transport(...)
r2s.step2_activation(...)
r2s.step3_photon_transport(...)
r2s.step3_photon_source(...)
r2s.step4_photon_transport(...)
For users looking for more control over the calculation, these lower-level
methods can be used in lieu of the :meth:`openmc.deplete.R2SManager.run` method.
@ -255,4 +270,3 @@ relevant tallies. This can be done with the aid of the
# Apply time correction factors
tally = d1s.apply_time_correction(dose_tally, factors, time_index)

38
docs/source/usersguide/depletion.rst
View file

@ -451,7 +451,6 @@ to transfer xenon from one material to another, you'd use::
integrator.add_transfer_rate(mat1, ['Xe'], 0.1, destination_material=mat2)
External Source Rates
=====================
External source rates define a fixed mass feed or removal of nuclides to or
from a depletable material. Unlike transfer rates, which are proportional to
@ -517,3 +516,40 @@ transfers between materials via ``destination_material``. See
:ref:`methods_depletion` for the augmented-matrix formulation used when both
features are active.
Comparing to Other Codes
Comparing depletion results from OpenMC with those from another code, such as
MCNP or Serpent, requires more than constructing equivalent transport models.
At each depletion step, differences in the transport solution, nuclear data,
reaction rate normalization, and numerical integration can all affect the
result. Small differences can also accumulate over successive depletion steps.
For a meaningful comparison, align as many of the following inputs and methods
as possible:
- Geometry and material definitions and associated physical properties such as
temperature
- Neutron cross section library (e.g., ENDF/B-VIII.0)
- Treatment of thermal scattering and unresolved resonance probability tables
- Neutron reactions accounted for in the depletion chain
- Decay data in the depletion chain
- Isomeric branching ratios for reactions in the depletion chain
- Fission product yields in the depletion chain
- Fission product yield interpolation method
(``CoupledOperator(fission_yield_mode=...)``)
- Reaction rate normalization, including fission Q values
(``CoupledOperator(fission_q=...)``)
- Depletion integration method (``PredictorIntegrator``, ``CECMIntegrator``,
etc.) and time-step sizes
When comparing to codes that use ACE format cross sections, it is recommended to
directly convert the ACE files to HDF5 format using functionality from the
:mod:`openmc.data` module (see :ref:`create_xs_library`). Some of the
LANL-distributed ACE libraries used with MCNP have also been converted to HDF5
format and are available for download at https://openmc.org/data.
Even after these choices have been aligned, exact agreement should not be
expected. Codes may use different approximations or numerical methods that
cannot be configured identically. When investigating a discrepancy, first
compare transport results and one-group reaction rates at the initial time, then
compare changes over subsequent timesteps.

10
docs/source/usersguide/processing.rst
View file

@ -97,7 +97,15 @@ VTK Mesh File Generation
------------------------
VTK files of OpenMC meshes can be created using the
:meth:`openmc.Mesh.write_data_to_vtk` method. Data can be applied to the
:meth:`openmc.Mesh.write_data_to_vtk` method. This method supports several VTK
formats depending on the mesh type. Structured meshes
(:class:`~openmc.RegularMesh`, :class:`~openmc.RectilinearMesh`,
:class:`~openmc.CylindricalMesh`, and :class:`~openmc.SphericalMesh`) can be
exported to legacy VTK format (``.vtk``). The :class:`~openmc.UnstructuredMesh`
class supports VTK unstructured grid formats (``.vtu``) as well as an HDF5-based
format (``.vtkhdf``) that does not require the ``vtk`` module to write.
Data can be applied to the
elements of the resulting mesh from mesh filter objects. This data can be
provided either as a flat array or, in the case of structured meshes
(:class:`~openmc.RegularMesh`, :class:`~openmc.RectilinearMesh`,

13
docs/source/usersguide/random_ray.rst
View file

@ -1105,10 +1105,15 @@ external source is present in the problem. Simulation settings (e.g., number of
rays, batches, etc.) will be identical for both calculations. At the
conclusion of the run, all results (e.g., tallies, plots, etc.) will be
derived from the adjoint flux rather than the forward flux but are not labeled
any differently. The initial forward flux solution will not be stored or
available in the final statepoint file. Those wishing to do analysis requiring
both the forward and adjoint solutions will need to run two separate
simulations and load both statepoint files.
any differently. When an initial forward solve is performed (i.e., when no
user-specified adjoint source is present), its output files are also written to
disk with a ``forward`` infix, so they are not overwritten by the subsequent
adjoint solve. This applies to the statepoint, ``tallies.out``, and any voxel
plots, e.g., ``statepoint.forward.N.h5`` and ``tallies.forward.out``; the
adjoint solve keeps the usual file names. This allows analyses requiring both
the forward and adjoint solutions to be performed from a single run. When
generating FW-CADIS weight windows, no weight window file is written for the
forward solve, as only the final adjoint-derived weight windows are meaningful.
.. note::
Use of the automated

52
docs/source/usersguide/settings.rst
View file

@ -291,6 +291,53 @@ example, the following would generate a photon source::
For a full list of all classes related to statistical distributions, see
:ref:`pythonapi_stats`.
Tokamak Plasma Sources
----------------------
For fusion applications, the :class:`openmc.TokamakSource` class provides a
native parametric neutron source for tokamak plasmas. Rather than specifying
spatial, angular, and energy distributions separately, the source is defined by
the plasma geometry (using a `Miller-style flux-surface parameterization
<https://doi.org/10.1063/1.872666>`_) and a radial emission profile. Source
sites are sampled directly from the plasma volume without rejection.
The plasma shape is described by the major radius :math:`R_0`, minor radius
:math:`a`, elongation :math:`\kappa`, triangularity :math:`\delta`, and
Shafranov shift :math:`\Delta`. The neutron birth profile is given as an
emission density :math:`S(r/a)` tabulated on a normalized minor-radius grid that
runs from 0 (magnetic axis) to 1 (last closed flux surface); only the shape of
the profile matters, since it is normalized internally. The emission density is
linearly interpolated between the supplied points and refined internally for
radial sampling. For example::
import numpy as np
r_over_a = np.linspace(0.0, 1.0, 50)
emission = (1.0 - r_over_a**2)**2 # peaked on-axis profile
source = openmc.TokamakSource(
major_radius=620.0, # cm
minor_radius=200.0, # cm
elongation=1.8,
triangularity=0.45,
shafranov_shift=10.0, # cm
r_over_a=r_over_a,
emission_density=emission,
energy=openmc.stats.muir(e0=14.08e6, m_rat=5.0, kt=20000.0),
)
settings.source = source
The ``energy`` argument accepts either a single
:class:`~openmc.stats.Univariate` distribution applied at all radii, or a
sequence with one distribution per ``r_over_a`` grid point to model a
radially-varying neutron spectrum (energies are then sampled by stochastic
interpolation between the two distributions bracketing the sampled radius). A
time distribution can be given with the ``time`` argument; by default, particles
are born at :math:`t=0`. The toroidal extent can be restricted with
``phi_start`` and ``phi_extent`` to model a sector of the plasma, and
``vertical_shift`` translates the plasma center along the z-axis.
File-based Sources
------------------
@ -792,6 +839,11 @@ collision_track.h5 file at the end of the simulation. The file contains
300 recorded collisions that occurred in materials with IDs 1 or 2, involving
fission or (n,2n) reactions on the nuclides U-238 or O-16, within cells
with IDs 5 and 12.
.. note::
Electron and positron collision-track events are not associated with a
specific nuclide. If a ``nuclides`` entry is specified, these events are omitted.
The file can be read using :func:`openmc.read_collision_track_file`.
The example below shows how to extract the data from the collision_track
feature and displays the fields stored in the file:

8
include/openmc/bank.h
View file

@ -34,18 +34,24 @@ extern vector<vector<double>> ifp_fission_lifetime_bank;
extern vector<int64_t> progeny_per_particle;
extern SharedArray<SourceSite> shared_secondary_bank_read;
extern SharedArray<SourceSite> shared_secondary_bank_write;
} // namespace simulation
//==============================================================================
// Non-member functions
//==============================================================================
void sort_fission_bank();
void sort_bank(SharedArray<SourceSite>& bank, bool is_fission_bank);
void free_memory_bank();
void init_fission_bank(int64_t max);
int64_t synchronize_global_secondary_bank(
SharedArray<SourceSite>& shared_secondary_bank);
} // namespace openmc
#endif // OPENMC_BANK_H

34
include/openmc/capi.h
View file

@ -9,14 +9,41 @@
extern "C" {
#endif
//! Run a stochastic volume calculation
//
//! \return Status (negative if an error occurred)
int openmc_calculate_volumes();
int openmc_cell_filter_get_bins(
int32_t index, const int32_t** cells, int32_t* n);
//! Get the fill for a cell
//
//! \param index Index in the cells array
//! \param type Type of the fill
//! \param indices Array of material indices for cell
//! \param n Length of indices array
//! \return Status (negative if an error occurred)
int openmc_cell_get_fill(
int32_t index, int* type, int32_t** indices, int32_t* n);
//! Get the ID of a cell
//
//! \param index Index in the cells array
//! \param id ID of the cell
//! \return Status (negative if an error occurred)
int openmc_cell_get_id(int32_t index, int32_t* id);
//! Get the temperature of a cell
//
//! \param index Index in the cells array
//! \param instance Which instance of the cell. If a null pointer is
//! passed, the temperature of the first instance is returned.
//! \param T temperature of the cell
//!\return Status (negative if an error occurred)
int openmc_cell_get_temperature(
int32_t index, const int32_t* instance, double* T);
int openmc_cell_get_density(
int32_t index, const int32_t* instance, double* rho);
int openmc_cell_get_translation(int32_t index, double xyz[]);
@ -123,8 +150,13 @@ int openmc_new_filter(const char* type, int32_t* index);
int openmc_next_batch(int* status);
int openmc_nuclide_name(int index, const char** name);
int openmc_plot_geometry();
// Deprecated; use openmc_slice_data.
int openmc_id_map(const void* slice, int32_t* data_out);
// Deprecated; use openmc_slice_data.
int openmc_property_map(const void* slice, double* data_out);
int openmc_slice_data(const double origin[3], const double u_span[3],
const double v_span[3], const size_t pixels[2], bool show_overlaps, int level,
int32_t filter_index, int32_t* geom_data, double* property_data);
int openmc_get_plot_index(int32_t id, int32_t* index);
int openmc_plot_get_id(int32_t index, int32_t* id);
int openmc_plot_set_id(int32_t index, int32_t id);
@ -275,7 +307,7 @@ int openmc_zernike_filter_set_params(
int openmc_particle_filter_get_bins(int32_t idx, int32_t bins[]);
//! Sets the mesh and energy grid for CMFD reweight
//! \param[in] meshtyally_id id of CMFD Mesh Tally
//! \param[in] meshtally_id id of CMFD Mesh Tally
//! \param[in] cmfd_indices indices storing spatial and energy dimensions of
//! CMFD problem \param[in] norm CMFD normalization factor
void openmc_initialize_mesh_egrid(

24
include/openmc/cell.h
View file

@ -145,6 +145,30 @@ private:
bool simple_; //!< Does the region contain only intersections?
};
//==============================================================================
// XML parsing helpers for <cell> nodes
//==============================================================================
//! Parse material IDs from a <cell> XML node.
//! \param node XML node containing a "material" attribute or child element
//! \param cell_id Cell ID used in error messages
//! \return Vector of material IDs (MATERIAL_VOID for "void")
vector<int32_t> parse_cell_material_xml(pugi::xml_node node, int32_t cell_id);
//! Parse temperatures in [K] from a <cell> XML node.
//! Validates that all values are non-negative and the list is non-empty.
//! \param node XML node containing a "temperature" attribute or child element
//! \param cell_id Cell ID used in error messages
//! \return Vector of temperatures in [K]
vector<double> parse_cell_temperature_xml(pugi::xml_node node, int32_t cell_id);
//! Parse densities in [g/cm³] from a <cell> XML node.
//! Validates that all values are positive and the list is non-empty.
//! \param node XML node containing a "density" attribute or child element
//! \param cell_id Cell ID used in error messages
//! \return Vector of densities in [g/cm³]
vector<double> parse_cell_density_xml(pugi::xml_node node, int32_t cell_id);
//==============================================================================
class Cell {

2
include/openmc/chain.h
View file

@ -101,6 +101,8 @@ extern vector<unique_ptr<ChainNuclide>> chain_nuclides;
void read_chain_file_xml();
void free_memory_chain();
} // namespace openmc
#endif // OPENMC_CHAIN_H

3
include/openmc/constants.h
View file

@ -35,7 +35,7 @@ constexpr array<int, 2> VERSION_VOXEL {2, 0};
constexpr array<int, 2> VERSION_MGXS_LIBRARY {1, 0};
constexpr array<int, 2> VERSION_PROPERTIES {1, 1};
constexpr array<int, 2> VERSION_WEIGHT_WINDOWS {1, 0};
constexpr array<int, 2> VERSION_COLLISION_TRACK {1, 1};
constexpr array<int, 2> VERSION_COLLISION_TRACK {1, 2};
// ============================================================================
// ADJUSTABLE PARAMETERS
@ -368,6 +368,7 @@ enum class SolverType { MONTE_CARLO, RANDOM_RAY };
enum class RandomRayVolumeEstimator { NAIVE, SIMULATION_AVERAGED, HYBRID };
enum class RandomRaySourceShape { FLAT, LINEAR, LINEAR_XY };
enum class RandomRaySampleMethod { PRNG, HALTON, S2 };
enum class RandomRaySolve { FORWARD, FORWARD_FOR_ADJOINT, ADJOINT };
//==============================================================================
// Geometry Constants

19
include/openmc/dagmc.h
View file

@ -94,6 +94,10 @@ private:
class DAGUniverse : public Universe {
public:
using MaterialOverrides = std::unordered_map<int32_t, vector<int32_t>>;
using TemperatureOverrides = std::unordered_map<int32_t, vector<double>>;
using DensityOverrides = std::unordered_map<int32_t, vector<double>>;
explicit DAGUniverse(pugi::xml_node node);
//! Create a new DAGMC universe
@ -112,6 +116,9 @@ public:
//! Initialize the DAGMC accel. data structures, indices, material
//! assignments, etc.
void initialize();
void initialize(const MaterialOverrides& material_overrides,
const TemperatureOverrides& temperature_overrides,
const DensityOverrides& density_overrides = {});
//! Reads UWUW materials and returns an ID map
void read_uwuw_materials();
@ -146,7 +153,8 @@ public:
//! Assign a material overriding normal assignement to a cell
//! \param[in] c The OpenMC cell to which the material is assigned
void override_assign_material(std::unique_ptr<DAGCell>& c) const;
void override_assign_material(std::unique_ptr<DAGCell>& c,
const MaterialOverrides& material_overrides) const;
//! Return the index into the model cells vector for a given DAGMC volume
//! handle in the universe
@ -187,7 +195,9 @@ private:
void set_id(); //!< Deduce the universe id from model::universes
void init_dagmc(); //!< Create and initialise DAGMC pointer
void init_metadata(); //!< Create and initialise dagmcMetaData pointer
void init_geometry(); //!< Create cells and surfaces from DAGMC entities
void init_geometry(const MaterialOverrides& material_overrides,
const TemperatureOverrides& temperature_overrides,
const DensityOverrides& density_overrides);
std::string
filename_; //!< Name of the DAGMC file used to create this universe
@ -201,11 +211,6 @@ private:
//!< generate new material IDs for the universe
bool has_graveyard_; //!< Indicates if the DAGMC geometry has a "graveyard"
//!< volume
std::unordered_map<int32_t, vector<int32_t>>
material_overrides_; //!< Map of material overrides
//!< keys correspond to the DAGMCCell id
//!< values are a list of material ids used
//!< for the override
};
//==============================================================================

65
include/openmc/distribution.h
View file

@ -407,6 +407,71 @@ private:
double integral_; //!< Integral of distribution
};
//==============================================================================
// DecaySpectrum — non-owning mixture of decay photon distributions
//==============================================================================
//! Energy distribution formed by mixing multiple decay photon spectra.
//!
//! Unlike the general Mixture distribution, this class holds non-owning
//! pointers to the component distributions (which live in
//! data::chain_nuclides). Each component is weighted by the activity
//! (atoms * decay_constant) of the corresponding nuclide.
class DecaySpectrum : public Distribution {
public:
//============================================================================
// Types, aliases
struct Sample {
double energy;
double weight;
int parent_nuclide;
};
//============================================================================
// Constructors
//! Construct from an XML node containing nuclide names and atom densities.
//!
//! Reads child ``<nuclide>`` elements with ``name`` and ``density``
//! attributes, resolves them against the loaded depletion chain, and
//! constructs the mixed distribution.
explicit DecaySpectrum(pugi::xml_node node);
//============================================================================
// Methods
//! Sample a value from the distribution and return the parent nuclide index
//! \param seed Pseudorandom number seed pointer
//! \return (Sampled energy, sample weight, chain nuclide index)
Sample sample_with_parent(uint64_t* seed) const;
//! Sample a value from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return (sampled value, sample weight)
std::pair<double, double> sample(uint64_t* seed) const override;
double integral() const override;
protected:
//! Sample a value (unbiased) from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample_unbiased(uint64_t* seed) const override;
private:
//! Initialize decay spectrum sampling data
//! \param nuclide_indices Indices of decay photon emitters in
//! data::chain_nuclides
//! \param atoms Number of atoms for each component.
void init(vector<int> nuclide_indices, const vector<double>& atoms);
vector<int> nuclide_indices_; //!< Indices of emitting nuclides in the chain
DiscreteIndex di_; //!< Discrete index for component selection
double integral_; //!< Total photon emission rate
};
} // namespace openmc
#endif // OPENMC_DISTRIBUTION_H

4
include/openmc/error.h
View file

@ -64,14 +64,14 @@ void write_message(
int level, const std::string& message, const Params&... fmt_args)
{
if (settings::verbosity >= level) {
write_message(fmt::format(message, fmt_args...));
write_message(fmt::format(fmt::runtime(message), fmt_args...));
}
}
template<typename... Params>
void write_message(const std::string& message, const Params&... fmt_args)
{
write_message(fmt::format(message, fmt_args...));
write_message(fmt::format(fmt::runtime(message), fmt_args...));
}
#ifdef OPENMC_MPI

13
include/openmc/event.h
View file

@ -112,6 +112,19 @@ void process_collision_events();
//! \param n_particles The number of particles in the particle buffer
void process_death_events(int64_t n_particles);
//! Process event queues until all are empty. Each iteration processes the
//! longest queue first to maximize vectorization efficiency.
void process_transport_events();
//! Initialize secondary particles from a shared secondary bank for
//! event-based transport
//
//! \param n_particles The number of particles to initialize
//! \param offset The offset index in the shared secondary bank
//! \param shared_secondary_bank The shared secondary bank to read from
void process_init_secondary_events(int64_t n_particles, int64_t offset,
const SharedArray<SourceSite>& shared_secondary_bank);
} // namespace openmc
#endif // OPENMC_EVENT_H

38
include/openmc/geometry.h
View file

@ -3,9 +3,12 @@
#include <cmath>
#include <cstdint>
#include <unordered_map>
#include <vector>
#include "openmc/array.h"
#include "openmc/constants.h"
#include "openmc/random_ray/source_region.h" // For hash_combine
#include "openmc/vector.h"
namespace openmc {
@ -13,6 +16,34 @@ namespace openmc {
class BoundaryInfo;
class GeometryState;
//==============================================================================
//! OverlapKey to store cell and universe data of a single overlap, along with
//! a functor for hashing an OverlapKey into an unordered_map.
//==============================================================================
struct OverlapKey {
int universe_id;
int cell1_id;
int cell2_id;
bool operator==(const OverlapKey& other) const
{
return universe_id == other.universe_id && cell1_id == other.cell1_id &&
cell2_id == other.cell2_id;
}
};
struct OverlapKeyHash {
std::size_t operator()(const OverlapKey& k) const
{
size_t seed = 0;
hash_combine(seed, k.universe_id);
hash_combine(seed, k.cell1_id);
hash_combine(seed, k.cell2_id);
return seed;
}
};
//==============================================================================
// Global variables
//==============================================================================
@ -24,6 +55,10 @@ extern "C" int n_coord_levels; //!< Number of CSG coordinate levels
extern vector<int64_t> overlap_check_count;
// Overlap data structures get cleared every slice_data run
extern vector<OverlapKey> overlap_keys;
extern std::unordered_map<OverlapKey, int, OverlapKeyHash> overlap_key_index;
} // namespace model
//==============================================================================
@ -38,8 +73,7 @@ inline bool coincident(double d1, double d2)
//==============================================================================
//! Check for overlapping cells at a particle's position.
//==============================================================================
bool check_cell_overlap(GeometryState& p, bool error = true);
int check_cell_overlap(GeometryState& p, bool error = true);
//==============================================================================
//! Get the cell instance for a particle at the specified universe level

26
include/openmc/math_functions.h
View file

@ -213,6 +213,20 @@ double exprel(double x);
//! \return log(1+x)/x without loss of precision near 0
double log1prel(double x);
//! Evaluate the cylindrical Bessel function of the first kind J_n(x)
//!
//! Uses std::cyl_bessel_j where available (e.g., libstdc++). On standard
//! library implementations lacking the C++17 special math functions (e.g.,
//! libc++ on Apple Clang/LLVM), falls back to the ascending power series,
//! which converges to machine precision for the small arguments (|x| <= 2)
//! used in OpenMC. Unlike std::cyl_bessel_j, negative arguments are handled
//! via the parity relation J_n(-x) = (-1)^n J_n(x).
//!
//! \param n Non-negative integer order of the Bessel function
//! \param x Real argument
//! \return J_n(x)
double cyl_bessel_j(int n, double x);
//! Helper function to get index and interpolation function on an incident
//! energy grid
//!
@ -233,5 +247,17 @@ void get_energy_index(
double standard_normal_cdf(double z);
//==============================================================================
//! Return true if two floating-point values are approximately equal within a
//! combined relative and absolute tolerance.
//!
//! \param a first floating point value
//! \param b second floating point value
//! \param rel_tol relative tolerance
//! \param abs_tol absolute tolerance
//! \return true if a and b are approximately equal, false otherwise
//==============================================================================
bool isclose(double a, double b, double rel_tol, double abs_tol);
} // namespace openmc
#endif // OPENMC_MATH_FUNCTIONS_H

3
include/openmc/nuclide.h
View file

@ -84,6 +84,9 @@ public:
double collapse_rate(int MT, double temperature, span<const double> energy,
span<const double> flux) const;
//! Return a ParticleType object representing this nuclide
ParticleType particle_type() const { return {Z_, A_, metastable_}; }
//============================================================================
// Data members
std::string name_; //!< Name of nuclide, e.g. "U235"

3
include/openmc/particle.h
View file

@ -71,7 +71,8 @@ public:
void event_advance();
void event_cross_surface();
void event_collide();
void event_revive_from_secondary();
void event_revive_from_secondary(const SourceSite& site);
void event_check_limit_and_revive();
void event_death();
//! pulse-height recording

40
include/openmc/particle_data.h
View file

@ -50,8 +50,11 @@ struct SourceSite {
// Extra attributes that don't show up in source written to file
int parent_nuclide {-1};
int64_t parent_id;
int64_t progeny_id;
int64_t parent_id {0};
int64_t progeny_id {0};
double wgt_born {1.0};
double wgt_ww_born {-1.0};
int64_t n_split {0};
};
struct CollisionTrackSite {
@ -533,14 +536,14 @@ private:
uint64_t seeds_[N_STREAMS];
int stream_;
vector<SourceSite> secondary_bank_;
vector<SourceSite> local_secondary_bank_;
// Keep track of how many secondary particles were created in the collision
// and what the starting index is in the secondary bank for this particle
int n_secondaries_ {0};
int secondary_bank_index_ {0};
int64_t current_work_;
int64_t current_work_ {0};
vector<double> flux_derivs_;
@ -563,7 +566,9 @@ private:
int n_event_ {0};
int n_split_ {0};
int64_t n_tracks_ {0}; //!< number of tracks in this particle history
int64_t n_split_ {0};
double ww_factor_ {0.0};
int64_t n_progeny_ {0};
@ -577,10 +582,8 @@ public:
// Methods and accessors
// Cross section caches
NuclideMicroXS& neutron_xs(int i)
{
return neutron_xs_[i];
} // Microscopic neutron cross sections
// Microscopic neutron cross sections
NuclideMicroXS& neutron_xs(int i) { return neutron_xs_[i]; }
const NuclideMicroXS& neutron_xs(int i) const { return neutron_xs_[i]; }
// Microscopic photon cross sections
@ -693,12 +696,14 @@ public:
int& stream() { return stream_; }
// secondary particle bank
SourceSite& secondary_bank(int i) { return secondary_bank_[i]; }
const SourceSite& secondary_bank(int i) const { return secondary_bank_[i]; }
decltype(secondary_bank_)& secondary_bank() { return secondary_bank_; }
decltype(secondary_bank_) const& secondary_bank() const
SourceSite& local_secondary_bank(int i) { return local_secondary_bank_[i]; }
const SourceSite& local_secondary_bank(int i) const
{
return secondary_bank_;
return local_secondary_bank_[i];
}
decltype(local_secondary_bank_)& local_secondary_bank()
{
return local_secondary_bank_;
}
// Number of secondaries created in a collision
@ -747,13 +752,16 @@ public:
int& n_event() { return n_event_; }
// Number of times variance reduction has caused a particle split
int n_split() const { return n_split_; }
int& n_split() { return n_split_; }
int64_t n_split() const { return n_split_; }
int64_t& n_split() { return n_split_; }
// Particle-specific factor for on-the-fly weight window adjustment
double ww_factor() const { return ww_factor_; }
double& ww_factor() { return ww_factor_; }
// Number of tracks in this particle history
int64_t& n_tracks() { return n_tracks_; }
// Number of progeny produced by this particle
int64_t& n_progeny() { return n_progeny_; }

124
include/openmc/plot.h
View file

@ -5,6 +5,7 @@
#include <sstream>
#include <unordered_map>
#include <unordered_set>
#include <vector>
#include "openmc/tensor.h"
#include "pugixml.hpp"
@ -18,6 +19,8 @@
#include "openmc/position.h"
#include "openmc/random_lcg.h"
#include "openmc/ray.h"
#include "openmc/tallies/filter.h"
#include "openmc/tallies/filter_match.h"
#include "openmc/xml_interface.h"
namespace openmc {
@ -148,11 +151,12 @@ public:
struct IdData {
// Constructor
IdData(size_t h_res, size_t v_res);
IdData(size_t h_res, size_t v_res, bool include_filter = false);
// Methods
void set_value(size_t y, size_t x, const GeometryState& p, int level);
void set_overlap(size_t y, size_t x);
void set_value(size_t y, size_t x, const Particle& p, int level,
Filter* filter = nullptr, FilterMatch* match = nullptr);
void set_overlap(size_t y, size_t x, int overlap_idx);
// Members
tensor::Tensor<int32_t> data_; //!< 2D array of cell & material ids
@ -160,16 +164,34 @@ struct IdData {
struct PropertyData {
// Constructor
PropertyData(size_t h_res, size_t v_res);
PropertyData(size_t h_res, size_t v_res, bool include_filter = false);
// Methods
void set_value(size_t y, size_t x, const GeometryState& p, int level);
void set_overlap(size_t y, size_t x);
void set_value(size_t y, size_t x, const Particle& p, int level,
Filter* filter = nullptr, FilterMatch* match = nullptr);
void set_overlap(size_t y, size_t x, int overlap_idx);
// Members
tensor::Tensor<double> data_; //!< 2D array of temperature & density data
};
struct RasterData {
// Constructor
RasterData(size_t h_res, size_t v_res, bool include_filter = false);
// Methods
void set_value(size_t y, size_t x, const Particle& p, int level,
Filter* filter = nullptr, FilterMatch* match = nullptr);
void set_overlap(size_t y, size_t x, int overlap_idx);
// Members
tensor::Tensor<int32_t>
id_data_; //!< [v_res, h_res, 3 or 4]: cell, instance, mat, [filter_bin]
tensor::Tensor<double>
property_data_; //!< [v_res, h_res, 2]: temperature, density
bool include_filter_; //!< Whether filter bin index is included
};
//===============================================================================
// Plot class
//===============================================================================
@ -177,7 +199,7 @@ struct PropertyData {
class SlicePlotBase {
public:
template<class T>
T get_map() const;
T get_map(int32_t filter_index = -1) const;
enum class PlotBasis { xy = 1, xz = 2, yz = 3 };
@ -188,70 +210,65 @@ public:
// Members
public:
Position origin_; //!< Plot origin in geometry
Position width_; //!< Plot width in geometry
PlotBasis basis_; //!< Plot basis (XY/XZ/YZ)
array<size_t, 3> pixels_; //!< Plot size in pixels
bool slice_color_overlaps_; //!< Show overlapping cells?
int slice_level_ {-1}; //!< Plot universe level
Position origin_; //!< Plot origin in geometry
Direction u_span_; //!< Full-width span vector in geometry
Direction v_span_; //!< Full-height span vector in geometry
array<size_t, 3> pixels_; //!< Plot size in pixels
bool show_overlaps_; //!< Show overlapping cells?
int slice_level_ {-1}; //!< Plot universe level
private:
};
template<class T>
T SlicePlotBase::get_map() const
T SlicePlotBase::get_map(int32_t filter_index) const
{
size_t width = pixels_[0];
size_t height = pixels_[1];
// get pixel size
double in_pixel = (width_[0]) / static_cast<double>(width);
double out_pixel = (width_[1]) / static_cast<double>(height);
// size data array
T data(width, height);
// setup basis indices and initial position centered on pixel
int in_i, out_i;
Position xyz = origin_;
switch (basis_) {
case PlotBasis::xy:
in_i = 0;
out_i = 1;
break;
case PlotBasis::xz:
in_i = 0;
out_i = 2;
break;
case PlotBasis::yz:
in_i = 1;
out_i = 2;
break;
default:
UNREACHABLE();
// Determine if filter is being used
bool include_filter = (filter_index >= 0);
Filter* filter = nullptr;
if (include_filter) {
filter = model::tally_filters[filter_index].get();
}
// set initial position
xyz[in_i] = origin_[in_i] - width_[0] / 2. + in_pixel / 2.;
xyz[out_i] = origin_[out_i] + width_[1] / 2. - out_pixel / 2.;
// size data array
T data(width, height, include_filter);
// arbitrary direction
Direction dir = {1. / std::sqrt(2.), 1. / std::sqrt(2.), 0.0};
// compute pixel steps and top-left pixel center
Direction u_step = u_span_ / static_cast<double>(width);
Direction v_step = v_span_ / static_cast<double>(height);
Position start =
origin_ - 0.5 * u_span_ + 0.5 * v_span_ + 0.5 * u_step - 0.5 * v_step;
// Validate that span vectors define a valid plane
Position cross = u_span_.cross(v_span_);
if (cross.norm() == 0.0) {
fatal_error("Slice span vectors are invalid (zero area).");
}
// Use an arbitrary direction that is not aligned with any coordinate axis.
// The direction has no physical meaning for plotting but is used by
// Surface::sense() to break ties when a pixel is coincident with a surface.
Direction dir = {1.0 / std::sqrt(2.0), 1.0 / std::sqrt(2.0), 0.0};
#pragma omp parallel
{
GeometryState p;
p.r() = xyz;
Particle p;
p.r() = start;
p.u() = dir;
p.coord(0).universe() = model::root_universe;
int level = slice_level_;
int j {};
FilterMatch match;
#pragma omp for
for (int y = 0; y < height; y++) {
p.r()[out_i] = xyz[out_i] - out_pixel * y;
Position row = start - v_step * static_cast<double>(y);
for (int x = 0; x < width; x++) {
p.r()[in_i] = xyz[in_i] + in_pixel * x;
p.r() = row + u_step * static_cast<double>(x);
p.n_coord() = 1;
// local variables
bool found_cell = exhaustive_find_cell(p);
@ -260,10 +277,13 @@ T SlicePlotBase::get_map() const
j = level;
}
if (found_cell) {
data.set_value(y, x, p, j);
data.set_value(y, x, p, j, filter, &match);
}
if (slice_color_overlaps_ && check_cell_overlap(p, false)) {
data.set_overlap(y, x);
if (show_overlaps_) {
int overlap_idx = check_cell_overlap(p, false);
if (overlap_idx >= 0) {
data.set_overlap(y, x, overlap_idx);
}
}
} // inner for
}
@ -297,6 +317,8 @@ public:
void print_info() const override;
PlotType type_; //!< Plot type (Slice/Voxel)
Position width_; //!< Axis-aligned width from plot.xml
PlotBasis basis_; //!< Basis from plot.xml for slice plots
int meshlines_width_; //!< Width of lines added to the plot
int index_meshlines_mesh_ {-1}; //!< Index of the mesh to draw on the plot
RGBColor meshlines_color_; //!< Color of meshlines on the plot

5
include/openmc/random_ray/flat_source_domain.h
View file

@ -76,7 +76,10 @@ public:
//----------------------------------------------------------------------------
// Static Data members
static bool volume_normalized_flux_tallies_;
static bool adjoint_; // If the user wants outputs based on the adjoint flux
// If the user wants outputs based on the adjoint flux
static bool adjoint_requested_;
// The solve currently being executed
static RandomRaySolve solve_;
static bool fw_cadis_local_;
static double
diagonal_stabilization_rho_; // Adjusts strength of diagonal stabilization

2
include/openmc/random_ray/random_ray_simulation.h
View file

@ -22,7 +22,7 @@ public:
void apply_fixed_sources_and_mesh_domains();
void prepare_fw_fixed_sources_adjoint();
void prepare_local_fixed_sources_adjoint();
void prepare_adjoint_simulation(bool fw_adjoint);
void prepare_adjoint_simulation(bool from_forward);
void simulate();
void output_simulation_results() const;
void instability_check(

6
include/openmc/random_ray/source_region.h
View file

@ -51,10 +51,10 @@ inline void hash_combine(size_t& seed, const size_t v)
// every iteration.
struct TallyTask {
int tally_idx;
int filter_idx;
int64_t filter_idx;
int score_idx;
int score_type;
TallyTask(int tally_idx, int filter_idx, int score_idx, int score_type)
TallyTask(int tally_idx, int64_t filter_idx, int score_idx, int score_type)
: tally_idx(tally_idx), filter_idx(filter_idx), score_idx(score_idx),
score_type(score_type)
{}
@ -690,7 +690,7 @@ private:
// Private Methods
// Helper function for indexing
inline int index(int64_t sr, int g) const { return sr * negroups_ + g; }
inline int64_t index(int64_t sr, int g) const { return sr * negroups_ + g; }
};
} // namespace openmc

3
include/openmc/settings.h
View file

@ -98,7 +98,8 @@ extern bool uniform_source_sampling; //!< sample sources uniformly?
extern bool ufs_on; //!< uniform fission site method on?
extern bool urr_ptables_on; //!< use unresolved resonance prob. tables?
extern bool use_decay_photons; //!< use decay photons for D1S
extern "C" bool weight_windows_on; //!< are weight windows are enabled?
extern bool use_shared_secondary_bank; //!< Use shared bank for secondaries
extern "C" bool weight_windows_on; //!< are weight windows are enabled?
extern bool weight_window_checkpoint_surface; //!< enable weight window check
//!< upon surface crossing?
extern bool weight_window_checkpoint_collision; //!< enable weight window check

56
include/openmc/shared_array.h
View file

@ -4,6 +4,8 @@
//! \file shared_array.h
//! \brief Shared array data structure
#include <algorithm> // for copy_n
#include "openmc/memory.h"
namespace openmc {
@ -30,14 +32,12 @@ public:
//! Default constructor.
SharedArray() = default;
//! Construct a zero size container with space to hold capacity number of
//! elements.
//! Construct a container with `size` elements and capacity equal to `size`.
//
//! \param capacity The number of elements for the container to allocate
//! space for
SharedArray(int64_t capacity) : capacity_(capacity)
//! \param size The number of elements to allocate and initialize
SharedArray(int64_t size) : size_(size), capacity_(size)
{
data_ = make_unique<T[]>(capacity);
data_ = make_unique<T[]>(size);
}
//==========================================================================
@ -97,8 +97,52 @@ public:
capacity_ = 0;
}
//! Push back an element to the array, with capacity and reallocation behavior
//! as if this were a vector. This does not perform any thread safety checks.
//! If the size exceeds the capacity, then the capacity will double just as
//! with a vector. Data will be reallocated and moved to a new pointer and
//! copied in before the new item is appended. Old data will be freed.
void thread_unsafe_append(const T& value)
{
if (size_ == capacity_) {
int64_t new_capacity = capacity_ == 0 ? 8 : 2 * capacity_;
unique_ptr<T[]> new_data = make_unique<T[]>(new_capacity);
std::copy_n(data_.get(), size_, new_data.get());
data_ = std::move(new_data);
capacity_ = new_capacity;
}
data_[size_++] = value;
}
//! Increase the size of the container by count elements without assigning
//! values to the new elements. Existing elements are preserved if the
//! container needs to grow. This does not perform any thread safety checks.
//
//! \param count The number of elements to append
//! \return The starting index of the appended range
int64_t extend_uninitialized(int64_t count)
{
int64_t offset = size_;
int64_t new_size = size_ + count;
if (new_size > capacity_) {
int64_t new_capacity = capacity_ == 0 ? 8 : capacity_;
while (new_capacity < new_size) {
new_capacity *= 2;
}
unique_ptr<T[]> new_data = make_unique<T[]>(new_capacity);
if (size_ > 0) {
std::copy_n(data_.get(), size_, new_data.get());
}
data_ = std::move(new_data);
capacity_ = new_capacity;
}
size_ = new_size;
return offset;
}
//! Return the number of elements in the container
int64_t size() { return size_; }
int64_t size() const { return size_; }
//! Resize the container to contain a specified number of elements. This is
//! useful in cases where the container is written to in a non-thread safe

33
include/openmc/simulation.h
View file

@ -49,6 +49,9 @@ extern const RegularMesh* ufs_mesh;
extern vector<double> k_generation;
extern vector<int64_t> work_index;
extern int64_t
simulation_tracks_completed; //!< Number of tracks completed on this rank
} // namespace simulation
//==============================================================================
@ -59,7 +62,7 @@ extern vector<int64_t> work_index;
void allocate_banks();
//! Determine number of particles to transport per process
void calculate_work();
void calculate_work(int64_t n_particles);
//! Initialize nuclear data before a simulation
void initialize_data();
@ -70,8 +73,9 @@ void initialize_batch();
//! Initialize a fission generation
void initialize_generation();
//! Full initialization of a particle history
void initialize_history(Particle& p, int64_t index_source);
//! Full initialization of a particle track
void initialize_particle_track(
Particle& p, int64_t index_source, bool is_secondary);
//! Finalize a batch
//!
@ -92,16 +96,35 @@ void broadcast_results();
void free_memory_simulation();
//! Simulate a single particle history (and all generated secondary particles,
//! if enabled), from birth to death
//! Compute unique particle ID from a 1-based source index
//! \param index_source 1-based source index within this rank's work
//! \return globally unique particle ID
int64_t compute_particle_id(int64_t index_source);
//! Compute the transport RNG seed from a particle ID
//! \param particle_id the particle's globally unique ID
//! \return seed value passed to init_particle_seeds()
int64_t compute_transport_seed(int64_t particle_id);
//! Simulate a single particle history from birth to death, inclusive of any
//! secondary particles. In shared secondary mode, only a single track is
//! transported and secondaries are deposited into a shared bank instead.
void transport_history_based_single_particle(Particle& p);
//! Simulate all particle histories using history-based parallelism
void transport_history_based();
//! Simulate all particles using history-based parallelism, with a shared
//! secondary bank
void transport_history_based_shared_secondary();
//! Simulate all particle histories using event-based parallelism
void transport_event_based();
//! Simulate all particles using event-based parallelism, with a shared
//! secondary bank
void transport_event_based_shared_secondary();
} // namespace openmc
#endif // OPENMC_SIMULATION_H

136
include/openmc/source.h
View file

@ -7,9 +7,12 @@
#include <atomic>
#include <limits>
#include <unordered_set>
#include <utility> // for pair
#include "pugixml.hpp"
#include "openmc/array.h"
#include "openmc/distribution.h"
#include "openmc/distribution_multi.h"
#include "openmc/distribution_spatial.h"
#include "openmc/memory.h"
@ -260,6 +263,139 @@ private:
vector<unique_ptr<IndependentSource>> sources_; //!< Source distributions
};
//==============================================================================
//! Parametric tokamak plasma neutron source
//!
//! This source samples neutron positions from a tokamak plasma geometry using
//! Miller-style flux surface parameterization with user-specified emission
//! profiles and energy distributions.
//!
//! Flux surface parameterization:
//! R = R0 + r*cos(alpha + delta*sin(alpha)) + Delta*(1 - (r/a)^2)
//! Z = kappa * r * sin(alpha)
//!
//! The sampling algorithm:
//! 1. Sample minor radius r from precomputed CDF of S(r) * Jacobian
//! 2. Sample poloidal angle alpha from conditional P(alpha|r) using mixture
//! of precomputed CDFs weighted by functions of r
//! 3. Sample energy and time from user-provided distribution(s)
//! 4. Sample isotropic direction
//! 5. Sample toroidal angle phi uniformly in [phi_start, phi_start +
//! phi_extent]
//! 6. Transform (r, alpha, phi) to Cartesian (x, y, z), applying the optional
//! vertical shift of the plasma center
//!
//! The user provides the emission density S(r) directly (e.g., from transport
//! codes like TRANSP, ASTRA, etc.), allowing full flexibility in reaction
//! physics calculations. S(r) is a profile in arbitrary units sampled on the
//! r_over_a grid; only its shape matters, since it is normalized internally.
//! Energy distributions can be specified as either a single distribution for
//! all r, or one distribution per radial point.
//==============================================================================
class TokamakSource : public Source {
public:
// Constructors
explicit TokamakSource(pugi::xml_node node);
//! Sample from the tokamak source distribution
//! \param[inout] seed Pseudorandom seed pointer
//! \return Sampled site
SourceSite sample(uint64_t* seed) const override;
private:
//==========================================================================
// Private methods
//! Precompute data structures for efficient sampling
void precompute_sampling_distributions();
//! Sample minor radius from marginal CDF
//! \param seed Pseudorandom seed pointer
//! \return Sampled r/a value
double sample_r_over_a(uint64_t* seed) const;
//! Sample poloidal angle given r using mixture of precomputed CDFs
//! \param r_norm Normalized minor radius r/a
//! \param seed Pseudorandom seed pointer
//! \return Sampled poloidal angle alpha [rad]
double sample_poloidal_angle(double r_norm, uint64_t* seed) const;
//! Compute the k-th mixture weight w_k(r) * I_hat_k for poloidal sampling
//! \param k Basis function index (0-5)
//! \param r Normalized minor radius r/a
//! \return Mixture weight for component k
double mixture_weight(int k, double r) const;
//! Sample energy from the distribution(s)
//! \param r_norm Normalized minor radius r/a (for distribution selection)
//! \param seed Pseudorandom seed pointer
//! \return (Sampled energy [eV], importance weight)
std::pair<double, double> sample_energy(double r_norm, uint64_t* seed) const;
//! Transform from flux coordinates (r, alpha, phi) to Cartesian (x, y, z)
//! \param r Minor radius [cm]
//! \param alpha Poloidal angle [rad]
//! \param phi Toroidal angle [rad]
//! \return Position in Cartesian coordinates [cm]
Position flux_to_cartesian(double r, double alpha, double phi) const;
//==========================================================================
// Data members
// Emission profile (input)
vector<double> r_over_a_; //!< Normalized minor radius grid points
vector<double> emission_density_; //!< Emission density S(r) at grid points
// Energy distribution(s): either 1 for all r, or one per r point
vector<unique_ptr<Distribution>> energy_dists_;
// Time distribution (defaults to a delta distribution at t=0)
UPtrDist time_;
// Angular distribution (isotropic)
UPtrAngle angle_;
// Tokamak geometry parameters
double major_radius_; //!< Major radius R0 [cm]
double minor_radius_; //!< Minor radius a [cm]
double elongation_; //!< Elongation kappa
double triangularity_; //!< Triangularity delta
double shafranov_shift_; //!< Shafranov shift Delta [cm]
double vertical_shift_; //!< Vertical shift of plasma center [cm]
// Normalized geometry parameters (precomputed for efficiency)
double epsilon_; //!< Inverse aspect ratio a/R0
double delta_tilde_; //!< Normalized Shafranov shift Delta/a
// Toroidal angle bounds
double phi_start_; //!< Starting toroidal angle [rad]
double phi_extent_; //!< Toroidal angle extent [rad]
// Precomputed distribution for radial sampling
unique_ptr<Tabular> radial_dist_;
// Coefficients of the radial geometric polynomial: A*r - B*r^2 - C*r^3
// Also used as the analytical normalization for poloidal mixture weights
double radial_poly_a_; //!< 1 + ε·Δ̃
double radial_poly_b_; //!< (3/8)·c₁·ε
double radial_poly_c_; //!< 2·ε·Δ̃
// Precomputed Bernstein basis functions for poloidal angle sampling.
// Using the factorization f(r_tilde, alpha) = R_tilde x J_tilde where:
// R_tilde = b0*(1-r)^2 + 2*b1*r*(1-r) + b2*r^2 (Bernstein quadratic)
// J_tilde = b3*(1-r) + b4*r (Bernstein linear)
// The product gives 6 non-negative basis functions g_k(alpha) with
// weights w_k(r_tilde) that are products of Bernstein polynomials.
// Distributions are tabulated on [0, pi] exploiting up-down symmetry.
static constexpr int N_POLOIDAL_BASIS = 6; //!< Number of basis functions
int n_alpha_; //!< Number of poloidal angle grid points
array<unique_ptr<Tabular>, N_POLOIDAL_BASIS>
poloidal_dists_; //!< Distributions for each basis function g_k(alpha)
array<double, N_POLOIDAL_BASIS>
poloidal_integrals_; //!< Integrals of g_k(alpha) over [0, pi]
};
//==============================================================================
// Functions
//==============================================================================

8
include/openmc/tallies/tally.h
View file

@ -97,9 +97,9 @@ public:
//! Given already-set filters, set the stride lengths
void set_strides();
int32_t strides(int i) const { return strides_[i]; }
int64_t strides(int i) const { return strides_[i]; }
int32_t n_filter_bins() const { return n_filter_bins_; }
int64_t n_filter_bins() const { return n_filter_bins_; }
bool multiply_density() const { return multiply_density_; }
@ -184,9 +184,9 @@ private:
vector<int32_t> filters_; //!< Filter indices in global filters array
//! Index strides assigned to each filter to support 1D indexing.
vector<int32_t> strides_;
vector<int64_t> strides_;
int32_t n_filter_bins_ {0};
int64_t n_filter_bins_ {0};
//! Whether to multiply by atom density for reaction rates
bool multiply_density_ {true};

2
include/openmc/tallies/tally_scoring.h
View file

@ -41,7 +41,7 @@ public:
FilterBinIter& operator++();
int index_ {1};
int64_t index_ {1};
double weight_ {1.};
vector<FilterMatch>& filter_matches_;

8
include/openmc/weight_windows.h
View file

@ -52,8 +52,12 @@ struct WeightWindow {
double weight_cutoff {DEFAULT_WEIGHT_CUTOFF};
int max_split {10};
//! Whether the weight window is in a valid state
bool is_valid() const { return lower_weight >= 0.0; }
//! Whether the weight window is in a valid state. A non-positive lower
//! bound indicates that no weight window information exists at this
//! location (generators mark such cells with -1, and a lower bound of zero
//! conventionally turns the weight window game off in a cell, as in MCNP
//! wwinp files), in which case no weight window game is played.
bool is_valid() const { return lower_weight > 0.0; }
//! Adjust the weight window by a constant factor
void scale(double factor)

230
openmc/dagmc.py
View file

@ -36,12 +36,6 @@ class DAGMCUniverse(openmc.UniverseBase):
auto_mat_ids : bool
Set IDs automatically on initialization (True) or report overlaps in ID
space between OpenMC and UWUW materials (False)
material_overrides : dict, optional
A dictionary of material overrides. The keys are material name strings
and the values are Iterables of openmc.Material objects. If a material
name is found in the DAGMC file, the material will be replaced with the
openmc.Material object in the value.
Attributes
----------
id : int
@ -78,15 +72,6 @@ class DAGMCUniverse(openmc.UniverseBase):
The number of surfaces in the model.
.. versionadded:: 0.13.2
material_overrides : dict
A dictionary of material overrides. Keys are cell IDs; values are
iterables of :class:`openmc.Material` objects. The material assignment
of each DAGMC cell ID key will be replaced with the
:class:`~openmc.Material` object in the value. If the value contains
multiple :class:`~openmc.Material` objects, each Material in the list
will be assigned to the corresponding instance of the cell.
.. versionadded:: 0.15.1
"""
def __init__(self,
@ -94,16 +79,12 @@ class DAGMCUniverse(openmc.UniverseBase):
universe_id=None,
name='',
auto_geom_ids=False,
auto_mat_ids=False,
material_overrides=None):
auto_mat_ids=False):
super().__init__(universe_id, name)
# Initialize class attributes
self.filename = filename
self.auto_geom_ids = auto_geom_ids
self.auto_mat_ids = auto_mat_ids
self._material_overrides = {}
if material_overrides is not None:
self.material_overrides = material_overrides
def __repr__(self):
string = super().__repr__()
@ -130,47 +111,17 @@ class DAGMCUniverse(openmc.UniverseBase):
@property
def material_overrides(self):
return self._material_overrides
raise AttributeError(
"DAGMCUniverse.material_overrides has been removed. Use "
"DAGMCCell objects added via add_cell() to manage per-cell "
"material assignments.")
@material_overrides.setter
def material_overrides(self, val):
cv.check_type('material overrides', val, Mapping)
for key, value in val.items():
self.add_material_override(key, value)
def replace_material_assignment(self, material_name: str, material: openmc.Material):
"""Replace a material assignment within the DAGMC universe.
Replace the material assignment of all cells filled with a material in
the DAGMC universe. The universe must be synchronized in an initialized
Model (see :meth:`~openmc.DAGMCUniverse.sync_dagmc_cells`) before
calling this method.
.. versionadded:: 0.15.1
Parameters
----------
material_name : str
Material name to replace
material : openmc.Material
Material to replace the material_name with
"""
if material_name not in self.material_names:
raise ValueError(
f"No material with name '{material_name}' found in the DAGMC universe")
if not self.cells:
raise RuntimeError("This DAGMC universe has not been synchronized "
"on an initialized Model.")
for cell in self.cells.values():
if cell.fill is None:
continue
if isinstance(cell.fill, openmc.Iterable):
cell.fill = list(map(lambda x: material if x.name == material_name else x, cell.fill))
else:
cell.fill = material if cell.fill.name == material_name else cell.fill
raise AttributeError(
"DAGMCUniverse.material_overrides has been removed. Use "
"DAGMCCell objects added via add_cell() to manage per-cell "
"material assignments.")
def add_material_override(self, key, overrides=None):
"""Add a material override to the universe.
@ -201,7 +152,10 @@ class DAGMCUniverse(openmc.UniverseBase):
if key not in self.cells:
raise ValueError(f"Cell ID '{key}' not found in DAGMC universe")
self._material_overrides[key] = overrides
if len(overrides) == 1:
self.cells[key].fill = overrides[0]
else:
self.cells[key].fill = list(overrides)
@property
def auto_geom_ids(self):
@ -290,12 +244,6 @@ class DAGMCUniverse(openmc.UniverseBase):
memo.add(self)
# Ensure that the material overrides are up-to-date
for cell in self.cells.values():
if cell.fill is None:
continue
self.add_material_override(cell, cell.fill)
# Set xml element values
dagmc_element = ET.Element('dagmc_universe')
dagmc_element.set('id', str(self.id))
@ -307,17 +255,10 @@ class DAGMCUniverse(openmc.UniverseBase):
if self.auto_mat_ids:
dagmc_element.set('auto_mat_ids', 'true')
dagmc_element.set('filename', str(self.filename))
if self._material_overrides:
mat_element = ET.Element('material_overrides')
for key in self._material_overrides:
cell_overrides = ET.Element('cell_override')
cell_overrides.set("id", str(key))
material_element = ET.Element('material_ids')
material_element.text = ' '.join(
str(t.id) for t in self._material_overrides[key])
cell_overrides.append(material_element)
mat_element.append(cell_overrides)
dagmc_element.append(mat_element)
if self.cells:
for cell in self.cells.values():
cell_element = cell.create_xml_subelement(xml_element, memo)
dagmc_element.append(cell_element)
xml_element.append(dagmc_element)
def bounding_region(
@ -442,7 +383,7 @@ class DAGMCUniverse(openmc.UniverseBase):
return out
@classmethod
def from_xml_element(cls, elem, mats = None):
def from_xml_element(cls, elem, mats=None):
"""Generate DAGMC universe from XML element
Parameters
@ -471,21 +412,56 @@ class DAGMCUniverse(openmc.UniverseBase):
out.auto_geom_ids = bool(get_text(elem, "auto_geom_ids"))
out.auto_mat_ids = bool(get_text(elem, "auto_mat_ids"))
el_mat_override = elem.find('material_overrides')
if el_mat_override is not None:
if mats is None:
raise ValueError("Material overrides found in DAGMC universe "
"but no materials were provided to populate "
"the mapping.")
out._material_overrides = {}
for elem in el_mat_override.findall('cell_override'):
cell_id = int(get_text(elem, 'id'))
mat_ids = get_elem_list(elem, "material_ids", str) or []
mat_objs = [mats[mat_id] for mat_id in mat_ids]
out._material_overrides[cell_id] = mat_objs
has_overrides = elem.find('material_overrides') is not None
has_cells = elem.find('cell') is not None
if has_overrides and has_cells:
raise ValueError(
"DAGMCUniverse cannot specify both <material_overrides> and "
"<cell> sub-elements. Use <cell> elements only.")
if has_overrides:
warnings.warn(
"DAGMCUniverse <material_overrides> is deprecated and will be "
"removed in a future version. Use nested <cell> elements "
"instead.", DeprecationWarning, stacklevel=2)
out._parse_legacy_material_overrides(elem, mats)
elif has_cells:
out._parse_cell_overrides(elem, mats)
return out
def _parse_legacy_material_overrides(self, elem, mats):
"""Parse the deprecated <material_overrides> XML format and populate
the universe with equivalent DAGMCCell objects."""
if mats is None:
raise ValueError(
"DAGMC material overrides found but no materials were "
"provided to populate the mapping.")
mo_elem = elem.find('material_overrides')
for co_elem in mo_elem.findall('cell_override'):
cell_id = int(get_text(co_elem, 'id'))
mat_ids = co_elem.find('material_ids').text.split()
fill_objs = [mats[mid] for mid in mat_ids]
fill = fill_objs[0] if len(fill_objs) == 1 else fill_objs
if cell_id in self.cells:
raise ValueError(
f"Duplicate DAGMC cell override specified for cell {cell_id}.")
self.add_cell(DAGMCCell(cell_id=cell_id, fill=fill))
def _parse_cell_overrides(self, elem, mats):
if mats is None:
raise ValueError("DAGMC cell overrides found in DAGMC universe but "
"no materials were provided to populate the "
"mapping.")
for cell_elem in elem.findall('cell'):
cell_id = int(get_text(cell_elem, 'id'))
if cell_id in self.cells:
raise ValueError(
f"Duplicate DAGMC cell override specified for cell {cell_id}.")
DAGMCCell.from_xml_element(cell_elem, mats, self)
def _partial_deepcopy(self):
"""Clone all of the openmc.DAGMCUniverse object's attributes except for
its cells, as they are copied within the clone function. This should
@ -565,7 +541,13 @@ class DAGMCUniverse(openmc.UniverseBase):
fill = [mats_per_id[mat.id] for mat in dag_cell.fill if mat]
else:
fill = mats_per_id[dag_cell.fill.id] if dag_cell.fill else None
self.add_cell(openmc.DAGMCCell(cell_id=dag_cell_id, fill=fill))
name = dag_cell.name
if dag_cell_id in self._cells:
self._cells[dag_cell_id].name = name
self._cells[dag_cell_id].fill = fill
else:
self.add_cell(
openmc.DAGMCCell(cell_id=dag_cell_id, name=name, fill=fill))
@add_plot_params
def plot(self, *args, **kwargs):
@ -594,6 +576,14 @@ class DAGMCCell(openmc.Cell):
DAG_parent_universe : int
The parent universe of the cell.
Notes
-----
DAGMC geometries are composed of triangulated surfaces, which means cell
volumes can in principle be computed exactly (e.g. via mesh-based
integration). Manually specifying :attr:`volume` overrides any such
calculation and may introduce inconsistencies if the value does not
accurately reflect the true geometric volume.
"""
def __init__(self, cell_id=None, name='', fill=None):
super().__init__(cell_id, name, fill, None)
@ -625,8 +615,62 @@ class DAGMCCell(openmc.Cell):
raise TypeError("plot is not available for DAGMC cells.")
def create_xml_subelement(self, xml_element, memo=None):
raise TypeError("create_xml_subelement is not available for DAGMC cells.")
if self.fill_type not in ('void', 'material', 'distribmat'):
raise TypeError("DAGMC cell overrides currently only support "
"material fills.")
if self.temperature is not None and self.fill_type not in (
'material', 'distribmat'
):
raise TypeError("DAGMC cell temperature overrides require a "
"material fill.")
if self.density is not None and self.fill_type not in ('material', 'distribmat'):
raise TypeError("DAGMC cell density overrides require a "
"material fill.")
if any(getattr(self, attr) is not None for attr in ('translation', 'rotation')):
raise TypeError("DAGMC cell overrides do not support translation "
"or rotation.")
return super().create_xml_subelement(xml_element, memo)
@classmethod
def from_xml_element(cls, elem, surfaces, materials, get_universe):
raise TypeError("from_xml_element is not available for DAGMC cells.")
def from_xml_element(cls, elem, mats, universe):
"""Generate a DAGMCCell from an XML <cell> override element.
Parameters
----------
elem : lxml.etree._Element
`<cell>` element containing a DAGMC cell property override
mats : dict
Dictionary mapping material ID strings to :class:`openmc.Material`
instances
universe : DAGMCUniverse
Universe to add the parsed cell to.
Returns
-------
DAGMCCell
DAGMCCell instance
"""
if not isinstance(universe, DAGMCUniverse):
raise TypeError(
f"universe must be a DAGMCUniverse instance, "
f"got {type(universe).__name__}.")
cell_id = int(get_text(elem, 'id'))
# Validate attributes that are unsupported for DAGMC cell overrides
for tag in ('region', 'fill', 'universe'):
if get_text(elem, tag) is not None:
raise ValueError(
f"DAGMC cell {cell_id} override cannot specify '{tag}'.")
for tag in ('translation', 'rotation'):
if get_text(elem, tag) is not None:
raise ValueError(
f"DAGMC cell {cell_id} override does not support "
f"'{tag}'.")
if get_elem_list(elem, 'material', str) is None:
raise ValueError(
f"DAGMC cell {cell_id} must specify a material override.")
return super().from_xml_element(
elem, surfaces={}, materials=mats,
get_universe=lambda _: universe)

7
openmc/data/angle_distribution.py
View file

@ -10,8 +10,8 @@ from openmc.mixin import EqualityMixin
from openmc.stats import Univariate, Tabular, Uniform, Legendre
from .function import INTERPOLATION_SCHEME
from .data import EV_PER_MEV
from .endf import get_head_record, get_cont_record, get_tab1_record, \
get_list_record, get_tab2_record
from .endf import as_evaluation, get_head_record, get_cont_record, \
get_tab1_record, get_list_record, get_tab2_record
class AngleDistribution(EqualityMixin):
@ -213,7 +213,7 @@ class AngleDistribution(EqualityMixin):
Parameters
----------
ev : openmc.data.endf.Evaluation
ev : openmc.data.endf.Evaluation or endf.Material
ENDF evaluation
mt : int
The MT value of the reaction to get angular distributions for
@ -224,6 +224,7 @@ class AngleDistribution(EqualityMixin):
Angular distribution
"""
ev = as_evaluation(ev)
file_obj = StringIO(ev.section[4, mt])
# Read HEAD record

60
openmc/data/data.py
View file

@ -1,14 +1,23 @@
from __future__ import annotations
import itertools
import json
import os
import re
from pathlib import Path
from math import sqrt, log
from typing import TYPE_CHECKING, Literal
from warnings import warn
from endf.data import (ATOMIC_NUMBER, ATOMIC_SYMBOL, ELEMENT_SYMBOL,
EV_PER_MEV, K_BOLTZMANN, gnds_name, zam)
import openmc
from openmc.checkvalue import PathLike
if TYPE_CHECKING:
from openmc.deplete import Chain
gnds_name.__module__ = __name__
zam.__module__ = __name__
@ -296,25 +305,47 @@ def atomic_weight(element):
raise ValueError(f"No naturally-occurring isotopes for element '{element}'.")
def half_life(isotope):
def half_life(
isotope: str,
chain_file: Literal[False] | None | PathLike | Chain = False
) -> float | None:
"""Return half-life of isotope in seconds or None if isotope is stable
Half-life values are from the `ENDF/B-VIII.0 decay sublibrary
<https://www.nndc.bnl.gov/endf-b8.0/download.html>`_.
By default, half-life values are from the `ENDF/B-VIII.0 decay sublibrary
<https://www.nndc.bnl.gov/endf-b8.0/download.html>`_. A depletion chain can
also be used as the source of half-life values.
.. versionadded:: 0.13.1
.. versionchanged:: 0.15.4
Added the ``chain_file`` argument.
Parameters
----------
isotope : str
Name of isotope, e.g., 'Pu239'
chain_file : False, None, PathLike, or openmc.deplete.Chain, optional
Source of half-life values. If ``False``, only ENDF/B-VIII.0 data is
used. If ``None``, the chain specified by
``openmc.config['chain_file']`` is used when available. If a path or
:class:`openmc.deplete.Chain` is given, that chain is used. For ``None``
or an explicit chain, nuclides absent from the chain fall back to
ENDF/B-VIII.0 data.
Returns
-------
float
Half-life of isotope in [s]
float or None
Half-life of isotope in [s], or None if the isotope is stable
"""
if chain_file is not False:
if chain_file is not None or openmc.config.get('chain_file') is not None:
# Local import avoids a circular dependency
from openmc.deplete.chain import _get_chain
chain = _get_chain(chain_file)
if isotope in chain:
return chain[isotope].half_life
global _HALF_LIFE
if not _HALF_LIFE:
# Load ENDF/B-VIII.0 data from JSON file
@ -324,7 +355,10 @@ def half_life(isotope):
return _HALF_LIFE.get(isotope.lower())
def decay_constant(isotope):
def decay_constant(
isotope: str,
chain_file: Literal[False] | None | PathLike | Chain = False
) -> float:
"""Return decay constant of isotope in [s^-1]
Decay constants are based on half-life values from the
@ -333,10 +367,20 @@ def decay_constant(isotope):
.. versionadded:: 0.13.1
.. versionchanged:: 0.15.4
Added the ``chain_file`` argument.
Parameters
----------
isotope : str
Name of isotope, e.g., 'Pu239'
chain_file : False, None, PathLike, or openmc.deplete.Chain, optional
Source of half-life values. If ``False``, only ENDF/B-VIII.0 data is
used. If ``None``, the chain specified by
``openmc.config['chain_file']`` is used when available. If a path or
:class:`openmc.deplete.Chain` is given, that chain is used. For ``None``
or an explicit chain, nuclides absent from the chain fall back to
ENDF/B-VIII.0 data.
Returns
-------
@ -348,7 +392,7 @@ def decay_constant(isotope):
openmc.data.half_life
"""
t = half_life(isotope)
t = half_life(isotope, chain_file)
return _LOG_TWO / t if t else 0.0
@ -496,5 +540,3 @@ def isotopes(element: str) -> list[tuple[str, float]]:
result.append(kv)
return result

25
openmc/data/decay.py
View file

@ -14,7 +14,8 @@ from openmc.mixin import EqualityMixin
from openmc.stats import Discrete, Tabular, Univariate, combine_distributions
from .data import gnds_name, zam
from .function import INTERPOLATION_SCHEME
from .endf import Evaluation, get_head_record, get_list_record, get_tab1_record
from .endf import (
as_evaluation, get_head_record, get_list_record, get_tab1_record)
# Gives name and (change in A, change in Z) resulting from decay
@ -75,7 +76,7 @@ class FissionProductYields(EqualityMixin):
Parameters
----------
ev_or_filename : str of openmc.data.endf.Evaluation
ev_or_filename : str, openmc.data.endf.Evaluation, or endf.Material
ENDF fission product yield evaluation to read from. If given as a
string, it is assumed to be the filename for the ENDF file.
@ -133,11 +134,7 @@ class FissionProductYields(EqualityMixin):
return energies, data
# Get evaluation if str is passed
if isinstance(ev_or_filename, Evaluation):
ev = ev_or_filename
else:
ev = Evaluation(ev_or_filename)
ev = as_evaluation(ev_or_filename)
# Assign basic nuclide properties
self.nuclide = {
@ -164,7 +161,7 @@ class FissionProductYields(EqualityMixin):
Parameters
----------
ev_or_filename : str or openmc.data.endf.Evaluation
ev_or_filename : str, openmc.data.endf.Evaluation, or endf.Material
ENDF fission product yield evaluation to read from. If given as a
string, it is assumed to be the filename for the ENDF file.
@ -292,7 +289,7 @@ class Decay(EqualityMixin):
Parameters
----------
ev_or_filename : str of openmc.data.endf.Evaluation
ev_or_filename : str, openmc.data.endf.Evaluation, or endf.Material
ENDF radioactive decay data evaluation to read from. If given as a
string, it is assumed to be the filename for the ENDF file.
@ -323,11 +320,7 @@ class Decay(EqualityMixin):
"""
def __init__(self, ev_or_filename):
# Get evaluation if str is passed
if isinstance(ev_or_filename, Evaluation):
ev = ev_or_filename
else:
ev = Evaluation(ev_or_filename)
ev = as_evaluation(ev_or_filename)
file_obj = StringIO(ev.section[8, 457])
@ -486,7 +479,7 @@ class Decay(EqualityMixin):
Parameters
----------
ev_or_filename : str or openmc.data.endf.Evaluation
ev_or_filename : str, openmc.data.endf.Evaluation, or endf.Material
ENDF radioactive decay data evaluation to read from. If given as a
string, it is assumed to be the filename for the ENDF file.
@ -649,5 +642,3 @@ def decay_energy(nuclide: str):
warn(f"Chain file '{chain_file}' does not have any decay energy.")
return _DECAY_ENERGY.get(nuclide, 0.0)

129
openmc/data/dose/dose.py
View file

@ -4,50 +4,76 @@ import numpy as np
import openmc.checkvalue as cv
_FILES = {
('icrp74', 'neutron'): Path('icrp74') / 'neutrons.txt',
('icrp74', 'photon'): Path('icrp74') / 'photons.txt',
('icrp116', 'electron'): Path('icrp116') / 'electrons.txt',
('icrp116', 'helium'): Path('icrp116') / 'helium_ions.txt',
('icrp116', 'mu-'): Path('icrp116') / 'negative_muons.txt',
('icrp116', 'pi-'): Path('icrp116') / 'negative_pions.txt',
('icrp116', 'neutron'): Path('icrp116') / 'neutrons.txt',
('icrp116', 'photon'): Path('icrp116') / 'photons.txt',
('icrp116', 'photon kerma'): Path('icrp116') / 'photons_kerma.txt',
('icrp116', 'mu+'): Path('icrp116') / 'positive_muons.txt',
('icrp116', 'pi+'): Path('icrp116') / 'positive_pions.txt',
('icrp116', 'positron'): Path('icrp116') / 'positrons.txt',
('icrp116', 'proton'): Path('icrp116') / 'protons.txt',
_FULL_GEOMETRIES = ('AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO')
_LIMITED_GEOMETRIES = ('AP', 'PA', 'ISO')
_TABLES = {
('icrp74', 'effective', 'neutron'): (
Path('icrp74') / 'neutrons.txt', _FULL_GEOMETRIES),
('icrp74', 'effective', 'photon'): (
Path('icrp74') / 'photons.txt', _FULL_GEOMETRIES),
('icrp116', 'effective', 'electron'): (
Path('icrp116') / 'electrons.txt', _LIMITED_GEOMETRIES),
('icrp116', 'effective', 'helium'): (
Path('icrp116') / 'helium_ions.txt', _LIMITED_GEOMETRIES),
('icrp116', 'effective', 'mu-'): (
Path('icrp116') / 'negative_muons.txt', _LIMITED_GEOMETRIES),
('icrp116', 'effective', 'pi-'): (
Path('icrp116') / 'negative_pions.txt', _LIMITED_GEOMETRIES),
('icrp116', 'effective', 'neutron'): (
Path('icrp116') / 'neutrons.txt', _FULL_GEOMETRIES),
('icrp116', 'effective', 'photon'): (
Path('icrp116') / 'photons.txt', _FULL_GEOMETRIES),
('icrp116', 'effective', 'photon kerma'): (
Path('icrp116') / 'photons_kerma.txt', _FULL_GEOMETRIES),
('icrp116', 'effective', 'mu+'): (
Path('icrp116') / 'positive_muons.txt', _LIMITED_GEOMETRIES),
('icrp116', 'effective', 'pi+'): (
Path('icrp116') / 'positive_pions.txt', _LIMITED_GEOMETRIES),
('icrp116', 'effective', 'positron'): (
Path('icrp116') / 'positrons.txt', _LIMITED_GEOMETRIES),
('icrp116', 'effective', 'proton'): (
Path('icrp116') / 'protons.txt', _FULL_GEOMETRIES),
('icrp74', 'ambient', 'neutron'): (
Path('icrp74') / 'neutrons_H10.txt', None),
('icrp74', 'ambient', 'photon'): (
Path('icrp74') / 'photons_H10.txt', None),
}
_DOSE_TABLES = {}
def _load_dose_icrp(data_source: str, particle: str):
"""Load effective dose tables from text files.
def _load_dose_table(data_source: str, dose_quantity: str, particle: str):
"""Load dose tables from text files.
Parameters
----------
data_source : {'icrp74', 'icrp116'}
The dose conversion data source to use
dose_quantity : {'effective', 'ambient'}
Dose quantity to load. 'ambient' corresponds to ambient dose
equivalent H*(10).
particle : {'neutron', 'photon', 'photon kerma', 'electron', 'positron'}
Incident particle
"""
path = Path(__file__).parent / _FILES[data_source, particle]
key = (data_source, dose_quantity, particle)
path = Path(__file__).parent / _TABLES[key][0]
data = np.loadtxt(path, skiprows=3, encoding='utf-8')
data[:, 0] *= 1e6 # Change energies to eV
_DOSE_TABLES[data_source, particle] = data
_DOSE_TABLES[key] = data
def dose_coefficients(particle, geometry='AP', data_source='icrp116'):
"""Return effective dose conversion coefficients.
def dose_coefficients(
particle, geometry='AP', data_source='icrp116', dose_quantity='effective'
):
"""Return dose conversion coefficients.
This function provides fluence (and air kerma) to effective or ambient dose
(H*(10)) conversion coefficients for various types of external exposures
based on values in ICRP publications. Corrected values found in a
corrigendum are used rather than the values in the original report.
Available libraries include `ICRP Publication 74
This function provides fluence (and air kerma) to effective dose or ambient
dose equivalent (H*(10)) conversion coefficients for various types of
external exposures based on values in ICRP publications. Corrected values
found in a corrigendum are used rather than the values in the original
report. Available libraries include `ICRP Publication 74
<https://doi.org/10.1016/S0146-6453(96)90010-X>` and `ICRP Publication 116
<https://doi.org/10.1016/j.icrp.2011.10.001>`.
@ -63,45 +89,58 @@ def dose_coefficients(particle, geometry='AP', data_source='icrp116'):
particle : {'neutron', 'photon', 'photon kerma', 'electron', 'positron'}
Incident particle
geometry : {'AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO'}
Irradiation geometry assumed. Refer to ICRP-116 (Section 3.2) for the
meaning of the options here.
Irradiation geometry assumed for effective dose coefficients. Refer to
ICRP-116 (Section 3.2) for the meaning of the options here. This
argument does not apply when ``dose_quantity`` is 'ambient'.
data_source : {'icrp74', 'icrp116'}
The data source for the effective dose conversion coefficients.
The data source for the dose conversion coefficients.
dose_quantity : {'effective', 'ambient'}
Dose quantity to return. 'effective' returns effective dose
coefficients; 'ambient' returns ambient dose equivalent (H*(10))
coefficients.
Returns
-------
energy : numpy.ndarray
Energies at which dose conversion coefficients are given
dose_coeffs : numpy.ndarray
Effective dose coefficients in [pSv cm^2] at provided energies. For
'photon kerma', the coefficients are given in [Sv/Gy].
Dose coefficients in [pSv cm^2] at provided energies. For 'photon
kerma', the coefficients are given in [Sv/Gy].
"""
cv.check_value('geometry', geometry, {'AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO'})
cv.check_value('geometry', geometry, _FULL_GEOMETRIES)
cv.check_value('data_source', data_source, {'icrp74', 'icrp116'})
cv.check_value('dose_quantity', dose_quantity, {'effective', 'ambient'})
if (data_source, particle) not in _FILES:
available_particles = sorted({p for (ds, p) in _FILES if ds == data_source})
key = (data_source, dose_quantity, particle)
if key not in _TABLES:
available_particles = sorted(
p for ds, dq, p in _TABLES
if ds == data_source and dq == dose_quantity
)
msg = (
f"'{particle}' has no dose data in data source {data_source}. "
f"Available particles for {data_source} are: {available_particles}"
f"'{particle}' has no {dose_quantity} dose data in data source "
f"{data_source}. Available particles for {data_source} "
f"with dose quantity {dose_quantity} are: {available_particles}"
)
raise ValueError(msg)
elif (data_source, particle) not in _DOSE_TABLES:
_load_dose_icrp(data_source, particle)
elif key not in _DOSE_TABLES:
_load_dose_table(data_source, dose_quantity, particle)
# Get all data for selected particle
data = _DOSE_TABLES[data_source, particle]
data = _DOSE_TABLES[key]
columns = _TABLES[key][1]
# Determine index for selected geometry
if particle in ('neutron', 'photon', 'proton', 'photon kerma'):
columns = ('AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO')
if columns is None:
if geometry != 'AP':
raise ValueError(
"Irradiation geometry is not defined for ambient dose "
"equivalent coefficients. Use the default geometry='AP'."
)
index = 0
else:
columns = ('AP', 'PA', 'ISO')
index = columns.index(geometry)
index = columns.index(geometry)
# Pull out energy and dose from table
energy = data[:, 0].copy()
dose_coeffs = data[:, index + 1].copy()
return energy, dose_coeffs

50
openmc/data/dose/icrp74/neutrons_H10.txt Normal file
View file

@ -0,0 +1,50 @@
Neutrons: Ambient per fluence, in units of pSv cm², for monoenergetic particles incident.
Energy (MeV) Dose
1.00E-09 6.60
1.00E-08 9.00
2.53E-08 10.6
1.00E-07 12.9
2.00E-07 13.5
5.00E-07 13.6
1.00E-06 13.3
2.00E-06 12.9
5.00E-06 12.0
1.00E-05 11.3
2.00E-05 10.6
5.00E-05 9.90
1.00E-04 9.40
2.00E-04 8.90
5.00E-04 8.30
1.00E-03 7.90
2.00E-03 7.70
5.00E-03 8.00
1.00E-02 10.5
2.00E-02 16.6
3.00E-02 23.7
5.00E-02 41.1
7.00E-02 60.0
1.00E-01 88.0
1.50E-01 132
2.00E-01 170
3.00E-01 233
5.00E-01 322
7.00E-01 375
9.00E-01 400
1 416
1.2 425
2 420
3 412
4 408
5 405
6 400
7 405
8 409
9 420
10 440
12 480
14 520
15 540
16 555
18 570
20 600

28
openmc/data/dose/icrp74/photons_H10.txt Normal file
View file

@ -0,0 +1,28 @@
Photons: Ambient dose (H*10) per fluence, in units of pSv cm²
Energy (MeV) Dose
0.010 0.061
0.015 0.83
0.020 1.05
0.030 0.81
0.040 0.64
0.050 0.55
0.060 0.51
0.080 0.53
0.100 0.61
0.150 0.89
0.200 1.20
0.300 1.80
0.400 2.38
0.500 2.93
0.600 3.44
0.800 4.38
1 5.20
1.5 6.90
2 8.60
3 11.1
4 13.4
5 15.5
6 17.6
8 21.6
10 25.6

34
openmc/data/endf.py
View file

@ -12,7 +12,12 @@ import re
from .data import gnds_name
from .function import Tabulated1D
from endf.material import _LIBRARY, _SUBLIBRARY, get_materials as get_evaluations
from endf.material import (
Material,
_LIBRARY,
_SUBLIBRARY,
get_materials as get_evaluations,
)
from endf.incident_neutron import SUM_RULES
from endf.records import (
float_endf,
@ -44,7 +49,7 @@ class Evaluation:
Parameters
----------
filename_or_obj : str or file-like
filename_or_obj : str, file-like, or endf.Material
Path to ENDF file to read or an open file positioned at the start of an
ENDF material
@ -64,17 +69,25 @@ class Evaluation:
"""
def __init__(self, filename_or_obj):
self.section = {}
self.info = {}
self.target = {}
self.projectile = {}
self.reaction_list = []
if isinstance(filename_or_obj, Material):
self.section = dict(filename_or_obj.section_text)
self.section_data = filename_or_obj.section_data
self.material = filename_or_obj.MAT
self._read_header()
return
if isinstance(filename_or_obj, (str, PurePath)):
fh = open(str(filename_or_obj), 'r')
need_to_close = True
else:
fh = filename_or_obj
need_to_close = False
self.section = {}
self.info = {}
self.target = {}
self.projectile = {}
self.reaction_list = []
# Skip TPID record. Evaluators sometimes put in TPID records that are
# ill-formated because they lack MF/MT values or put them in the wrong
@ -199,3 +212,10 @@ class Evaluation:
self.target['mass_number'],
self.target['isomeric_state'])
def as_evaluation(ev_or_filename):
"""Return an object supporting OpenMC's legacy Evaluation interface."""
if isinstance(ev_or_filename, Evaluation):
return ev_or_filename
else:
return Evaluation(ev_or_filename)

7
openmc/data/fission_energy.py
View file

@ -5,7 +5,8 @@ from io import StringIO
import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
from .data import EV_PER_MEV
from .endf import get_cont_record, get_list_record, get_tab1_record, Evaluation
from .endf import (
as_evaluation, get_cont_record, get_list_record, get_tab1_record)
from .function import Function1D, Tabulated1D, Polynomial, sum_functions
@ -195,7 +196,7 @@ class FissionEnergyRelease(EqualityMixin):
Parameters
----------
ev : openmc.data.endf.Evaluation
ev : openmc.data.endf.Evaluation or endf.Material
ENDF evaluation
incident_neutron : openmc.data.IncidentNeutron
Corresponding incident neutron dataset
@ -206,7 +207,7 @@ class FissionEnergyRelease(EqualityMixin):
Fission energy release data
"""
cv.check_type('evaluation', ev, Evaluation)
ev = as_evaluation(ev)
# Check to make sure this ENDF file matches the expected isomer.
if ev.target['atomic_number'] != incident_neutron.atomic_number:

10
openmc/data/neutron.py
View file

@ -13,7 +13,8 @@ from . import HDF5_VERSION, HDF5_VERSION_MAJOR
from .ace import Library, Table, get_table, get_metadata
from .data import ATOMIC_SYMBOL, K_BOLTZMANN, EV_PER_MEV, gnds_name
from .endf import (
Evaluation, SUM_RULES, get_head_record, get_tab1_record, get_evaluations)
Evaluation, SUM_RULES, as_evaluation, get_head_record, get_tab1_record,
get_evaluations)
from .fission_energy import FissionEnergyRelease
from .function import Tabulated1D, Sum, ResonancesWithBackground
from .njoy import make_ace, make_pendf
@ -652,7 +653,7 @@ class IncidentNeutron(EqualityMixin):
Parameters
----------
ev_or_filename : openmc.data.endf.Evaluation or str
ev_or_filename : openmc.data.endf.Evaluation, endf.Material, or str
ENDF evaluation to read from. If given as a string, it is assumed to
be the filename for the ENDF file.
@ -666,10 +667,7 @@ class IncidentNeutron(EqualityMixin):
Incident neutron continuous-energy data
"""
if isinstance(ev_or_filename, Evaluation):
ev = ev_or_filename
else:
ev = Evaluation(ev_or_filename)
ev = as_evaluation(ev_or_filename)
atomic_number = ev.target['atomic_number']
mass_number = ev.target['mass_number']

10
openmc/data/njoy.py
View file

@ -259,14 +259,14 @@ broadr / %%%%%%%%%%%%%%%%%%%%%%% Doppler broaden XS %%%%%%%%%%%%%%%%%%%%%%%%%%%%
_TEMPLATE_HEATR = """
heatr / %%%%%%%%%%%%%%%%%%%%%%%%% Add heating kerma %%%%%%%%%%%%%%%%%%%%%%%%%%%%
{nendf} {nheatr_in} {nheatr} /
{mat} 4 0 0 0 /
{mat} 4 0 0 0 0 {ed}/
302 318 402 444 /
"""
_TEMPLATE_HEATR_LOCAL = """
heatr / %%%%%%%%%%%%%%%%% Add heating kerma (local photons) %%%%%%%%%%%%%%%%%%%%
{nendf} {nheatr_in} {nheatr_local} /
{mat} 4 0 0 1 /
{mat} 4 0 0 1 0 {ed}/
302 318 402 444 /
"""
@ -411,7 +411,7 @@ def make_pendf(filename, pendf='pendf', **kwargs):
def make_ace(filename, temperatures=None, acer=True, xsdir=None,
output_dir=None, pendf=False, error=0.001, broadr=True,
heatr=True, gaspr=True, purr=True, evaluation=None,
smoothing=True, **kwargs):
smoothing=True, displacement_energy=None, **kwargs):
"""Generate incident neutron ACE file from an ENDF file
File names can be passed to
@ -463,6 +463,9 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None,
indicates which evaluation should be used.
smoothing : bool, optional
If the smoothing option (ACER card 6) is on (True) or off (False).
displacement_energy : float, optional
Threshold displacement energy in [eV]. Used in HEATR to calculate
damage-energy cross section. When None, use NJOY defaults.
**kwargs
Keyword arguments passed to :func:`openmc.data.njoy.run`
@ -515,6 +518,7 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None,
# heatr
if heatr:
ed = displacement_energy if displacement_energy is not None else 0
nheatr_in = nlast
nheatr_local = nheatr_in + 1
tapeout[nheatr_local] = (output_dir / "heatr_local") if heatr is True \

22
openmc/data/photon.py
View file

@ -15,7 +15,8 @@ from openmc.mixin import EqualityMixin
from . import HDF5_VERSION, HDF5_VERSION_MAJOR
from .ace import Table, get_metadata, get_table
from .data import ATOMIC_SYMBOL, EV_PER_MEV
from .endf import Evaluation, get_head_record, get_tab1_record, get_list_record
from .endf import (
as_evaluation, get_head_record, get_tab1_record, get_list_record)
from .function import Tabulated1D
@ -272,7 +273,7 @@ class AtomicRelaxation(EqualityMixin):
Parameters
----------
ev_or_filename : str or openmc.data.endf.Evaluation
ev_or_filename : str, openmc.data.endf.Evaluation, or endf.Material
ENDF atomic relaxation evaluation to read from. If given as a
string, it is assumed to be the filename for the ENDF file.
@ -282,10 +283,7 @@ class AtomicRelaxation(EqualityMixin):
Atomic relaxation data
"""
if isinstance(ev_or_filename, Evaluation):
ev = ev_or_filename
else:
ev = Evaluation(ev_or_filename)
ev = as_evaluation(ev_or_filename)
# Atomic relaxation data is always MF=28, MT=533
if (28, 533) not in ev.section:
@ -606,10 +604,10 @@ class IncidentPhoton(EqualityMixin):
Parameters
----------
photoatomic : str or openmc.data.endf.Evaluation
photoatomic : str, openmc.data.endf.Evaluation, or endf.Material
ENDF photoatomic data evaluation to read from. If given as a string,
it is assumed to be the filename for the ENDF file.
relaxation : str or openmc.data.endf.Evaluation, optional
relaxation : str, openmc.data.endf.Evaluation, or endf.Material, optional
ENDF atomic relaxation data evaluation to read from. If given as a
string, it is assumed to be the filename for the ENDF file.
@ -619,10 +617,7 @@ class IncidentPhoton(EqualityMixin):
Photon interaction data
"""
if isinstance(photoatomic, Evaluation):
ev = photoatomic
else:
ev = Evaluation(photoatomic)
ev = as_evaluation(photoatomic)
Z = ev.target['atomic_number']
data = cls(Z)
@ -1071,7 +1066,7 @@ class PhotonReaction(EqualityMixin):
Parameters
----------
ev : openmc.data.endf.Evaluation
ev : openmc.data.endf.Evaluation or endf.Material
ENDF photo-atomic interaction data evaluation
mt : int
The MT value of the reaction to get data for
@ -1082,6 +1077,7 @@ class PhotonReaction(EqualityMixin):
Photon reaction data
"""
ev = as_evaluation(ev)
rx = cls(mt)
# Read photon cross section

7
openmc/data/reaction.py
View file

@ -13,8 +13,8 @@ from .angle_distribution import AngleDistribution
from .angle_energy import AngleEnergy
from .correlated import CorrelatedAngleEnergy
from .data import ATOMIC_SYMBOL, K_BOLTZMANN, EV_PER_MEV
from .endf import get_head_record, get_tab1_record, get_list_record, \
get_tab2_record, get_cont_record
from .endf import as_evaluation, get_head_record, get_tab1_record, \
get_list_record, get_tab2_record, get_cont_record
from .energy_distribution import EnergyDistribution, LevelInelastic, \
DiscretePhoton
from .function import Tabulated1D, Polynomial
@ -1151,7 +1151,7 @@ class Reaction(EqualityMixin):
Parameters
----------
ev : openmc.data.endf.Evaluation
ev : openmc.data.endf.Evaluation or endf.Material
ENDF evaluation
mt : int
The MT value of the reaction to get data for
@ -1162,6 +1162,7 @@ class Reaction(EqualityMixin):
Reaction data
"""
ev = as_evaluation(ev)
rx = Reaction(mt)
# Integrated cross section

7
openmc/data/resonance.py
View file

@ -7,7 +7,9 @@ import pandas as pd
import openmc.checkvalue as cv
from .data import NEUTRON_MASS
from .endf import get_head_record, get_cont_record, get_tab1_record, get_list_record
from .endf import (
as_evaluation, get_head_record, get_cont_record, get_tab1_record,
get_list_record)
try:
from .reconstruct import wave_number, penetration_shift, reconstruct_mlbw, \
reconstruct_slbw, reconstruct_rm
@ -77,7 +79,7 @@ class Resonances:
Parameters
----------
ev : openmc.data.endf.Evaluation
ev : openmc.data.endf.Evaluation or endf.Material
ENDF evaluation
Returns
@ -86,6 +88,7 @@ class Resonances:
Resonance data
"""
ev = as_evaluation(ev)
file_obj = io.StringIO(ev.section[2, 151])
# Determine whether discrete or continuous representation

3
openmc/data/resonance_covariance.py
View file

@ -74,7 +74,7 @@ class ResonanceCovariances(Resonances):
Parameters
----------
ev : openmc.data.endf.Evaluation
ev : openmc.data.endf.Evaluation or endf.Material
ENDF evaluation
resonances : openmc.data.Resonance object
openmc.data.Resonanance object generated from the same evaluation
@ -86,6 +86,7 @@ class ResonanceCovariances(Resonances):
Resonance covariance data
"""
ev = endf.as_evaluation(ev)
file_obj = io.StringIO(ev.section[32, 151])
# Determine whether discrete or continuous representation

7
openmc/data/thermal.py
View file

@ -1042,7 +1042,7 @@ class ThermalScattering(EqualityMixin):
Parameters
----------
ev_or_filename : openmc.data.endf.Evaluation or str
ev_or_filename : openmc.data.endf.Evaluation, endf.Material, or str
ENDF evaluation to read from. If given as a string, it is assumed to
be the filename for the ENDF file.
divide_incoherent_elastic : bool
@ -1056,10 +1056,7 @@ class ThermalScattering(EqualityMixin):
Thermal scattering data
"""
if isinstance(ev_or_filename, endf.Evaluation):
ev = ev_or_filename
else:
ev = endf.Evaluation(ev_or_filename)
ev = endf.as_evaluation(ev_or_filename)
# Read incoherent inelastic data
assert (7, 4) in ev.section, 'No MF=7, MT=4 found in thermal scattering'

12
openmc/deplete/chain.py
View file

@ -319,16 +319,16 @@ class Chain:
String arguments in ``decay_files``, ``fpy_files``, and
``neutron_files`` will be treated as file names to be read.
Alternatively, :class:`openmc.data.endf.Evaluation` instances
can be included in these arguments.
Alternatively, :class:`openmc.data.endf.Evaluation` or
``endf.Material`` instances can be included in these arguments.
Parameters
----------
decay_files : list of str or openmc.data.endf.Evaluation
decay_files : list of str, openmc.data.endf.Evaluation, or endf.Material
List of ENDF decay sub-library files
fpy_files : list of str or openmc.data.endf.Evaluation
fpy_files : list of str, openmc.data.endf.Evaluation, or endf.Material
List of ENDF neutron-induced fission product yield sub-library files
neutron_files : list of str or openmc.data.endf.Evaluation
neutron_files : list of str, openmc.data.endf.Evaluation, or endf.Material
List of ENDF neutron reaction sub-library files
reactions : iterable of str, optional
Transmutation reactions to include in the depletion chain, e.g.,
@ -363,7 +363,7 @@ class Chain:
print('Processing neutron sub-library files...')
reactions = {}
for f in neutron_files:
evaluation = openmc.data.endf.Evaluation(f)
evaluation = openmc.data.endf.as_evaluation(f)
name = evaluation.gnds_name
reactions[name] = {}
for mf, mt, nc, mod in evaluation.reaction_list:

11
openmc/deplete/d1s.py
View file

@ -186,8 +186,6 @@ def apply_time_correction(
# Apply TCF, broadcasting to the correct dimensions
tcf.shape = (1, -1, 1, 1, 1)
new_tally._sum = tally_sum * tcf
new_tally._sum_sq = tally_sum_sq * (tcf*tcf)
new_tally._mean = tally_mean * tcf
new_tally._std_dev = tally_std_dev * tcf
@ -196,6 +194,8 @@ def apply_time_correction(
if sum_nuclides:
# Sum over parent nuclides (note that when combining different bins for
# parent nuclide, we can't work directly on sum_sq)
new_tally._sum = None
new_tally._sum_sq = None
new_tally._mean = new_tally.mean.sum(axis=1).reshape(shape)
new_tally._std_dev = np.linalg.norm(new_tally.std_dev, axis=1).reshape(shape)
new_tally._derived = True
@ -203,9 +203,10 @@ def apply_time_correction(
# Remove ParentNuclideFilter
new_tally.filters.pop(i_filter)
else:
# Change shape back to (filter combinations, nuclides, scores)
new_tally._sum.shape = shape
new_tally._sum_sq.shape = shape
# Apply TCF and change shape back to (filter combinations, nuclides,
# scores)
new_tally._sum = (tally_sum * tcf).reshape(shape)
new_tally._sum_sq = (tally_sum_sq * (tcf*tcf)).reshape(shape)
new_tally._mean.shape = shape
new_tally._std_dev.shape = shape

4
openmc/deplete/independent_operator.py
View file

@ -339,8 +339,12 @@ class IndependentOperator(OpenMCOperator):
for i_nuc in nuc_index:
nuc = self.nuc_ind_map[i_nuc]
if nuc not in xs._index_nuc:
continue
for i_rx in react_index:
rx = self.rx_ind_map[i_rx]
if rx not in xs._index_rx:
continue
# Determine reaction rate by multiplying xs in [b] by flux
# in [n-cm/src] to give [(reactions/src)*b-cm/atom]

66
openmc/deplete/microxs.py
View file

@ -84,7 +84,12 @@ def get_microxs_and_flux(
reactions listed in the depletion chain file are used.
energies : iterable of float or str
Energy group boundaries in [eV] or the name of the group structure.
If left as None energies will default to [0.0, 100e6]
If left as None, no energy filter is applied to the flux tally. When
`reaction_rate_mode` is "direct", these boundaries define the output
flux and microscopic cross section energy group structure. When
`reaction_rate_mode` is "flux", these boundaries define the multigroup
flux tally used to collapse continuous-energy cross sections; returned
fluxes and microscopic cross sections are one-group.
reaction_rate_mode : {"direct", "flux"}, optional
The "direct" method tallies reaction rates directly (per energy
group). The "flux" method tallies a multigroup flux spectrum and then
@ -110,7 +115,9 @@ def get_microxs_and_flux(
reaction_rate_opts : dict, optional
When `reaction_rate_mode="flux"`, allows selecting a subset of
nuclide/reaction pairs to be computed via direct reaction-rate tallies
(per energy group). Supported keys: "nuclides", "reactions".
over one energy bin spanning the full `energies` range. Supported keys:
"nuclides", "reactions". If "reactions" are specified without
"nuclides", all selected nuclides are used.
Returns
-------
@ -139,10 +146,14 @@ def get_microxs_and_flux(
nuclides = [nuc.name for nuc in chain.nuclides
if nuc.name in nuclides_with_data]
# Set up the reaction rate and flux tallies
# Set up the reaction rate and flux tallies. When energies are omitted, no
# energy filter is needed for the transport calculation. A one-group energy
# range is still needed later if flux collapse is requested.
collapse_energies = energies
if energies is None:
energies = [0.0, 100.0e6]
if isinstance(energies, str):
energy_filter = None
collapse_energies = [0.0, 100.0e6]
elif isinstance(energies, str):
energy_filter = openmc.EnergyFilter.from_group_structure(energies)
else:
energy_filter = openmc.EnergyFilter(energies)
@ -172,8 +183,11 @@ def get_microxs_and_flux(
rr_reactions = list(reactions)
elif reaction_rate_mode == 'flux' and reaction_rate_opts:
opts = reaction_rate_opts or {}
rr_nuclides = list(opts.get('nuclides', []))
rr_reactions = list(opts.get('reactions', []))
if rr_reactions:
rr_nuclides = list(opts.get('nuclides', nuclides))
else:
rr_nuclides = list(opts.get('nuclides', []))
# Keep only requested pairs within overall sets
if rr_nuclides:
rr_nuclides = [n for n in rr_nuclides if n in set(nuclides)]
@ -182,7 +196,7 @@ def get_microxs_and_flux(
# Use 1-group energy filter for RR in flux mode
has_rr = bool(rr_nuclides and rr_reactions)
if has_rr and reaction_rate_mode == 'flux':
if has_rr and reaction_rate_mode == 'flux' and energy_filter is not None:
rr_energy_filter = openmc.EnergyFilter(
[energy_filter.values[0], energy_filter.values[-1]])
else:
@ -194,14 +208,18 @@ def get_microxs_and_flux(
model.tallies = []
for i, domain_filter in enumerate(domain_filters):
flux_tally = openmc.Tally(name=f'MicroXS flux {i}')
flux_tally.filters = [domain_filter, energy_filter]
flux_tally.filters = [domain_filter]
if energy_filter is not None:
flux_tally.filters.append(energy_filter)
flux_tally.scores = ['flux']
model.tallies.append(flux_tally)
flux_tallies.append(flux_tally)
if has_rr:
rr_tally = openmc.Tally(name=f'MicroXS RR {i}')
rr_tally.filters = [domain_filter, rr_energy_filter]
rr_tally.filters = [domain_filter]
if rr_energy_filter is not None:
rr_tally.filters.append(rr_energy_filter)
rr_tally.nuclides = rr_nuclides
rr_tally.multiply_density = False
rr_tally.scores = rr_reactions
@ -215,7 +233,8 @@ def get_microxs_and_flux(
model.export_to_model_xml()
comm.barrier()
# Reinitialize with tallies
openmc.lib.init(intracomm=comm)
output = run_kwargs.get('output', True) if run_kwargs else True
openmc.lib.init(intracomm=comm, output=output)
with TemporaryDirectory() as temp_dir:
# Indicate to run in temporary directory unless being executed through
@ -255,8 +274,12 @@ def get_microxs_and_flux(
all_flux_arrays = []
for flux_tally in flux_tallies:
# Get flux values and make energy groups last dimension
flux = flux_tally.get_reshaped_data() # (domains, groups, 1, 1)
flux = np.moveaxis(flux, 1, -1) # (domains, 1, 1, groups)
flux = flux_tally.get_reshaped_data()
if energy_filter is None:
flux = flux[..., np.newaxis] # (domains, 1, 1, groups)
else:
# (domains, groups, 1, 1) -> (domains, 1, 1, groups)
flux = np.moveaxis(flux, 1, -1)
all_flux_arrays.append(flux)
fluxes.extend(flux.squeeze((1, 2)))
@ -266,8 +289,15 @@ def get_microxs_and_flux(
for flux_arr, rr_tally in zip(all_flux_arrays, rr_tallies):
flux = flux_arr
# Get reaction rates and make energy groups last dimension
reaction_rates = rr_tally.get_reshaped_data() # (domains, groups, nuclides, reactions)
reaction_rates = np.moveaxis(reaction_rates, 1, -1) # (domains, nuclides, reactions, groups)
reaction_rates = rr_tally.get_reshaped_data()
if rr_energy_filter is None:
# (domains, nuclides, reactions) ->
# (domains, nuclides, reactions, groups)
reaction_rates = reaction_rates[..., np.newaxis]
else:
# (domains, groups, nuclides, reactions) ->
# (domains, nuclides, reactions, groups)
reaction_rates = np.moveaxis(reaction_rates, 1, -1)
# If RR is 1-group, sum flux over groups
if reaction_rate_mode == "flux":
@ -279,16 +309,20 @@ def get_microxs_and_flux(
direct_micros.extend(
MicroXS(xs_i, rr_nuclides, rr_reactions) for xs_i in xs)
# If using flux mode, compute flux-collapsed microscopic XS
if reaction_rate_mode == 'flux':
# Compute flux-collapsed microscopic XS
flux_micros = [MicroXS.from_multigroup_flux(
energies=energies,
energies=collapse_energies,
multigroup_flux=flux_i,
chain_file=chain_file,
nuclides=nuclides,
reactions=reactions
) for flux_i in fluxes]
# We need to return one-group fluxes to match the microscopic cross
# sections, which are always one-group by virtue of the collapse
fluxes = [flux.sum(keepdims=True) for flux in fluxes]
# Decide which micros to use and merge if needed
if reaction_rate_mode == 'flux' and rr_tallies:
micros = [m1.merge(m2) for m1, m2 in zip(flux_micros, direct_micros)]

411
openmc/deplete/r2s.py
View file

@ -1,13 +1,16 @@
from __future__ import annotations
from collections.abc import Sequence
from contextlib import nullcontext
import copy
from datetime import datetime
import json
from numbers import Integral
from pathlib import Path
import numpy as np
import openmc
from . import IndependentOperator, PredictorIntegrator
from .chain import Chain
from .microxs import get_microxs_and_flux, write_microxs_hdf5, read_microxs_hdf5
from .results import Results
from ..checkvalue import PathLike
@ -149,11 +152,12 @@ class R2SManager:
photon_time_indices: Sequence[int] | None = None,
output_dir: PathLike | None = None,
bounding_boxes: dict[int, openmc.BoundingBox] | None = None,
chain_file: PathLike | None = None,
chain_file: PathLike | Chain | None = None,
micro_kwargs: dict | None = None,
mat_vol_kwargs: dict | None = None,
run_kwargs: dict | None = None,
operator_kwargs: dict | None = None,
by_parent_nuclide: bool = False,
):
"""Run the R2S calculation.
@ -177,7 +181,7 @@ class R2SManager:
timesteps. For example, if two timesteps are specified, the array of
times would contain three entries, and [2] would indicate computing
photon results at the last time. A value of None indicates to run
photon transport for each time.
photon transport at each time that has a decay photon source.
output_dir : PathLike, optional
Path to directory where R2S calculation outputs will be saved. If
not provided, a timestamped directory 'r2s_YYYY-MM-DDTHH-MM-SS' is
@ -187,9 +191,10 @@ class R2SManager:
Dictionary mapping cell IDs to bounding boxes used for spatial
source sampling in cell-based R2S calculations. Required if method
is 'cell-based'.
chain_file : PathLike, optional
Path to the depletion chain XML file to use during activation. If
not provided, the default configured chain file will be used.
chain_file : PathLike or openmc.deplete.Chain, optional
Path to the depletion chain XML file or depletion chain object to
use during activation. If not provided, the default configured
chain file will be used.
micro_kwargs : dict, optional
Additional keyword arguments passed to
:func:`openmc.deplete.get_microxs_and_flux` during the neutron
@ -204,6 +209,11 @@ class R2SManager:
operator_kwargs : dict, optional
Additional keyword arguments passed to
:class:`openmc.deplete.IndependentOperator`.
by_parent_nuclide : bool, optional
Whether to score photon tallies separately for each parent
radionuclide. A :class:`~openmc.ParentNuclideFilter` is added to
tallies that do not already contain one, with bins determined from
the prepared decay photon sources.
Returns
-------
@ -216,6 +226,8 @@ class R2SManager:
# consistency (different ranks may have slightly different times)
stamp = datetime.now().strftime('%Y-%m-%dT%H-%M-%S')
output_dir = Path(comm.bcast(f'r2s_{stamp}'))
else:
output_dir = Path(output_dir)
# Set run_kwargs for the neutron transport step
if micro_kwargs is None:
@ -226,22 +238,39 @@ class R2SManager:
operator_kwargs = {}
run_kwargs.setdefault('output', False)
micro_kwargs.setdefault('run_kwargs', run_kwargs)
# If a chain file is provided, prefer it for steps 1 and 2
if chain_file is not None:
micro_kwargs.setdefault('chain_file', chain_file)
operator_kwargs.setdefault('chain_file', chain_file)
self.step1_neutron_transport(
output_dir / 'neutron_transport', mat_vol_kwargs, micro_kwargs
)
self.step2_activation(
timesteps, source_rates, timestep_units, output_dir / 'activation',
operator_kwargs=operator_kwargs
)
self.step3_photon_transport(
photon_time_indices, bounding_boxes, output_dir / 'photon_transport',
mat_vol_kwargs=mat_vol_kwargs, run_kwargs=run_kwargs
# DecaySpectrum distributions are resolved in the C++ solver using
# OPENMC_CHAIN_FILE. If a Chain object was passed, write an XML
# representation alongside the R2S outputs.
if isinstance(chain_file, Chain):
output_dir.mkdir(parents=True, exist_ok=True)
chain_path = output_dir / 'chain.xml'
if comm.rank == 0:
chain_file.export_to_xml(chain_path)
comm.barrier()
else:
chain_path = chain_file
chain_context = (
openmc.config.patch('chain_file', chain_path)
if chain_path is not None else nullcontext()
)
with chain_context:
self.step1_neutron_transport(
output_dir / 'neutron_transport', mat_vol_kwargs, micro_kwargs
)
self.step2_activation(
timesteps, source_rates, timestep_units,
output_dir / 'activation', operator_kwargs=operator_kwargs
)
self.step3_photon_source(
photon_time_indices, bounding_boxes, output_dir / 'photon_transport',
mat_vol_kwargs=mat_vol_kwargs,
)
self.step4_photon_transport(
output_dir / 'photon_transport', run_kwargs=run_kwargs,
by_parent_nuclide=by_parent_nuclide,
)
return output_dir
@ -329,7 +358,7 @@ class R2SManager:
# Run neutron transport and get fluxes and micros. Run via openmc.lib to
# maintain a consistent parallelism strategy with the activation step.
with TemporarySession():
with TemporarySession(output=False):
self.results['fluxes'], self.results['micros'] = get_microxs_and_flux(
self.neutron_model, domains, **micro_kwargs)
@ -420,23 +449,20 @@ class R2SManager:
# Get depletion results
self.results['depletion_results'] = Results(output_path)
def step3_photon_transport(
def step3_photon_source(
self,
time_indices: Sequence[int] | None = None,
bounding_boxes: dict[int, openmc.BoundingBox] | None = None,
output_dir: PathLike = 'photon_transport',
mat_vol_kwargs: dict | None = None,
run_kwargs: dict | None = None,
):
"""Run the photon transport step.
"""Create decay photon sources.
This step performs photon transport calculations using decay photon
sources created from the activated materials. For each specified time,
it creates appropriate photon sources and runs a transport calculation.
In mesh-based mode, the sources are created using the mesh material
volumes, while in cell-based mode, they are created using bounding boxes
for each cell. This step will populate the 'photon_tallies' key in the
results dictionary.
This step creates decay photon sources from the activated materials for
each specified time. In mesh-based mode, the sources are created using
mesh material volumes, while in cell-based mode, they are created using
bounding boxes for each cell. This step will populate the
'photon_sources' key in the results dictionary.
Parameters
----------
@ -446,40 +472,54 @@ class R2SManager:
timesteps. For example, if two timesteps are specified, the array of
times would contain three entries, and [2] would indicate computing
photon results at the last time. A value of None indicates to run
photon transport for each time.
photon transport at each time that has a decay photon source.
bounding_boxes : dict[int, openmc.BoundingBox], optional
Dictionary mapping cell IDs to bounding boxes used for spatial
source sampling in cell-based R2S calculations. Required if method
is 'cell-based'.
output_dir : PathLike, optional
Path to directory where photon transport outputs will be saved.
Path to directory where photon source outputs will be saved.
mat_vol_kwargs : dict, optional
Additional keyword arguments passed to
:meth:`openmc.MeshBase.material_volumes`.
run_kwargs : dict, optional
Additional keyword arguments passed to :meth:`openmc.Model.run`
during the photon transport step. By default, output is disabled.
"""
# Do not retain sources from an earlier successful call if this source
# preparation attempt fails.
self.results.pop('photon_sources', None)
# TODO: Automatically determine bounding box for each cell
if bounding_boxes is None and self.method == 'cell-based':
raise ValueError("bounding_boxes must be provided for cell-based "
"R2S calculations.")
# Set default run arguments if not provided
if run_kwargs is None:
run_kwargs = {}
run_kwargs.setdefault('output', False)
# Write out JSON file with tally IDs that can be used for loading
# results
output_dir = Path(output_dir)
output_dir.mkdir(parents=True, exist_ok=True)
# Get default time indices if not provided
# Determine and validate time indices before preparing source data.
n_steps = len(self.results['depletion_results'])
implicit_time_indices = time_indices is None
if time_indices is None:
n_steps = len(self.results['depletion_results'])
time_indices = list(range(n_steps))
else:
time_indices = list(time_indices)
if not time_indices:
raise ValueError('time_indices must contain at least one index')
normalized_indices = []
for index in time_indices:
if isinstance(index, bool) or not isinstance(index, Integral):
raise TypeError('time_indices must contain only integers')
index = int(index)
if index < -n_steps or index >= n_steps:
raise IndexError(
f'Photon time index {index} is out of range for '
f'{n_steps} depletion results')
normalized_index = index % n_steps
normalized_indices.append(normalized_index)
# Remove duplicates while preserving order
time_indices = list(dict.fromkeys(normalized_indices))
# Check whether the photon model is different
neutron_univ = self.neutron_model.geometry.root_universe
@ -505,6 +545,106 @@ class R2SManager:
self.results['mesh_material_volumes_photon'] = photon_mmv_list
# Get dictionary of cells in the photon model
if different_photon_model:
photon_cells = self.photon_model.geometry.get_all_cells()
# Determine eligible work items upfront (independent of time index).
if self.method == 'mesh-based':
work_items = self._get_mesh_work_items()
else:
work_items = []
for cell, original_mat in zip(
self.domains, self.results['activation_materials']):
if different_photon_model:
if cell.id not in photon_cells or \
cell.fill.id != photon_cells[cell.id].fill.id:
continue
work_items.append((cell, original_mat, bounding_boxes[cell.id]))
# Create decay photon sources for each time index
photon_sources = {
time_index: self._create_photon_sources(time_index, work_items)
for time_index in time_indices
}
# Determine if any times have no decay photon sources. If the user
# didn't specify any specific time indices, remove those times from the
# photon_sources dictionary. If the user did specify time indices, raise
# an error if any of those times have no decay photon sources.
empty_indices = [
time_index for time_index, sources in photon_sources.items()
if not sources
]
if implicit_time_indices:
for time_index in empty_indices:
del photon_sources[time_index]
if not photon_sources:
raise RuntimeError(
'No decay photon sources were found at any depletion time')
elif empty_indices:
indices = ', '.join(str(index) for index in empty_indices)
raise RuntimeError(
f'No decay photon source was found for requested time '
f'indices: {indices}')
self.results['photon_sources'] = photon_sources
def step4_photon_transport(
self,
output_dir: PathLike = 'photon_transport',
run_kwargs: dict | None = None,
by_parent_nuclide: bool = False,
):
"""Run photon transport using prepared decay photon sources.
This step runs a photon transport calculation for each source list
created by :meth:`step3_photon_source`. It will populate the
'photon_tallies' key in the results dictionary.
Parameters
----------
output_dir : PathLike, optional
Path to directory where photon transport outputs will be saved.
run_kwargs : dict, optional
Additional keyword arguments passed to :meth:`openmc.Model.run`.
By default, output is disabled.
by_parent_nuclide : bool, optional
Whether to score photon tallies separately for each parent
radionuclide. A :class:`~openmc.ParentNuclideFilter` is added to
tallies that do not already contain one, with bins determined from
the prepared decay photon sources.
"""
if 'photon_sources' not in self.results:
raise RuntimeError(
'Photon sources must be created with step3_photon_source '
'before running photon transport.')
photon_sources = self.results['photon_sources']
if not photon_sources:
raise RuntimeError(
'No decay photon sources are available for transport')
if by_parent_nuclide:
radionuclides = sorted({
nuclide
for sources in photon_sources.values()
for source in sources
for nuclide in source.energy.nuclides
})
if radionuclides:
parent_filter = openmc.ParentNuclideFilter(radionuclides)
for tally in self.photon_model.tallies:
if not tally.contains_filter(openmc.ParentNuclideFilter):
tally.filters.append(parent_filter)
if run_kwargs is None:
run_kwargs = {}
run_kwargs.setdefault('output', False)
output_dir = Path(output_dir)
output_dir.mkdir(parents=True, exist_ok=True)
if comm.rank == 0:
tally_ids = [tally.id for tally in self.photon_model.tallies]
with open(output_dir / 'tally_ids.json', 'w') as f:
@ -512,49 +652,11 @@ class R2SManager:
self.results['photon_tallies'] = {}
# Get dictionary of cells in the photon model
if different_photon_model:
photon_cells = self.photon_model.geometry.get_all_cells()
# Ensure photon transport is enabled in settings.
self.photon_model.settings.photon_transport = True
for time_index in time_indices:
# Create decay photon source
if self.method == 'mesh-based':
self.photon_model.settings.source = \
self.get_decay_photon_source_mesh(time_index)
else:
sources = []
results = self.results['depletion_results']
for cell, original_mat in zip(self.domains, self.results['activation_materials']):
# Skip if the cell is not in the photon model or the
# material has changed
if different_photon_model:
if cell.id not in photon_cells or \
cell.fill.id != photon_cells[cell.id].fill.id:
continue
# Get bounding box for the cell
bounding_box = bounding_boxes[cell.id]
# Get activated material composition
activated_mat = results[time_index].get_material(str(original_mat.id))
# Create decay photon source source
space = openmc.stats.Box(*bounding_box)
energy = activated_mat.get_decay_photon_energy()
strength = energy.integral() if energy is not None else 0.0
source = openmc.IndependentSource(
space=space,
energy=energy,
particle='photon',
strength=strength,
constraints={'domains': [cell]}
)
sources.append(source)
self.photon_model.settings.source = sources
# Convert time_index (which may be negative) to a normal index
if time_index < 0:
time_index = len(self.results['depletion_results']) + time_index
for time_index, sources in photon_sources.items():
self.photon_model.settings.source = sources
# Run photon transport calculation
photon_dir = Path(output_dir) / f'time_{time_index}'
@ -567,58 +669,30 @@ class R2SManager:
sp.tallies[tally.id] for tally in self.photon_model.tallies
]
def get_decay_photon_source_mesh(
self,
time_index: int = -1
) -> list[openmc.IndependentSource]:
"""Create decay photon source for a mesh-based calculation.
def _get_mesh_work_items(self):
"""Enumerate mesh-based work items across all meshes.
For each mesh element-material combination across all meshes, an
:class:`~openmc.IndependentSource` is created with a
:class:`~openmc.stats.Box` spatial distribution based on the bounding
box of the material within the mesh element. A material constraint is
also applied so that sampled source sites are limited to the correct
region.
When the photon transport model is different from the neutron model, the
photon MeshMaterialVolumes is used to determine whether an (element,
material) combination exists in the photon model.
Parameters
----------
time_index : int, optional
Time index for the decay photon source. Default is -1 (last time).
Returns a list of (index_mat, mat_id, bbox) tuples for each eligible
mesh element--material combination, where index_mat is the index into
the activation materials list, mat_id is the material ID, and bbox is
the bounding box for that mesh element--material combination.
Returns
-------
list of openmc.IndependentSource
A list of IndependentSource objects for the decay photons, one for
each mesh element-material combination with non-zero source strength.
list of tuple
Each tuple is (index_mat, mat_id, bbox).
"""
mat_dict = self.neutron_model._get_all_materials()
# List to hold all sources
sources = []
# Index in the overall list of activated materials
index_mat = 0
# Get various results from previous steps
mmv_list = self.results['mesh_material_volumes']
materials = self.results['activation_materials']
results = self.results['depletion_results']
photon_mmv_list = self.results.get('mesh_material_volumes_photon')
work_items = []
index_mat = 0
for mesh_idx, mat_vols in enumerate(mmv_list):
photon_mat_vols = photon_mmv_list[mesh_idx] \
if photon_mmv_list is not None else None
# Total number of mesh elements for this mesh
n_elements = mat_vols.num_elements
for index_elem in range(n_elements):
# Determine which materials exist in the photon model for this element
if photon_mat_vols is not None:
photon_materials = {
mat_id
@ -626,36 +700,77 @@ class R2SManager:
if mat_id is not None
}
for mat_id, _, bbox in mat_vols.by_element(index_elem, include_bboxes=True):
# Skip void volume
for mat_id, _, bbox in mat_vols.by_element(
index_elem, include_bboxes=True):
if mat_id is None:
continue
# Skip if this material doesn't exist in photon model
if photon_mat_vols is not None and mat_id not in photon_materials:
if photon_mat_vols is not None \
and mat_id not in photon_materials:
index_mat += 1
continue
# Get activated material composition
original_mat = materials[index_mat]
activated_mat = results[time_index].get_material(str(original_mat.id))
# Create decay photon source
energy = activated_mat.get_decay_photon_energy()
if energy is not None:
strength = energy.integral()
space = openmc.stats.Box(*bbox)
sources.append(openmc.IndependentSource(
space=space,
energy=energy,
particle='photon',
strength=strength,
constraints={'domains': [mat_dict[mat_id]]}
))
# Increment index of activated material
work_items.append((index_mat, mat_id, bbox))
index_mat += 1
return work_items
def _create_photon_sources(self, time_index, work_items):
"""Create decay photon sources for a set of regions.
Builds :class:`openmc.IndependentSource` objects with
:class:`openmc.stats.DecaySpectrum` energy distributions that will be
serialized to XML and resolved against the depletion chain by the C++
solver.
Parameters
----------
time_index : int
Index into depletion results.
work_items : list of tuple
For mesh-based: list of (index_mat, mat_id, bbox).
For cell-based: list of (cell, original_mat, bbox).
Returns
-------
list of openmc.IndependentSource
Photon sources for each activated region.
"""
step_result = self.results['depletion_results'][time_index]
materials = self.results['activation_materials']
mesh_based = self.method == 'mesh-based'
if mesh_based:
mat_dict = self.neutron_model._get_all_materials()
sources = []
for item in work_items:
if mesh_based:
index_mat, domain_id, bbox = item
original_mat = materials[index_mat]
domain = mat_dict[domain_id]
else:
cell, original_mat, bbox = item
domain = cell
activated_mat = step_result.get_material(str(original_mat.id))
nuclides = activated_mat.get_nuclide_atom_densities()
if not nuclides:
continue
# Eliminate nuclides with zero density
nuclides = {nuclide: density for nuclide, density in nuclides.items()
if density > 0}
energy = openmc.stats.DecaySpectrum(nuclides, activated_mat.volume)
energy.clip(inplace=True)
if not energy.nuclides:
continue
sources.append(openmc.IndependentSource(
space=openmc.stats.Box(bbox.lower_left, bbox.upper_right),
energy=energy,
particle='photon',
constraints={'domains': [domain]},
))
return sources
def load_results(self, path: PathLike):

35
openmc/deplete/results.py
View file

@ -1,12 +1,17 @@
from __future__ import annotations
import numbers
import bisect
import math
from collections.abc import Iterable
from typing import Literal
from warnings import warn
import h5py
import numpy as np
import openmc
from .chain import Chain, _get_chain
from .stepresult import StepResult, VERSION_RESULTS
import openmc.checkvalue as cv
from openmc.data import atomic_mass, AVOGADRO
@ -103,7 +108,8 @@ class Results(list):
mat: Material | str,
units: str = "Bq/cm3",
by_nuclide: bool = False,
volume: float | None = None
volume: float | None = None,
chain_file: Literal[False] | None | PathLike | Chain = None
) -> tuple[np.ndarray, np.ndarray | list[dict]]:
"""Get activity of material over time.
@ -115,22 +121,31 @@ class Results(list):
Material object or material id to evaluate
units : {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3', 'Bq/m3'}
Specifies the type of activity to return, options include total
activity [Bq], specific [Bq/g, Bq/kg] or volumetric activity [Bq/cm3].
activity [Bq], specific [Bq/g, Bq/kg] or volumetric activity
[Bq/cm3].
by_nuclide : bool
Specifies if the activity should be returned for the material as a
whole or per nuclide. Default is False.
volume : float, optional
Volume of the material. If not passed, defaults to using the
:attr:`Material.volume` attribute.
chain_file : False, None, PathLike, or openmc.deplete.Chain, optional
Source of half-life values. If ``False``, only ENDF/B-VIII.0 data is
used. If ``None``, the chain specified by
``openmc.config['chain_file']`` is used when available. If a path or
:class:`openmc.deplete.Chain` is given, that chain is used. For
``None`` or an explicit chain, nuclides absent from the chain fall
back to ENDF/B-VIII.0 data.
.. versionadded:: 0.15.4
Returns
-------
times : numpy.ndarray
Array of times in [s]
activities : numpy.ndarray or List[dict]
Array of total activities if by_nuclide = False (default)
or list of dictionaries of activities by nuclide if
by_nuclide = True.
Array of total activities if by_nuclide = False (default) or list of
dictionaries of activities by nuclide if by_nuclide = True.
"""
if isinstance(mat, Material):
@ -140,6 +155,13 @@ class Results(list):
else:
raise TypeError('mat should be of type openmc.Material or str')
if chain_file is not False:
if chain_file is None:
if openmc.config.get('chain_file') is not None:
chain_file = _get_chain(None)
else:
chain_file = _get_chain(chain_file)
times = np.empty_like(self, dtype=float)
if by_nuclide:
activities = [None] * len(self)
@ -149,7 +171,8 @@ class Results(list):
# Evaluate activity for each depletion time
for i, result in enumerate(self):
times[i] = result.time[0]
activities[i] = result.get_material(mat_id).get_activity(units, by_nuclide, volume)
activities[i] = result.get_material(mat_id).get_activity(
units, by_nuclide, volume, chain_file=chain_file)
return times, activities

4
openmc/examples.py
View file

@ -150,7 +150,7 @@ def pwr_core() -> openmc.Model:
rpv_steel.add_nuclide('Ni60', 0.0026776, 'wo')
rpv_steel.add_nuclide('Mn55', 0.01, 'wo')
rpv_steel.add_nuclide('Cr52', 0.002092475, 'wo')
rpv_steel.add_nuclide('C0', 0.0025, 'wo')
rpv_steel.add_element('C', 0.0025, 'wo')
rpv_steel.add_nuclide('Cu63', 0.0013696, 'wo')
lower_rad_ref = openmc.Material(6, name='Lower radial reflector')
@ -1618,4 +1618,4 @@ def random_ray_three_region_cube_with_detectors() -> openmc.Model:
model.settings = settings
model.tallies = tallies
return model
return model

5
openmc/lib/core.py
View file

@ -31,7 +31,10 @@ class _SourceSite(Structure):
('particle', c_int32),
('parent_nuclide', c_int),
('parent_id', c_int64),
('progeny_id', c_int64)]
('progeny_id', c_int64),
('wgt_born', c_double),
('wgt_ww_born', c_double),
('n_split', c_int64)]
# Define input type for numpy arrays that will be passed into C++ functions
# Must be an int or double array, with single dimension that is contiguous

404
openmc/lib/plot.py
View file

@ -9,6 +9,7 @@ from .core import _FortranObjectWithID
from .error import _error_handler
import numpy as np
import warnings
class _Position(Structure):
@ -51,217 +52,260 @@ class _Position(Structure):
return f"({self.x}, {self.y}, {self.z})"
class _PlotBase(Structure):
"""A structure defining a 2-D geometry slice with underlying c-types
def _extract_slice_data_args(plot):
"""Convert a legacy plot-like object into slice_data keyword arguments."""
try:
kwargs = {
'origin': tuple(plot.origin),
'width': (plot.width, plot.height),
'basis': plot.basis,
'pixels': (plot.h_res, plot.v_res),
'show_overlaps': getattr(plot, 'color_overlaps', False),
'level': getattr(plot, 'level', -1),
}
except AttributeError as exc:
raise TypeError(
"plot must be a legacy plot-like object with origin, width, "
"height, basis, h_res, and v_res attributes."
) from exc
return kwargs
C-Type Attributes
-----------------
origin_ : openmc.lib.plot._Position
A position defining the origin of the plot.
width_ : openmc.lib.plot._Position
The width of the plot along the x, y, and z axes, respectively
basis_ : c_int
The axes basis of the plot view.
pixels_ : c_size_t[3]
The resolution of the plot in the horizontal and vertical dimensions
color_overlaps_ : c_bool
Whether to assign unique IDs (-3) to overlapping regions.
level_ : c_int
The universe level for the plot view
Attributes
_dll.openmc_slice_data.argtypes = [
POINTER(c_double * 3), # origin
POINTER(c_double * 3), # u_span
POINTER(c_double * 3), # v_span
POINTER(c_size_t * 2), # pixels
c_bool, # show_overlaps
c_int, # level
c_int32, # filter_index
POINTER(c_int32), # geom_data
POINTER(c_double), # property_data (can be None)
]
_dll.openmc_slice_data.restype = c_int
_dll.openmc_slice_data.errcheck = _error_handler
def slice_data(origin, width=None, basis='xy', u_span=None, v_span=None,
pixels=None, show_overlaps=False, level=None, filter=None,
include_properties=True):
"""Generate a 2D raster of geometry and property data for plotting.
Parameters
----------
origin : tuple or list of ndarray
Origin (center) of the plot
width : float
The horizontal dimension of the plot in geometry units (cm)
height : float
The vertical dimension of the plot in geometry units (cm)
basis : string
One of {'xy', 'xz', 'yz'} indicating the horizontal and vertical
axes of the plot.
h_res : int
The horizontal resolution of the plot in pixels
v_res : int
The vertical resolution of the plot in pixels
level : int
The universe level for the plot (default: -1 -> all universes shown)
origin : sequence of float
Center position of the plot [x, y, z]
width : sequence of float
Width of the plot [horizontal, vertical]. Mutually exclusive with
u_span/v_span.
basis : {'xy', 'xz', 'yz'} or int
Plot basis. Ignored if u_span/v_span are provided.
u_span : sequence of float, optional
Full-width span vector for the horizontal axis (3 values). Mutually
exclusive with width.
v_span : sequence of float, optional
Full-height span vector for the vertical axis (3 values). Mutually
exclusive with width.
pixels : sequence of int
Number of pixels [horizontal, vertical]
show_overlaps : bool, optional
Whether to detect overlapping cells
level : int, optional
Universe level (None for deepest)
filter : openmc.lib.Filter, optional
Filter for bin index lookup
include_properties : bool, optional
Whether to compute temperature/density
Returns
-------
geom_data : numpy.ndarray
Array of shape (v_res, h_res, 3) or (v_res, h_res, 4) with int32 dtype.
Contains [cell_id, cell_instance, material_id] when no filter is provided,
or [cell_id, cell_instance, material_id, filter_bin] when a filter is provided.
property_data : numpy.ndarray or None
Array of shape (v_res, h_res, 2) with float64 dtype containing
[temperature, density], or None if include_properties=False
"""
_fields_ = [('origin_', _Position),
('width_', _Position),
('basis_', c_int),
('pixels_', 3*c_size_t),
('color_overlaps_', c_bool),
('level_', c_int)]
# Set deepest level as default
if level is None:
level = -1
if not isinstance(level, int):
raise TypeError("level must be an integer.")
def __init__(self):
self.level_ = -1
self.basis_ = 1
self.color_overlaps_ = False
if pixels is None:
raise ValueError("pixels must be specified.")
if len(pixels) != 2:
raise ValueError("pixels must be a length-2 sequence.")
@property
def origin(self):
return self.origin_
if width is not None and (u_span is not None or v_span is not None):
raise ValueError("width is mutually exclusive with u_span/v_span.")
@origin.setter
def origin(self, origin):
self.origin_.x = origin[0]
self.origin_.y = origin[1]
self.origin_.z = origin[2]
@property
def width(self):
return self.width_.x
@width.setter
def width(self, width):
self.width_.x = width
@property
def height(self):
return self.width_.y
@height.setter
def height(self, height):
self.width_.y = height
@property
def basis(self):
if self.basis_ == 1:
return 'xy'
elif self.basis_ == 2:
return 'xz'
elif self.basis_ == 3:
return 'yz'
raise ValueError(f"Plot basis {self.basis_} is invalid")
@basis.setter
def basis(self, basis):
if u_span is not None or v_span is not None:
if u_span is None or v_span is None:
raise ValueError("Both u_span and v_span must be provided.")
u_span = np.asarray(u_span, dtype=float)
v_span = np.asarray(v_span, dtype=float)
if u_span.shape != (3,) or v_span.shape != (3,):
raise ValueError("u_span and v_span must be length-3 sequences.")
u_norm = np.linalg.norm(u_span)
v_norm = np.linalg.norm(v_span)
if u_norm == 0.0 or v_norm == 0.0:
raise ValueError("u_span and v_span must be non-zero vectors.")
dot = float(np.dot(u_span, v_span))
ortho_tol = 1.0e-10 * u_norm * v_norm
if abs(dot) > ortho_tol:
raise ValueError("u_span and v_span must be orthogonal.")
else:
if width is None:
raise ValueError("width must be provided when u_span/v_span are not set.")
if len(width) != 2:
raise ValueError("width must be a length-2 sequence.")
basis_map = {'xy': 1, 'xz': 2, 'yz': 3}
if isinstance(basis, str):
valid_bases = ('xy', 'xz', 'yz')
basis = basis.lower()
if basis not in valid_bases:
if basis not in basis_map:
raise ValueError(f"{basis} is not a valid plot basis.")
if basis == 'xy':
self.basis_ = 1
elif basis == 'xz':
self.basis_ = 2
elif basis == 'yz':
self.basis_ = 3
return
if isinstance(basis, int):
valid_bases = (1, 2, 3)
if basis not in valid_bases:
basis = basis_map[basis]
elif isinstance(basis, int):
if basis not in basis_map.values():
raise ValueError(f"{basis} is not a valid plot basis.")
self.basis_ = basis
return
else:
raise ValueError(f"{basis} is not a valid plot basis.")
raise ValueError(f"{basis} of type {type(basis)} is an invalid plot basis")
if basis == 1:
u_span = np.array([width[0], 0.0, 0.0], dtype=float)
v_span = np.array([0.0, width[1], 0.0], dtype=float)
elif basis == 2:
u_span = np.array([width[0], 0.0, 0.0], dtype=float)
v_span = np.array([0.0, 0.0, width[1]], dtype=float)
else:
u_span = np.array([0.0, width[0], 0.0], dtype=float)
v_span = np.array([0.0, 0.0, width[1]], dtype=float)
@property
def h_res(self):
return self.pixels_[0]
origin = np.asarray(origin, dtype=float)
if origin.shape != (3,):
raise ValueError("origin must be a length-3 sequence.")
@h_res.setter
def h_res(self, h_res):
self.pixels_[0] = h_res
# Prepare ctypes arrays
origin_arr = (c_double * 3)(*origin)
u_span_arr = (c_double * 3)(*u_span)
v_span_arr = (c_double * 3)(*v_span)
pixels_arr = (c_size_t * 2)(*pixels)
@property
def v_res(self):
return self.pixels_[1]
# Get internal filter index from filter ID if filter is provided
if filter is not None:
filter_index = c_int32()
_dll.openmc_get_filter_index(filter.id, filter_index)
filter_index = filter_index.value
else:
filter_index = -1
@v_res.setter
def v_res(self, v_res):
self.pixels_[1] = v_res
# Allocate output arrays with dynamic size based on filter
n_geom_fields = 4 if filter is not None else 3
geom_data = np.zeros((pixels[1], pixels[0], n_geom_fields), dtype=np.int32)
if include_properties:
property_data = np.zeros((pixels[1], pixels[0], 2), dtype=np.float64)
prop_ptr = property_data.ctypes.data_as(POINTER(c_double))
else:
property_data = None
prop_ptr = None
@property
def level(self):
return int(self.level_)
_dll.openmc_slice_data(
origin_arr,
u_span_arr,
v_span_arr,
pixels_arr,
show_overlaps,
level,
filter_index,
geom_data.ctypes.data_as(POINTER(c_int32)),
prop_ptr
)
@level.setter
def level(self, level):
self.level_ = level
@property
def color_overlaps(self):
return self.color_overlaps_
@color_overlaps.setter
def color_overlaps(self, color_overlaps):
self.color_overlaps_ = color_overlaps
def __repr__(self):
out_str = ["-----",
"Plot:",
"-----",
f"Origin: {self.origin}",
f"Width: {self.width}",
f"Height: {self.height}",
f"Basis: {self.basis}",
f"HRes: {self.h_res}",
f"VRes: {self.v_res}",
f"Color Overlaps: {self.color_overlaps}",
f"Level: {self.level}"]
return '\n'.join(out_str)
_dll.openmc_id_map.argtypes = [POINTER(_PlotBase), POINTER(c_int32)]
_dll.openmc_id_map.restype = c_int
_dll.openmc_id_map.errcheck = _error_handler
return geom_data, property_data
def id_map(plot):
"""Deprecated compatibility wrapper for geometry ID maps.
This function is kept for compatibility and will be removed in a future
release. Use `slice_data(..., include_properties=False)` instead.
"""
Generate a 2-D map of cell and material IDs. Used for in-memory image
generation.
warnings.warn(
"openmc.lib.id_map is deprecated and will be removed in a future "
"release; use openmc.lib.slice_data(..., include_properties=False).",
FutureWarning,
)
Parameters
----------
plot : openmc.lib.plot._PlotBase
Object describing the slice of the model to be generated
Returns
-------
id_map : numpy.ndarray
A NumPy array with shape (vertical pixels, horizontal pixels, 3) of
OpenMC property ids with dtype int32. The last dimension of the array
contains, in order, cell IDs, cell instances, and material IDs.
"""
img_data = np.zeros((plot.v_res, plot.h_res, 3),
dtype=np.dtype('int32'))
_dll.openmc_id_map(plot, img_data.ctypes.data_as(POINTER(c_int32)))
return img_data
_dll.openmc_property_map.argtypes = [POINTER(_PlotBase), POINTER(c_double)]
_dll.openmc_property_map.restype = c_int
_dll.openmc_property_map.errcheck = _error_handler
kwargs = _extract_slice_data_args(plot)
geom_data, _ = slice_data(include_properties=False, **kwargs)
return geom_data[:, :, :3]
def property_map(plot):
"""
Generate a 2-D map of cell temperatures and material densities. Used for
in-memory image generation.
"""Deprecated compatibility wrapper for temperature/density maps.
Parameters
----------
plot : openmc.lib.plot._PlotBase
Object describing the slice of the model to be generated
This function is kept for compatibility and will be removed in a future
release. Use `slice_data(..., include_properties=True)` instead.
"""
warnings.warn(
"openmc.lib.property_map is deprecated and will be removed in a "
"future release; use openmc.lib.slice_data(..., "
"include_properties=True).",
FutureWarning,
)
kwargs = _extract_slice_data_args(plot)
_, prop_data = slice_data(include_properties=True, **kwargs)
return prop_data
_dll.openmc_slice_data_overlap_count.argtypes = [POINTER(c_size_t)]
_dll.openmc_slice_data_overlap_count.restype = c_int
_dll.openmc_slice_data_overlap_count.errcheck = _error_handler
_dll.openmc_slice_data_overlap_info.argtypes = [c_size_t, POINTER(c_int32)]
_dll.openmc_slice_data_overlap_info.restype = c_int
_dll.openmc_slice_data_overlap_info.errcheck = _error_handler
# Python wrappings for overlap functions
def slice_data_overlap_count() -> int:
"""Return the number of unique overlaps from the last slice plot.
Returns
-------
property_map : numpy.ndarray
A NumPy array with shape (vertical pixels, horizontal pixels, 2) of
OpenMC property ids with dtype float
int
Number of unique overlapping cell pairs detected by the most recent
:func:`slice_data` call with overlap checking enabled.
"""
prop_data = np.zeros((plot.v_res, plot.h_res, 2))
_dll.openmc_property_map(plot, prop_data.ctypes.data_as(POINTER(c_double)))
return prop_data
count = c_size_t()
_dll.openmc_slice_data_overlap_count(count)
return count.value
def slice_data_overlap_info() -> np.ndarray:
"""Return identifying information for overlaps from the last slice plot.
Returns
-------
numpy.ndarray
Array of shape ``(n, 3)`` with int32 dtype, where ``n`` is the number
of unique overlaps detected by the most recent :func:`slice_data` call
with overlap checking enabled. Each row contains ``[universe_id,
cell1_id, cell2_id]``.
"""
n = slice_data_overlap_count()
overlap_info = np.empty((n, 3), dtype=np.int32)
if n > 0:
_dll.openmc_slice_data_overlap_info(
n,
overlap_info.ctypes.data_as(POINTER(c_int32)),
)
return overlap_info
_dll.openmc_get_plot_index.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_get_plot_index.restype = c_int

50
openmc/material.py
View file

@ -8,7 +8,7 @@ from pathlib import Path
import re
import sys
import tempfile
from typing import Sequence, Dict
from typing import TYPE_CHECKING, Literal, Sequence, Dict
import warnings
import lxml.etree as ET
@ -28,6 +28,9 @@ from openmc.data.data import _get_element_symbol, JOULE_PER_EV
from openmc.data.function import Tabulated1D
from openmc.data import mass_energy_absorption_coefficient, dose_coefficients
if TYPE_CHECKING:
from openmc.deplete import Chain
# Units for density supported by OpenMC
DENSITY_UNITS = ('g/cm3', 'g/cc', 'kg/m3', 'atom/b-cm', 'atom/cm3', 'sum',
@ -296,6 +299,8 @@ class Material(IDManagerMixin):
mass += nuc.percent
# Compute and return the molar mass
if moles == 0.0:
raise ValueError("Material has no nuclides; cannot compute molar mass")
return mass / moles
@property
@ -1383,8 +1388,13 @@ class Material(IDManagerMixin):
return densities
def get_activity(self, units: str = 'Bq/cm3', by_nuclide: bool = False,
volume: float | None = None) -> dict[str, float] | float:
def get_activity(
self,
units: str = 'Bq/cm3',
by_nuclide: bool = False,
volume: float | None = None,
chain_file: Literal[False] | None | PathLike | Chain = None
) -> dict[str, float] | float:
"""Return the activity of the material or each nuclide within.
.. versionadded:: 0.13.1
@ -1403,13 +1413,22 @@ class Material(IDManagerMixin):
:attr:`Material.volume` attribute.
.. versionadded:: 0.13.3
chain_file : False, None, PathLike, or openmc.deplete.Chain, optional
Source of half-life values. If ``False``, only ENDF/B-VIII.0 data is
used. If ``None``, the chain specified by
``openmc.config['chain_file']`` is used when available. If a path or
:class:`openmc.deplete.Chain` is given, that chain is used. For
``None`` or an explicit chain, nuclides absent from the chain fall
back to ENDF/B-VIII.0 data.
.. versionadded:: 0.15.4
Returns
-------
Union[dict, float]
If by_nuclide is True then a dictionary whose keys are nuclide
names and values are activity is returned. Otherwise the activity
of the material is returned as a float.
If by_nuclide is True then a dictionary whose keys are nuclide names
and values are activity is returned. Otherwise the activity of the
material is returned as a float.
"""
cv.check_value('units', units, {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3', 'Bq/m3', 'Ci', 'Ci/m3'})
@ -1418,6 +1437,9 @@ class Material(IDManagerMixin):
if volume is None:
volume = self.volume
if units in {'Bq', 'Ci'} and volume is None:
raise ValueError(f"Volume must be set in order to compute activity in '{units}'.")
if units == 'Bq':
multiplier = volume
elif units == 'Bq/cm3':
@ -1435,7 +1457,8 @@ class Material(IDManagerMixin):
activity = {}
for nuclide, atoms_per_bcm in self.get_nuclide_atom_densities().items():
inv_seconds = openmc.data.decay_constant(nuclide)
inv_seconds = openmc.data.decay_constant(
nuclide, chain_file=chain_file)
activity[nuclide] = inv_seconds * 1e24 * atoms_per_bcm * multiplier
return activity if by_nuclide else sum(activity.values())
@ -1474,6 +1497,8 @@ class Material(IDManagerMixin):
if units == 'W':
multiplier = volume if volume is not None else self.volume
if multiplier is None:
raise ValueError("Volume must be set in order to compute total decay heat.")
elif units == 'W/cm3':
multiplier = 1
elif units == 'W/m3':
@ -2282,7 +2307,7 @@ class Materials(cv.CheckedList):
multigroup_fluxes: Sequence[Sequence[float]]
Energy-dependent multigroup flux values, where each sublist corresponds
to a specific material. Will be normalized so that it sums to 1.
energy_group_structures': Sequence[Sequence[float] | str]
energy_group_structures: Sequence[Sequence[float] | str]
Energy group boundaries in [eV] or the name of the group structure.
timesteps : iterable of float or iterable of tuple
Array of timesteps. Note that values are not cumulative. The units are
@ -2317,6 +2342,11 @@ class Materials(cv.CheckedList):
for mat in self:
mat.depletable = True
if len(multigroup_fluxes) != len(self):
raise ValueError("multigroup_fluxes length must match number of materials")
if len(energy_group_structures) != len(self):
raise ValueError("energy_group_structures length must match number of materials")
chain = _get_chain(chain_file)
# Create MicroXS objects for all materials
@ -2327,6 +2357,10 @@ class Materials(cv.CheckedList):
for material, flux, energy in zip(
self, multigroup_fluxes, energy_group_structures
):
if material.volume is None:
raise ValueError(
f"Material {material.id} has no volume; cannot deplete"
)
temperature = material.temperature or 293.6
micro_xs = openmc.deplete.MicroXS.from_multigroup_flux(
energies=energy,

234
openmc/mesh.py
View file

@ -3,7 +3,7 @@ import warnings
from abc import ABC, abstractmethod
from collections.abc import Iterable, Sequence, Mapping
from functools import wraps
from math import pi, sqrt, atan2
import math
from numbers import Integral, Real
from pathlib import Path
from typing import Protocol
@ -11,12 +11,10 @@ from typing import Protocol
import h5py
import lxml.etree as ET
import numpy as np
from pathlib import Path
import openmc
import openmc.checkvalue as cv
from openmc.checkvalue import PathLike
from openmc.utility_funcs import change_directory
from .bounding_box import BoundingBox
from ._xml import get_elem_list, get_text
from .mixin import IDManagerMixin
@ -268,10 +266,11 @@ class MeshBase(IDManagerMixin, ABC):
return string
def _volume_dim_check(self):
if self.n_dimension != 3 or \
any([d == 0 for d in self.dimension]):
raise RuntimeError(f'Mesh {self.id} is not 3D. '
'Volumes cannot be provided.')
if any(d == 0 for d in self.dimension):
raise RuntimeError(
f'Mesh {self.id} has a zero-size dimension. '
'Volumes cannot be provided.'
)
@classmethod
def from_hdf5(cls, group: h5py.Group):
@ -290,7 +289,7 @@ class MeshBase(IDManagerMixin, ABC):
"""
mesh_type = 'regular' if 'type' not in group.keys() else group['type'][()].decode()
mesh_id = int(group.name.split('/')[-1].lstrip('mesh '))
mesh_name = '' if not 'name' in group else group['name'][()].decode()
mesh_name = '' if 'name' not in group else group['name'][()].decode()
if mesh_type == 'regular':
return RegularMesh.from_hdf5(group, mesh_id, mesh_name)
@ -1902,7 +1901,7 @@ class CylindricalMesh(StructuredMesh):
self,
r_grid: Sequence[float],
z_grid: Sequence[float],
phi_grid: Sequence[float] = (0, 2*pi),
phi_grid: Sequence[float] = (0, 2*math.pi),
origin: Sequence[float] = (0., 0., 0.),
mesh_id: int | None = None,
name: str = '',
@ -1959,7 +1958,7 @@ class CylindricalMesh(StructuredMesh):
cv.check_length('mesh phi_grid', grid, 2)
cv.check_increasing('mesh phi_grid', grid)
grid = np.asarray(grid, dtype=float)
if np.any((grid < 0.0) | (grid > 2*pi)):
if np.any((grid < 0.0) | (grid > 2*math.pi)):
raise ValueError("phi_grid values must be in [0, 2π].")
self._phi_grid = grid
@ -2027,9 +2026,9 @@ class CylindricalMesh(StructuredMesh):
return string
def get_indices_at_coords(
self,
coords: Sequence[float]
) -> tuple[int, int, int]:
self,
coords: Sequence[float]
) -> tuple[int, int, int]:
"""Finds the index of the mesh element at the specified coordinates.
.. versionadded:: 0.15.0
@ -2045,7 +2044,7 @@ class CylindricalMesh(StructuredMesh):
The r, phi, z indices
"""
r_value_from_origin = sqrt((coords[0]-self.origin[0])**2 + (coords[1]-self.origin[1])**2)
r_value_from_origin = math.hypot(coords[0]-self.origin[0], coords[1]-self.origin[1])
if r_value_from_origin < self.r_grid[0] or r_value_from_origin > self.r_grid[-1]:
raise ValueError(
@ -2069,13 +2068,13 @@ class CylindricalMesh(StructuredMesh):
delta_x = coords[0] - self.origin[0]
delta_y = coords[1] - self.origin[1]
# atan2 returns values in -pi to +pi range
phi_value = atan2(delta_y, delta_x)
phi_value = math.atan2(delta_y, delta_x)
if delta_x < 0 and delta_y < 0:
# returned phi_value anticlockwise and negative
phi_value += 2 * pi
phi_value += 2 * math.pi
if delta_x > 0 and delta_y < 0:
# returned phi_value anticlockwise and negative
phi_value += 2 * pi
phi_value += 2 * math.pi
phi_grid_values = np.array(self.phi_grid)
@ -2110,7 +2109,7 @@ class CylindricalMesh(StructuredMesh):
dimension: Sequence[int] = (10, 10, 10),
mesh_id: int | None = None,
name: str = '',
phi_grid_bounds: Sequence[float] = (0.0, 2*pi),
phi_grid_bounds: Sequence[float] = (0.0, 2*math.pi),
enclose_domain: bool = False,
) -> CylindricalMesh:
"""Create CylindricalMesh from a bounding box.
@ -2338,8 +2337,8 @@ class SphericalMesh(StructuredMesh):
def __init__(
self,
r_grid: Sequence[float],
phi_grid: Sequence[float] = (0, 2*pi),
theta_grid: Sequence[float] = (0, pi),
phi_grid: Sequence[float] = (0, 2*math.pi),
theta_grid: Sequence[float] = (0, math.pi),
origin: Sequence[float] = (0., 0., 0.),
mesh_id: int | None = None,
name: str = '',
@ -2396,7 +2395,7 @@ class SphericalMesh(StructuredMesh):
cv.check_length('mesh theta_grid', grid, 2)
cv.check_increasing('mesh theta_grid', grid)
grid = np.asarray(grid, dtype=float)
if np.any((grid < 0.0) | (grid > pi)):
if np.any((grid < 0.0) | (grid > math.pi)):
raise ValueError("theta_grid values must be in [0, π].")
self._theta_grid = grid
@ -2410,7 +2409,7 @@ class SphericalMesh(StructuredMesh):
cv.check_length('mesh phi_grid', grid, 2)
cv.check_increasing('mesh phi_grid', grid)
grid = np.asarray(grid, dtype=float)
if np.any((grid < 0.0) | (grid > 2*pi)):
if np.any((grid < 0.0) | (grid > 2*math.pi)):
raise ValueError("phi_grid values must be in [0, 2π].")
self._phi_grid = grid
@ -2482,8 +2481,8 @@ class SphericalMesh(StructuredMesh):
dimension: Sequence[int] = (10, 10, 10),
mesh_id: int | None = None,
name: str = '',
phi_grid_bounds: Sequence[float] = (0.0, 2*pi),
theta_grid_bounds: Sequence[float] = (0.0, pi),
phi_grid_bounds: Sequence[float] = (0.0, 2*math.pi),
theta_grid_bounds: Sequence[float] = (0.0, math.pi),
enclose_domain: bool = False,
) -> SphericalMesh:
"""Create SphericalMesh from a bounding box.
@ -2644,10 +2643,82 @@ class SphericalMesh(StructuredMesh):
arr[..., 2] = z + origin[2]
return arr
def get_indices_at_coords(self, coords: Sequence[float]) -> tuple:
raise NotImplementedError(
"get_indices_at_coords is not yet implemented for SphericalMesh"
)
def get_indices_at_coords(
self,
coords: Sequence[float]
) -> tuple[int, int, int]:
"""Find the mesh cell indices containing the specified coordinates.
.. versionadded:: 0.15.4
Parameters
----------
coords : Sequence[float]
Cartesian coordinates of the point as (x, y, z).
Returns
-------
tuple[int, int, int]
The r, theta, phi indices.
Raises
------
ValueError
If the coordinates fall outside the mesh grid boundaries.
"""
dx = coords[0] - self.origin[0]
dy = coords[1] - self.origin[1]
dz = coords[2] - self.origin[2]
r_value = math.hypot(dx, dy, dz)
if r_value < self.r_grid[0] or r_value > self.r_grid[-1]:
raise ValueError(
f'The r value {r_value} computed from the specified '
f'coordinates is outside the r grid range '
f'[{self.r_grid[0]}, {self.r_grid[-1]}].'
)
r_index = int(min(
np.searchsorted(self.r_grid, r_value, side='right') - 1,
len(self.r_grid) - 2
))
if r_value == 0.0:
theta_value = 0.0
phi_value = 0.0
else:
theta_value = math.acos(dz / r_value)
phi_value = math.atan2(dy, dx)
if phi_value < 0:
phi_value += 2 * math.pi
if theta_value < self.theta_grid[0] or theta_value > self.theta_grid[-1]:
raise ValueError(
f'The theta value {theta_value} computed from the specified '
f'coordinates is outside the theta grid range '
f'[{self.theta_grid[0]}, {self.theta_grid[-1]}].'
)
theta_index = int(min(
np.searchsorted(self.theta_grid, theta_value, side='right') - 1,
len(self.theta_grid) - 2
))
if phi_value < self.phi_grid[0] or phi_value > self.phi_grid[-1]:
raise ValueError(
f'The phi value {phi_value} computed from the specified '
f'coordinates is outside the phi grid range '
f'[{self.phi_grid[0]}, {self.phi_grid[-1]}].'
)
phi_index = int(min(
np.searchsorted(self.phi_grid, phi_value, side='right') - 1,
len(self.phi_grid) - 2
))
return (r_index, theta_index, phi_index)
def require_statepoint_data(func):
@ -2737,7 +2808,8 @@ class UnstructuredMesh(MeshBase):
_UNSUPPORTED_ELEM = -1
_LINEAR_TET = 0
_LINEAR_HEX = 1
_VTK_TETRA = 10
_VTK_TET = 10
_VTK_HEX = 12
def __init__(self, filename: PathLike, library: str, mesh_id: int | None = None,
name: str = '', length_multiplier: float = 1.0,
@ -3048,7 +3120,9 @@ class UnstructuredMesh(MeshBase):
n_skipped += 1
continue
else:
raise RuntimeError(f"Invalid element type {elem_type} found")
raise RuntimeError(
f"Invalid element type {elem_type} found in mesh {self.id}"
)
for i, c in enumerate(conn):
if c == -1:
@ -3105,35 +3179,42 @@ class UnstructuredMesh(MeshBase):
datasets: dict | None = None,
volume_normalization: bool = True,
):
def append_dataset(dset, array):
"""Convenience function to append data to an HDF5 dataset"""
origLen = dset.shape[0]
dset.resize(origLen + array.shape[0], axis=0)
dset[origLen:] = array
# This writer supports linear tetrahedra and linear hexahedra elements
conn_list = [] # flattened connectivity ids
cell_sizes = [] # number of points per cell
vtk_types = [] # VTK cell types per cell (uint8)
n_skipped = 0
if self.library != "moab":
raise NotImplementedError("VTKHDF output is only supported for MOAB meshes")
# the self.connectivity contains arrays of length 8 to support hex
# elements as well, in the case of tetrahedra mesh elements, the
# last 4 values are -1 and are removed
trimmed_connectivity = []
for cell in self.connectivity:
# Find the index of the first -1 value, if any
first_negative_index = np.where(cell == -1)[0]
if first_negative_index.size > 0:
# Slice the array up to the first -1 value
trimmed_connectivity.append(cell[: first_negative_index[0]])
for conn, etype in zip(self.connectivity, self.element_types):
if etype == self._LINEAR_TET:
ids = conn[:4]
vtk_types.append(self._VTK_TET)
elif etype == self._LINEAR_HEX:
ids = conn[:8]
vtk_types.append(self._VTK_HEX)
elif etype == self._UNSUPPORTED_ELEM:
n_skipped += 1
continue
else:
# No -1 values, append the whole cell
trimmed_connectivity.append(cell)
trimmed_connectivity = np.array(trimmed_connectivity, dtype="int32").flatten()
raise RuntimeError(
f"Invalid element type {etype} found in mesh {self.id}"
)
conn_list.extend(ids.tolist())
cell_sizes.append(len(ids))
# MOAB meshes supports tet elements only so we know it has 4 points per cell
points_per_cell = 4
if n_skipped > 0:
warnings.warn(
f"{n_skipped} elements were not written because "
"they are not of type linear tet/hex"
)
# offsets are the indices of the first point of each cell in the array of points
offsets = np.arange(0, self.n_elements * points_per_cell + 1, points_per_cell)
connectivity = np.asarray(conn_list, dtype=np.int64)
# Offsets must be length (numCells + 1) with a leading 0 and
# cumulative end-indices thereafter, per VTK's layout
cell_sizes_arr = np.asarray(cell_sizes, dtype=np.int64)
offsets = np.zeros(cell_sizes_arr.size + 1, dtype=np.int64)
np.cumsum(cell_sizes_arr, out=offsets[1:])
for name, data in datasets.items():
if data.shape != self.dimension:
@ -3155,42 +3236,27 @@ class UnstructuredMesh(MeshBase):
dtype=h5py.string_dtype("ascii", len(ascii_type)),
)
# create hdf5 file structure
root.create_dataset("NumberOfPoints", (0,), maxshape=(None,), dtype="i8")
root.create_dataset("Types", (0,), maxshape=(None,), dtype="uint8")
root.create_dataset("Points", (0, 3), maxshape=(None, 3), dtype="f")
root.create_dataset(
"NumberOfConnectivityIds", (0,), maxshape=(None,), dtype="i8"
)
root.create_dataset("NumberOfCells", (0,), maxshape=(None,), dtype="i8")
root.create_dataset("Offsets", (0,), maxshape=(None,), dtype="i8")
root.create_dataset("Connectivity", (0,), maxshape=(None,), dtype="i8")
# Create HDF5 file structure compliant with VTKHDF UnstructuredGrid
n_points = int(len(self.vertices))
n_cells = int(len(cell_sizes))
n_conn_ids = int(len(connectivity))
append_dataset(root["NumberOfPoints"], np.array([len(self.vertices)]))
append_dataset(root["Points"], self.vertices)
append_dataset(
root["NumberOfConnectivityIds"],
np.array([len(trimmed_connectivity)]),
)
append_dataset(root["Connectivity"], trimmed_connectivity)
append_dataset(root["NumberOfCells"], np.array([self.n_elements]))
append_dataset(root["Offsets"], offsets)
append_dataset(
root["Types"], np.full(self.n_elements, self._VTK_TETRA, dtype="uint8")
)
root.create_dataset("NumberOfPoints", data=(n_points,), dtype="i8")
root.create_dataset("NumberOfCells", data=(n_cells,), dtype="i8")
root.create_dataset("NumberOfConnectivityIds", data=(n_conn_ids,), dtype="i8")
root.create_dataset("Points", data=self.vertices.astype(np.float64, copy=False), dtype="f8")
root.create_dataset("Types", data=np.asarray(vtk_types, dtype=np.uint8), dtype="uint8")
root.create_dataset("Offsets", data=offsets.astype("i8"), dtype="i8")
root.create_dataset("Connectivity", data=connectivity.astype("i8"), dtype="i8")
cell_data_group = root.create_group("CellData")
for name, data in datasets.items():
cell_data_group.create_dataset(
name, (0,), maxshape=(None,), dtype="float64", chunks=True
)
if volume_normalization:
data /= self.volumes
append_dataset(cell_data_group[name], data)
cell_data_group.create_dataset(
name, data=data, dtype="float64", chunks=True
)
@classmethod
def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str):

8
openmc/mgxs/library.py
View file

@ -1290,6 +1290,14 @@ class Library:
'are ignored since multiplicity or nu-scatter matrices '\
'were not tallied for ' + xsdata_name
warn(msg, RuntimeWarning)
if 'scatter matrix' in self.mgxs_types:
scatt_mgxs = self.get_mgxs(domain, 'scatter matrix')
elif 'consistent scatter matrix' in self.mgxs_types:
scatt_mgxs = self.get_mgxs(domain, 'consistent scatter matrix')
else:
raise ValueError(f'No scatter matrix found for {xsdata_name}.')
xsdata.set_scatter_matrix_mgxs(scatt_mgxs, temperature=temperature,
xs_type=xs_type,
nuclide=[nuclide],

238
openmc/model/model.py
View file

@ -21,7 +21,8 @@ import openmc
import openmc._xml as xml
from openmc.dummy_comm import DummyCommunicator
from openmc.executor import _process_CLI_arguments
from openmc.checkvalue import check_type, check_value, PathLike
from openmc.checkvalue import (check_type, check_value, check_greater_than,
check_length, PathLike)
from openmc.exceptions import InvalidIDError
from openmc.plots import add_plot_params, _BASIS_INDICES, id_map_to_rgb
from openmc.utility_funcs import change_directory
@ -33,6 +34,15 @@ class ModelModifier(Protocol):
...
def _check_pixels(pixels: int | Sequence[int]) -> None:
if isinstance(pixels, Integral):
check_greater_than('pixels', pixels, 0)
else:
check_length('pixels', pixels, 2)
for p in pixels:
check_greater_than('pixels', p, 0)
class Model:
"""Model container.
@ -265,22 +275,22 @@ class Model:
denom_tally = openmc.Tally(name='IFP denominator')
denom_tally.scores = ['ifp-denominator']
self.tallies.append(denom_tally)
# TODO: This should also be incorporated into lower-level calls in
# TODO: This should also be incorporated into lower-level calls in
# settings.py, but it requires information about the tallies currently
# on the active Model
def _assign_fw_cadis_tally_IDs(self):
# Verify that all tallies assigned as targets on WeightWindowGenerators
# exist within model.tallies. If this is the case, convert the .targets
# Verify that all tallies assigned as targets on WeightWindowGenerators
# exist within model.tallies. If this is the case, convert the .targets
# attribute of each WeightWindowGenerator to a sequence of tally IDs.
if len(self.settings.weight_window_generators) == 0:
return
# List of valid tally IDs
reference_tally_ids = np.asarray([tal.id for tal in self.tallies])
for wwg in self.settings.weight_window_generators:
# Only proceeds if the "targets" attribute is an openmc.Tallies,
# Only proceeds if the "targets" attribute is an openmc.Tallies,
# which means it hasn't been checked against model.tallies.
if isinstance(wwg.targets, openmc.Tallies):
id_vec = []
@ -291,8 +301,8 @@ class Model:
if tal == reference_tal:
id_next = reference_tal.id
break
if id_next == None:
if id_next is None:
raise RuntimeError(
f'Local FW-CADIS target tally {tal.id} not found on model.tallies!')
else:
@ -467,6 +477,8 @@ class Model:
This method iterates over all DAGMC universes in the geometry and
synchronizes their cells with the current material assignments. Requires
that the model has been initialized via :meth:`Model.init_lib`.
Synchronized DAGMC cells can then be edited and exported as nested
`<cell>` overrides inside each `<dagmc_universe>` element.
.. versionadded:: 0.15.1
@ -1048,6 +1060,8 @@ class Model:
pixels: int | Sequence[int],
basis: str
):
_check_pixels(pixels)
x, y, _ = _BASIS_INDICES[basis]
bb = self.bounding_box
@ -1125,28 +1139,150 @@ class Model:
array contains cell IDs, cell instances, and material IDs (in that
order).
"""
ids, _ = self.slice_data(
origin=origin,
width=width,
pixels=pixels,
basis=basis,
show_overlaps=color_overlaps,
level=-1,
include_properties=False,
**init_kwargs,
)
return ids
def slice_data(
self,
origin: Sequence[float] | None = None,
width: Sequence[float] | None = None,
pixels: int | Sequence[int] = 40000,
basis: str = 'xy',
u_span: Sequence[float] | None = None,
v_span: Sequence[float] | None = None,
show_overlaps: bool = False,
level: int = -1,
filter: openmc.Filter | None = None,
include_properties: bool = True,
**init_kwargs
) -> tuple[np.ndarray, np.ndarray | None]:
"""Generate geometry and property data for a 2D plot slice.
This method combines the functionality of :meth:`id_map` and property
mapping into a single call, avoiding duplicate geometry lookups. It also
supports filter bin index lookup for tally visualization.
.. versionadded:: 0.16.0
Parameters
----------
origin : Sequence[float], optional
Origin of the plot. If unspecified, this argument defaults to the
center of the bounding box if the bounding box does not contain inf
values for the provided basis, otherwise (0.0, 0.0, 0.0).
width : Sequence[float], optional
Width of the plot. If unspecified, this argument defaults to the
width of the bounding box if the bounding box does not contain inf
values for the provided basis, otherwise (10.0, 10.0).
pixels : int | Sequence[int], optional
If an iterable of ints is provided then this directly sets the
number of pixels to use in each basis direction. If a single int is
provided then this sets the total number of pixels in the plot and
the number of pixels in each basis direction is calculated from this
total and the image aspect ratio based on the width argument.
basis : {'xy', 'yz', 'xz'}, optional
Basis of the plot.
u_span : Sequence[float], optional
Full-width span vector for an oriented slice (3 values). Mutually
exclusive with width.
v_span : Sequence[float], optional
Full-height span vector for an oriented slice (3 values). Mutually
exclusive with width.
show_overlaps : bool, optional
Whether to identify and assign unique IDs (-3) to overlapping
regions. If False, overlapping regions will be assigned the ID of
the lowest-numbered cell that occupies that region. Defaults to
False.
level : int, optional
Universe level to plot (-1 for deepest). Defaults to -1.
filter : openmc.Filter, optional
If provided, the information for each pixel also includes an index
in the filter corresponding to the pixel position.
include_properties : bool, optional
Whether to include temperature/density data. Defaults to True.
**init_kwargs
Keyword arguments passed to :meth:`Model.init_lib`.
Returns
-------
geom_data : numpy.ndarray
Shape (v_res, h_res, 3) or (v_res, h_res, 4) int32 array. Contains
[cell_id, cell_instance, material_id] when no filter, or [cell_id,
cell_instance, material_id, filter_bin] with filter.
property_data : numpy.ndarray or None
Shape (v_res, h_res, 2) float64 array with [temperature, density],
or None if include_properties=False.
"""
import openmc.lib
origin, width, pixels = self._set_plot_defaults(
origin, width, pixels, basis)
_check_pixels(pixels)
# initialize the openmc.lib.plot._PlotBase object
plot_obj = openmc.lib.plot._PlotBase()
plot_obj.origin = origin
plot_obj.width = width[0]
plot_obj.height = width[1]
plot_obj.h_res = pixels[0]
plot_obj.v_res = pixels[1]
plot_obj.basis = basis
plot_obj.color_overlaps = color_overlaps
if width is not None and (u_span is not None or v_span is not None):
raise ValueError("width is mutually exclusive with u_span/v_span.")
if u_span is not None or v_span is not None:
if u_span is None or v_span is None:
raise ValueError("Both u_span and v_span must be provided.")
if origin is None:
origin = (0.0, 0.0, 0.0)
if isinstance(pixels, int):
u_norm = np.linalg.norm(u_span)
v_norm = np.linalg.norm(v_span)
aspect_ratio = u_norm / v_norm
pixels_y = math.sqrt(pixels / aspect_ratio)
pixels = (int(pixels / pixels_y), int(pixels_y))
else:
origin, width, pixels = self._set_plot_defaults(
origin, width, pixels, basis)
# Silence output by default. Also set arguments to start in volume
# calculation mode to avoid loading cross sections
init_kwargs.setdefault('output', False)
init_kwargs.setdefault('args', ['-c'])
# If filter does not already appear in the model, temporarily add a
# tally with the filter
original_length = len(self.tallies)
if filter is not None:
filter_ids = {f.id for t in self.tallies for f in t.filters}
if filter.id not in filter_ids:
# Create temporary tally while preserving ID assignment
next_id = openmc.Tally.next_id
temp_tally = openmc.Tally()
temp_tally.filters = [filter]
temp_tally.scores = ['flux']
self.tallies.append(temp_tally)
openmc.Tally.used_ids.remove(temp_tally.id)
openmc.Tally.next_id = next_id
with openmc.lib.TemporarySession(self, **init_kwargs):
return openmc.lib.id_map(plot_obj)
geom_data, property_data = openmc.lib.slice_data(
origin=origin,
width=width,
basis=basis,
u_span=u_span,
v_span=v_span,
pixels=pixels,
show_overlaps=show_overlaps,
level=level,
filter=filter,
include_properties=include_properties,
)
# If filter was temporarily added, remove it
if len(self.tallies) > original_length:
self.tallies.pop()
return geom_data, property_data
@add_plot_params
def plot(
@ -1179,7 +1315,11 @@ class Model:
import matplotlib.pyplot as plt
check_type('n_samples', n_samples, int | None)
if n_samples is not None:
check_greater_than('n_samples', n_samples, 0, equality=True)
check_type('plane_tolerance', plane_tolerance, Real)
check_greater_than('plane_tolerance', plane_tolerance, 0.0)
if legend_kwargs is None:
legend_kwargs = {}
legend_kwargs.setdefault('bbox_to_anchor', (1.05, 1))
@ -1214,13 +1354,14 @@ class Model:
"openmc.config before plotting.")
break
# Get ID map from the C API
id_map = self.id_map(
# Get plot IDs from the C API
id_map, _ = self.slice_data(
origin=origin,
width=width,
pixels=pixels,
basis=basis,
color_overlaps=show_overlaps
show_overlaps=show_overlaps,
include_properties=False,
)
# Generate colors if not provided
@ -1746,10 +1887,10 @@ class Model:
@staticmethod
def _auto_generate_mgxs_lib(
model: openmc.model.model,
model: openmc.model.Model,
groups: openmc.mgxs.EnergyGroups,
correction: str | none,
directory: pathlike,
correction: str | None,
directory: PathLike,
) -> openmc.mgxs.Library:
"""
Automatically generate a multi-group cross section libray from a model
@ -1956,7 +2097,7 @@ class Model:
# Set materials on the model
model.materials = [material]
if temperature != None:
if temperature is not None:
model.materials[-1].temperature = temperature
# Settings
@ -1983,7 +2124,7 @@ class Model:
mgxs_lib = Model._auto_generate_mgxs_lib(
model, groups, correction, directory)
if temperature != None:
if temperature is not None:
return mgxs_lib.get_xsdata(domain=material, xsdata_name=name,
temperature=temperature)
else:
@ -2056,12 +2197,12 @@ class Model:
)
temp_settings = {}
if temperature_settings == None:
if temperature_settings is None:
temp_settings = self.settings.temperature
else:
temp_settings = temperature_settings
if temperatures == None:
if temperatures is None:
mgxs_sets = []
for material in self.materials:
xs_data = Model._isothermal_infinite_media_mgxs(
@ -2234,7 +2375,7 @@ class Model:
model = openmc.Model()
model.geometry = stoch_geom
if temperature != None:
if temperature is not None:
for material in model.geometry.get_all_materials().values():
material.temperature = temperature
@ -2258,7 +2399,7 @@ class Model:
model, groups, correction, directory)
# Fetch all of the isothermal results.
if temperature != None:
if temperature is not None:
return {
mat.name : mgxs_lib.get_xsdata(domain=mat, xsdata_name=mat.name,
temperature=temperature)
@ -2344,12 +2485,12 @@ class Model:
)
temp_settings = {}
if temperature_settings == None:
if temperature_settings is None:
temp_settings = self.settings.temperature
else:
temp_settings = temperature_settings
if temperatures == None:
if temperatures is None:
mgxs_sets = Model._isothermal_stochastic_slab_mgxs(
geo,
groups,
@ -2442,7 +2583,7 @@ class Model:
model = copy.deepcopy(input_model)
model.tallies = openmc.Tallies()
if temperature != None:
if temperature is not None:
for material in model.geometry.get_all_materials().values():
material.temperature = temperature
@ -2458,7 +2599,7 @@ class Model:
model, groups, correction, directory)
# Fetch all of the isothermal results.
if temperature != None:
if temperature is not None:
return {
mat.name : mgxs_lib.get_xsdata(domain=mat, xsdata_name=mat.name,
temperature=temperature)
@ -2513,12 +2654,12 @@ class Model:
entries in openmc.Settings.temperature_settings.
"""
temp_settings = {}
if temperature_settings == None:
if temperature_settings is None:
temp_settings = self.settings.temperature
else:
temp_settings = temperature_settings
if temperatures == None:
if temperatures is None:
mgxs_sets = Model._isothermal_materialwise_mgxs(
self,
groups,
@ -2649,12 +2790,15 @@ class Model:
self.settings.run_mode = original_run_mode
break
# Make sure all materials have a name, and that the name is a valid HDF5
# dataset name
# Temporarily replace each material's name with a unique, valid HDF5
# dataset name (its name plus ID) for use as its MGXS library entry
# and macroscopic. The ID keeps the name unique even when materials
# share a name; the original names are restored at the end.
original_names = [material.name for material in self.materials]
for material in self.materials:
if not material.name or not material.name.strip():
material.name = f"material {material.id}"
material.name = re.sub(r'[^a-zA-Z0-9]', '_', material.name)
base = material.name if material.name and material.name.strip() \
else "material"
material.name = re.sub(r'[^a-zA-Z0-9]', '_', base) + f"_{material.id}"
# If needed, generate the needed MGXS data library file
if not Path(mgxs_path).is_file() or overwrite_mgxs_library:
@ -2686,6 +2830,10 @@ class Model:
self.settings.energy_mode = 'multi-group'
# Restore the user's original material names.
for material, name in zip(self.materials, original_names):
material.name = name
def convert_to_random_ray(self):
"""Convert a multigroup model to use random ray.

14
openmc/plots.py
View file

@ -383,7 +383,7 @@ def id_map_to_rgb(
"""
# Initialize RGB array with white background (values between 0 and 1 for matplotlib)
img = np.ones(id_map.shape, dtype=float)
# Get the appropriate index based on color_by
if color_by == 'cell':
id_index = 0 # Cell IDs are in the first channel
@ -391,14 +391,14 @@ def id_map_to_rgb(
id_index = 2 # Material IDs are in the third channel
else:
raise ValueError("color_by must be either 'cell' or 'material'")
# Get all unique IDs in the plot
unique_ids = np.unique(id_map[:, :, id_index])
# Generate default colors if not provided
if colors is None:
colors = {}
# Convert colors dict to use IDs as keys
color_map = {}
for key, color in colors.items():
@ -406,13 +406,13 @@ def id_map_to_rgb(
color_map[key.id] = color
else:
color_map[key] = color
# Generate random colors for IDs not in color_map
rng = np.random.RandomState(1)
for uid in unique_ids:
if uid > 0 and uid not in color_map:
color_map[uid] = rng.randint(0, 256, (3,))
# Apply colors to each pixel
for uid in unique_ids:
if uid == -1: # Background/void
@ -432,7 +432,7 @@ def id_map_to_rgb(
rgb = color
mask = id_map[:, :, id_index] == uid
img[mask] = np.array(rgb) / 255.0
return img
class PlotBase(IDManagerMixin):

43
openmc/settings.py
View file

@ -67,6 +67,10 @@ class Settings:
(ex: ["(n,fission)", 2, "(n,2n)"] ). (list of str or int)
:deposited_E_threshold: Number to define the minimum deposited energy during
per collision to trigger banking. (float)
create_delayed_neutrons : bool
Whether delayed neutrons are created in fission.
.. versionadded:: 0.13.3
create_fission_neutrons : bool
Indicate whether fission neutrons should be created or not.
cutoff : dict
@ -256,8 +260,15 @@ class Settings:
The type of calculation to perform (default is 'eigenvalue')
seed : int
Seed for the linear congruential pseudorandom number generator
stride : int
Number of random numbers allocated for each source particle history
shared_secondary_bank : bool
Whether to use a shared secondary particle bank. When enabled,
secondary particles are collected into a global bank, sorted for
reproducibility, and load-balanced across MPI ranks between
generations. If not specified, the shared secondary bank is
enabled automatically for fixed-source simulations with weight
windows active, and disabled otherwise.
.. versionadded:: 0.15.4
source : Iterable of openmc.SourceBase
Distribution of source sites in space, angle, and energy
source_rejection_fraction : float
@ -277,6 +288,8 @@ class Settings:
Options for writing state points. Acceptable keys are:
:batches: list of batches at which to write statepoint files
stride : int
Number of random numbers allocated for each source particle history
surf_source_read : dict
Options for reading surface source points. Acceptable keys are:
@ -373,10 +386,6 @@ class Settings:
.. versionadded:: 0.14.0
create_delayed_neutrons : bool
Whether delayed neutrons are created in fission.
.. versionadded:: 0.13.3
weight_windows_on : bool
Whether weight windows are enabled
@ -482,6 +491,7 @@ class Settings:
self._weight_window_generators = cv.CheckedList(
WeightWindowGenerator, 'weight window generators')
self._weight_windows_on = None
self._shared_secondary_bank = None
self._weight_windows_file = None
self._weight_window_checkpoints = {}
self._max_history_splits = None
@ -1306,6 +1316,15 @@ class Settings:
cv.check_type('weight windows on', value, bool)
self._weight_windows_on = value
@property
def shared_secondary_bank(self) -> bool:
return self._shared_secondary_bank
@shared_secondary_bank.setter
def shared_secondary_bank(self, value: bool):
cv.check_type('shared secondary bank', value, bool)
self._shared_secondary_bank = value
@property
def weight_window_checkpoints(self) -> dict:
return self._weight_window_checkpoints
@ -1924,6 +1943,11 @@ class Settings:
elem = ET.SubElement(root, "weight_windows_on")
elem.text = str(self._weight_windows_on).lower()
def _create_shared_secondary_bank_subelement(self, root):
if self._shared_secondary_bank is not None:
elem = ET.SubElement(root, "shared_secondary_bank")
elem.text = str(self._shared_secondary_bank).lower()
def _create_weight_window_generators_subelement(self, root, mesh_memo=None):
if not self.weight_window_generators:
return
@ -2430,6 +2454,11 @@ class Settings:
if text is not None:
self.weight_windows_on = text in ('true', '1')
def _shared_secondary_bank_from_xml_element(self, root):
text = get_text(root, 'shared_secondary_bank')
if text is not None:
self.shared_secondary_bank = text in ('true', '1')
def _weight_windows_file_from_xml_element(self, root):
text = get_text(root, 'weight_windows_file')
if text is not None:
@ -2597,6 +2626,7 @@ class Settings:
self._create_write_initial_source_subelement(element)
self._create_weight_windows_subelement(element, mesh_memo)
self._create_weight_windows_on_subelement(element)
self._create_shared_secondary_bank_subelement(element)
self._create_weight_window_generators_subelement(element, mesh_memo)
self._create_weight_windows_file_element(element)
self._create_weight_window_checkpoints_subelement(element)
@ -2714,6 +2744,7 @@ class Settings:
settings._write_initial_source_from_xml_element(elem)
settings._weight_windows_from_xml_element(elem, meshes)
settings._weight_windows_on_from_xml_element(elem)
settings._shared_secondary_bank_from_xml_element(elem)
settings._weight_windows_file_from_xml_element(elem)
settings._weight_window_generators_from_xml_element(elem, meshes)
settings._weight_window_checkpoints_from_xml_element(elem)

430
openmc/source.py
View file

@ -1,7 +1,7 @@
from __future__ import annotations
from abc import ABC, abstractmethod
from collections.abc import Iterable, Sequence
from numbers import Real
from numbers import Integral, Real
from pathlib import Path
import warnings
from typing import Any
@ -46,7 +46,7 @@ class SourceBase(ABC):
Attributes
----------
type : {'independent', 'file', 'compiled', 'mesh'}
type : {'independent', 'file', 'compiled', 'mesh', 'tokamak'}
Indicator of source type.
strength : float
Strength of the source
@ -203,6 +203,8 @@ class SourceBase(ABC):
return FileSource.from_xml_element(elem)
elif source_type == 'mesh':
return MeshSource.from_xml_element(elem, meshes)
elif source_type == 'tokamak':
return TokamakSource.from_xml_element(elem)
else:
raise ValueError(
f'Source type {source_type} is not recognized')
@ -896,6 +898,430 @@ class FileSource(SourceBase):
return cls(**kwargs)
class TokamakSource(SourceBase):
r"""A source representing neutron emission from a tokamak plasma.
This source samples neutron positions from a tokamak plasma geometry using
Miller-style flux surface parameterization. The user provides an emission
profile S(r/a) as a function of normalized minor radius, along with one or
more energy distributions.
The flux surface parameterization is
.. math::
\begin{aligned}
R &= R_0 + r \cos\left(\alpha + \delta \sin\alpha\right)
+ \Delta \left[1 - \left(\frac{r}{a}\right)^2\right] \\
Z &= Z_\mathrm{shift} + \kappa r \sin\alpha
\end{aligned}
where :math:`R_0` is major radius, :math:`a` is minor radius,
:math:`\kappa` is elongation, :math:`\delta` is triangularity,
:math:`\Delta` is the Shafranov shift, and :math:`Z_\mathrm{shift}` is
the vertical shift.
.. versionadded:: 0.15.4
Parameters
----------
major_radius : float
Major radius R0 in [cm]
minor_radius : float
Minor radius a in [cm]
elongation : float
Plasma elongation κ (must be > 0)
triangularity : float
Plasma triangularity δ (must be in [-1, 1])
shafranov_shift : float
Shafranov shift Δ in [cm] (must be >= 0 and < a/2)
r_over_a : numpy.ndarray
Normalized minor radius grid points, must start at 0 and end at 1
emission_density : numpy.ndarray
Emission density S(r) at each r/a point (arbitrary units, must be >= 0).
Values are linearly interpolated between grid points and refined on an
internal grid for radial sampling. Must have the same length as
``r_over_a`` and contain at least one positive value.
energy : openmc.stats.Univariate or Sequence[openmc.stats.Univariate]
Energy distribution(s). Either a single distribution used at all radii,
or one distribution per ``r_over_a`` grid point. When one distribution
per grid point is given, the energy of a sampled particle is drawn from
one of the two distributions bracketing its sampled radius, selected
stochastically with probability proportional to the proximity of the
radius to each grid point (stochastic interpolation).
time : openmc.stats.Univariate, optional
Time distribution of the source. If None, particles are born at
:math:`t=0`, matching the default behavior of
:class:`openmc.IndependentSource`.
phi_start : float
Starting toroidal angle in [rad] (default: 0)
phi_extent : float
Toroidal angle extent in [rad] (default: 2π)
n_alpha : int
Number of poloidal angle grid points for CDF sampling (default: 101)
vertical_shift : float
Vertical shift of the plasma center in [cm] (default: 0)
strength : float
Strength of the source (default: 1.0)
constraints : dict
Constraints on sampled source particles. See :class:`SourceBase` for
valid keys and values.
Attributes
----------
major_radius : float
Major radius R0 in [cm]
minor_radius : float
Minor radius a in [cm]
elongation : float
Plasma elongation κ
triangularity : float
Plasma triangularity δ
shafranov_shift : float
Shafranov shift Δ in [cm]
r_over_a : numpy.ndarray
Normalized minor radius grid points
emission_density : numpy.ndarray
Emission density S(r) at each r/a point
energy : list of openmc.stats.Univariate
Energy distribution(s)
time : openmc.stats.Univariate or None
Time distribution of the source
phi_start : float
Starting toroidal angle in [rad]
phi_extent : float
Toroidal angle extent in [rad]
n_alpha : int
Number of poloidal angle grid points
vertical_shift : float
Vertical shift of the plasma center in [cm]
strength : float
Strength of the source
type : str
Indicator of source type: 'tokamak'
constraints : dict
Constraints on sampled source particles
"""
def __init__(
self,
major_radius: float,
minor_radius: float,
elongation: float,
triangularity: float,
shafranov_shift: float,
r_over_a: Sequence[float],
emission_density: Sequence[float],
energy: Univariate | Sequence[Univariate],
time: Univariate | None = None,
phi_start: float = 0.0,
phi_extent: float = 2.0 * np.pi,
n_alpha: int = 101,
vertical_shift: float = 0.0,
strength: float = 1.0,
constraints: dict[str, Any] | None = None
):
super().__init__(strength=strength, constraints=constraints)
self.major_radius = major_radius
self.minor_radius = minor_radius
self.elongation = elongation
self.triangularity = triangularity
self.shafranov_shift = shafranov_shift
self.r_over_a = r_over_a
self.emission_density = emission_density
self.phi_start = phi_start
self.phi_extent = phi_extent
self.n_alpha = n_alpha
self.vertical_shift = vertical_shift
self.energy = energy
self.time = time
self._validate()
def _validate(self):
"""Validate relationships between tokamak source parameters."""
if self.minor_radius >= self.major_radius:
raise ValueError(
f"minor_radius ({self.minor_radius}) must be smaller than "
f"major_radius ({self.major_radius})")
if self.shafranov_shift >= 0.5 * self.minor_radius:
raise ValueError(
f"shafranov_shift ({self.shafranov_shift}) must be smaller "
f"than half the minor_radius ({0.5 * self.minor_radius})")
if len(self.emission_density) != len(self.r_over_a):
raise ValueError(
f"emission_density (length {len(self.emission_density)}) must "
f"have the same length as r_over_a (length {len(self.r_over_a)})")
if not np.any(self.emission_density > 0.0):
raise ValueError("emission_density must contain a positive value")
if len(self.energy) not in (1, len(self.r_over_a)):
raise ValueError(
f"Number of energy distributions ({len(self.energy)}) must be "
f"either 1 or equal to the number of r_over_a grid points "
f"({len(self.r_over_a)})")
@property
def type(self) -> str:
return "tokamak"
@property
def major_radius(self) -> float:
return self._major_radius
@major_radius.setter
def major_radius(self, value: float):
cv.check_type('major radius', value, Real)
cv.check_greater_than('major radius', value, 0.0)
self._major_radius = value
@property
def minor_radius(self) -> float:
return self._minor_radius
@minor_radius.setter
def minor_radius(self, value: float):
cv.check_type('minor radius', value, Real)
cv.check_greater_than('minor radius', value, 0.0)
self._minor_radius = value
@property
def elongation(self) -> float:
return self._elongation
@elongation.setter
def elongation(self, value: float):
cv.check_type('elongation', value, Real)
cv.check_greater_than('elongation', value, 0.0)
self._elongation = value
@property
def triangularity(self) -> float:
return self._triangularity
@triangularity.setter
def triangularity(self, value: float):
cv.check_type('triangularity', value, Real)
cv.check_greater_than('triangularity', value, -1.0, equality=True)
cv.check_less_than('triangularity', value, 1.0, equality=True)
self._triangularity = value
@property
def shafranov_shift(self) -> float:
return self._shafranov_shift
@shafranov_shift.setter
def shafranov_shift(self, value: float):
cv.check_type('Shafranov shift', value, Real)
cv.check_greater_than('Shafranov shift', value, 0.0, equality=True)
self._shafranov_shift = value
@property
def r_over_a(self) -> np.ndarray:
return self._r_over_a
@r_over_a.setter
def r_over_a(self, value: Sequence[float]):
value = np.asarray(value, dtype=float)
if value.ndim != 1 or len(value) < 2:
raise ValueError("r_over_a must be a 1-D array with at least 2 points")
if value[0] != 0.0:
raise ValueError("r_over_a must start at 0")
if value[-1] != 1.0:
raise ValueError("r_over_a must end at 1")
if not np.all(np.diff(value) > 0):
raise ValueError("r_over_a must be strictly increasing")
self._r_over_a = value
@property
def emission_density(self) -> np.ndarray:
return self._emission_density
@emission_density.setter
def emission_density(self, value: Sequence[float]):
value = np.asarray(value, dtype=float)
if value.ndim != 1:
raise ValueError("emission_density must be a 1-D array")
if np.any(value < 0):
raise ValueError("emission_density values cannot be negative")
self._emission_density = value
@property
def energy(self) -> list[Univariate]:
return self._energy
@energy.setter
def energy(self, value: Univariate | Sequence[Univariate]):
if isinstance(value, Univariate):
self._energy = [value]
else:
cv.check_iterable_type('energy distributions', value, Univariate)
self._energy = list(value)
@property
def time(self) -> Univariate | None:
return self._time
@time.setter
def time(self, value: Univariate | None):
if value is not None:
cv.check_type('time distribution', value, Univariate)
self._time = value
@property
def phi_start(self) -> float:
return self._phi_start
@phi_start.setter
def phi_start(self, value: float):
cv.check_type('phi_start', value, Real)
self._phi_start = value
@property
def phi_extent(self) -> float:
return self._phi_extent
@phi_extent.setter
def phi_extent(self, value: float):
cv.check_type('phi_extent', value, Real)
cv.check_greater_than('phi_extent', value, 0.0)
cv.check_less_than('phi_extent', value, 2.0 * np.pi, equality=True)
self._phi_extent = value
@property
def n_alpha(self) -> int:
return self._n_alpha
@n_alpha.setter
def n_alpha(self, value: int):
cv.check_type('n_alpha', value, Integral)
cv.check_greater_than('n_alpha', value, 2)
if value < 51:
warnings.warn(
"n_alpha values below 51 may introduce noticeable "
"discretization bias in tokamak source sampling", stacklevel=2)
self._n_alpha = value
@property
def vertical_shift(self) -> float:
return self._vertical_shift
@vertical_shift.setter
def vertical_shift(self, value: float):
cv.check_type('vertical shift', value, Real)
self._vertical_shift = value
def populate_xml_element(self, element):
"""Add necessary tokamak source information to an XML element
Returns
-------
element : lxml.etree._Element
XML element containing source data
"""
self._validate()
# Geometry parameters
ET.SubElement(element, "major_radius").text = str(self.major_radius)
ET.SubElement(element, "minor_radius").text = str(self.minor_radius)
ET.SubElement(element, "elongation").text = str(self.elongation)
ET.SubElement(element, "triangularity").text = str(self.triangularity)
ET.SubElement(element, "shafranov_shift").text = str(self.shafranov_shift)
# Toroidal angle bounds
ET.SubElement(element, "phi_start").text = str(self.phi_start)
ET.SubElement(element, "phi_extent").text = str(self.phi_extent)
# Poloidal sampling resolution
ET.SubElement(element, "n_alpha").text = str(self.n_alpha)
# Vertical shift
if self.vertical_shift != 0.0:
ET.SubElement(element, "vertical_shift").text = str(self.vertical_shift)
# Emission profile
ET.SubElement(element, "r_over_a").text = ' '.join(str(r) for r in self.r_over_a)
ET.SubElement(element, "emission_density").text = ' '.join(str(s) for s in self.emission_density)
# Energy distribution(s)
for dist in self.energy:
element.append(dist.to_xml_element('energy'))
# Time distribution
if self.time is not None:
element.append(self.time.to_xml_element('time'))
@classmethod
def from_xml_element(cls, elem: ET.Element) -> TokamakSource:
"""Generate tokamak source from an XML element
Parameters
----------
elem : lxml.etree._Element
XML element
Returns
-------
openmc.TokamakSource
Source generated from XML element
"""
# Read geometry parameters
major_radius = float(get_text(elem, 'major_radius'))
minor_radius = float(get_text(elem, 'minor_radius'))
elongation = float(get_text(elem, 'elongation'))
triangularity = float(get_text(elem, 'triangularity'))
shafranov_shift = float(get_text(elem, 'shafranov_shift'))
# Read optional parameters
phi_start_text = get_text(elem, 'phi_start')
phi_start = float(phi_start_text) if phi_start_text else 0.0
phi_extent_text = get_text(elem, 'phi_extent')
phi_extent = float(phi_extent_text) if phi_extent_text else 2.0 * np.pi
n_alpha_text = get_text(elem, 'n_alpha')
n_alpha = int(n_alpha_text) if n_alpha_text else 101
vertical_shift_text = get_text(elem, 'vertical_shift')
vertical_shift = float(vertical_shift_text) if vertical_shift_text else 0.0
# Read emission profile
r_over_a = np.array([float(x) for x in get_text(elem, 'r_over_a').split()])
emission_density = np.array([float(x) for x in get_text(elem, 'emission_density').split()])
# Read energy distributions
energy = [Univariate.from_xml_element(e) for e in elem.findall('energy')]
if len(energy) == 1:
energy = energy[0]
# Read time distribution
time_elem = elem.find('time')
time = Univariate.from_xml_element(time_elem) if time_elem is not None else None
# Read constraints and strength
constraints = cls._get_constraints(elem)
strength_text = get_text(elem, 'strength')
strength = float(strength_text) if strength_text else 1.0
return cls(
major_radius=major_radius,
minor_radius=minor_radius,
elongation=elongation,
triangularity=triangularity,
shafranov_shift=shafranov_shift,
r_over_a=r_over_a,
emission_density=emission_density,
energy=energy,
time=time,
phi_start=phi_start,
phi_extent=phi_extent,
n_alpha=n_alpha,
vertical_shift=vertical_shift,
strength=strength,
constraints=constraints
)
class SourceParticle:

309
openmc/stats/univariate.py
View file

@ -2,8 +2,10 @@ from __future__ import annotations
from abc import ABC, abstractmethod
from collections import defaultdict
from collections.abc import Iterable, Sequence
from functools import cache
from math import sqrt, pi, exp, log
from numbers import Real
from pathlib import Path
from warnings import warn
import lxml.etree as ET
@ -14,6 +16,7 @@ import scipy
import openmc.checkvalue as cv
from openmc.data import atomic_mass, NEUTRON_MASS
import openmc.data
from .._xml import get_elem_list, get_text
from ..mixin import EqualityMixin
@ -124,6 +127,8 @@ class Univariate(EqualityMixin, ABC):
return Legendre.from_xml_element(elem)
elif distribution == 'mixture':
return Mixture.from_xml_element(elem)
elif distribution == 'decay_spectrum':
return DecaySpectrum.from_xml_element(elem)
@abstractmethod
def _sample_unbiased(self, n_samples: int = 1, seed: int | None = None):
@ -2196,6 +2201,310 @@ class Mixture(Univariate):
return new_dist
class DecaySpectrum(Univariate):
"""Energy distribution from decay photon spectra of a mixture of nuclides.
This distribution stores nuclide names, their atom densities, and the volume
of the region. When written to XML and read by the C++ solver, the nuclide
names are resolved against the depletion chain to obtain the decay photon
energy spectra and decay constants. The resulting distribution is a mixture
of per-nuclide photon spectra weighted by absolute activity. The volume is
necessary so that the C++ solver can compute the total photon emission rate
in [photons/s], which is used as the source strength.
.. versionadded:: 0.15.4
Parameters
----------
nuclides : dict
Dictionary mapping nuclide name (str) to atom density (float) in units
of [atom/b-cm].
volume : float
Volume of the source region in [cm³]. Used together with atom densities
to compute the absolute photon emission rate.
Attributes
----------
nuclides : dict
Dictionary mapping nuclide name to atom density in [atom/b-cm].
volume : float
Volume of the source region in [cm³].
"""
def __init__(self, nuclides: dict[str, float], volume: float):
super().__init__(bias=None)
self._dist_cache = None
self._dist_cache_key = None
self.nuclides = nuclides
self.volume = volume
def __len__(self):
return len(self.nuclides)
@property
def nuclides(self):
return self._nuclides
@nuclides.setter
def nuclides(self, nuclides):
cv.check_type('nuclides', nuclides, dict)
for name, density in nuclides.items():
cv.check_type('nuclide name', name, str)
cv.check_type(f'atom density for {name}', density, Real)
cv.check_greater_than(f'atom density for {name}', density, 0.0, True)
self._nuclides = dict(nuclides)
self._dist_cache = None
self._dist_cache_key = None
@property
def volume(self):
return self._volume
@volume.setter
def volume(self, volume):
cv.check_type('volume', volume, Real)
cv.check_greater_than('volume', volume, 0.0)
self._volume = float(volume)
self._dist_cache = None
self._dist_cache_key = None
@staticmethod
def _chain_file_cache_key():
"""Return a hashable key for the active depletion chain."""
chain_file = openmc.config.get('chain_file')
if chain_file is None:
return None
path = Path(chain_file).resolve()
try:
stat = path.stat()
except OSError:
return (path, None, None)
return (path, stat.st_mtime, stat.st_size)
def to_distribution(self):
"""Convert to a concrete distribution using decay chain data.
Builds a combined photon energy distribution by looking up each nuclide
in the depletion chain via :func:`openmc.data.decay_photon_energy` and
weighting by absolute atom count (``density * 1e24 * volume``). The
result is cached on the object; the cache is invalidated automatically
when :attr:`nuclides` or :attr:`volume` are reassigned.
Requires ``openmc.config['chain_file']`` to be set.
Returns
-------
openmc.stats.Univariate or None
Combined photon energy distribution, or ``None`` if no nuclide in
:attr:`nuclides` has a photon source in the chain.
"""
chain_key = self._chain_file_cache_key()
if self._dist_cache is not None and self._dist_cache_key == chain_key:
return self._dist_cache
dists = []
weights = []
for name, density in self.nuclides.items():
dist = openmc.data.decay_photon_energy(name)
if dist is not None:
dists.append(dist)
weights.append(density * 1e24 * self.volume)
if not dists:
return None
self._dist_cache = combine_distributions(dists, weights)
self._dist_cache_key = chain_key
return self._dist_cache
def to_xml_element(self, element_name: str):
"""Return XML representation of the decay photon distribution
Parameters
----------
element_name : str
XML element name
Returns
-------
element : lxml.etree._Element
XML element containing decay photon distribution data
"""
element = ET.Element(element_name)
element.set("type", "decay_spectrum")
element.set("volume", str(self.volume))
nuclides = ET.SubElement(element, "nuclides")
nuclides.text = ' '.join(self.nuclides)
parameters = ET.SubElement(element, "parameters")
parameters.text = ' '.join(str(density) for density in self.nuclides.values())
return element
@classmethod
def from_xml_element(cls, elem: ET.Element):
"""Generate decay photon distribution from an XML element
Parameters
----------
elem : lxml.etree._Element
XML element
Returns
-------
openmc.stats.DecaySpectrum
Decay photon distribution generated from XML element
"""
volume = float(elem.get('volume'))
names = get_elem_list(elem, 'nuclides', str)
densities = get_elem_list(elem, 'parameters', float)
nuclides = dict(zip(names, densities))
return cls(nuclides, volume)
def _sample_unbiased(self, n_samples=1, seed=None):
dist = self.to_distribution()
if dist is None:
raise RuntimeError(
"DecaySpectrum._sample_unbiased requires chain data but none "
"was found. Ensure openmc.config['chain_file'] is set and the "
"chain contains photon sources for the nuclides present."
)
return dist.sample(n_samples, seed)[0]
def integral(self):
"""Return integral of the distribution
Returns the total photon emission rate in [photons/s] by delegating to
:meth:`to_distribution`. Returns ``0.0`` when no chain data is
available (e.g., ``openmc.config['chain_file']`` is not set).
Returns
-------
float
Total photon emission rate in [photons/s], or ``0.0`` if chain
data is unavailable.
"""
try:
dist = self.to_distribution()
except Exception:
return 0.0
if dist is None:
return 0.0
return dist.integral()
@staticmethod
@cache
def _photon_integral(nuclide: str, chain_key) -> float | None:
"""Return the per-atom photon emission integral for a nuclide"""
dist = openmc.data.decay_photon_energy(nuclide)
return dist.integral() if dist is not None else None
def clip(self, tolerance: float = 1e-9, inplace: bool = False):
"""Remove nuclides with negligible contribution to photon emission.
Nuclides that are stable or have no photon source in the depletion
chain are removed unconditionally. The remaining nuclides are ranked
by their photon emission rate (proportional to
``atom_density * decay_constant * photon_yield``) and the least
important are discarded until the cumulative discarded fraction of the
total emission rate exceeds *tolerance*.
Requires ``openmc.config['chain_file']`` to be set.
Parameters
----------
tolerance : float
Maximum fraction of total photon emission rate that may be
discarded.
inplace : bool
Whether to modify the current object in-place or return a new one.
Returns
-------
openmc.stats.DecaySpectrum
Distribution with negligible nuclides removed.
"""
# Compute per-nuclide emission rate; drop non-emitters
emitting_names = []
emitting_densities = []
rates = []
chain_key = self._chain_file_cache_key()
for name, density in self.nuclides.items():
integral = DecaySpectrum._photon_integral(name, chain_key)
if integral is None:
continue
emitting_names.append(name)
emitting_densities.append(density)
rates.append(density * self.volume * integral)
if not emitting_names:
new_nuclides = {}
else:
indices = _intensity_clip(rates, tolerance=tolerance)
new_nuclides = {
emitting_names[i]: emitting_densities[i] for i in indices
}
if inplace:
self._nuclides = new_nuclides
self._dist_cache = None
self._dist_cache_key = None
return self
return type(self)(new_nuclides, self.volume)
@property
def support(self):
return (0.0, np.inf)
def evaluate(self, x):
"""Evaluate the probability density at a given value.
Delegates to the combined distribution built from chain data. Raises
``NotImplementedError`` if the combined distribution is a
:class:`~openmc.stats.Mixture` (which does not support
``evaluate()``).
Parameters
----------
x : float
Value at which to evaluate the PDF.
Returns
-------
float
Probability density at *x*.
"""
dist = self.to_distribution()
if dist is None:
raise RuntimeError(
"DecaySpectrum.evaluate requires chain data. Ensure "
"openmc.config['chain_file'] is set."
)
return dist.evaluate(x)
def mean(self):
"""Return the mean of the distribution.
Delegates to the combined distribution built from chain data.
Returns
-------
float
Mean photon energy in [eV].
"""
dist = self.to_distribution()
if dist is None:
raise RuntimeError(
"DecaySpectrum.mean requires chain data. Ensure "
"openmc.config['chain_file'] is set."
)
return dist.mean()
def combine_distributions(
dists: Sequence[Discrete | Tabular | Mixture],
probs: Sequence[float]

3
pyproject.toml
View file

@ -40,6 +40,7 @@ dependencies = [
[project.optional-dependencies]
depletion-mpi = ["mpi4py"]
docs = [
"breathe",
"sphinx",
"sphinxcontrib-katex",
"sphinx-numfig",
@ -69,7 +70,7 @@ include = ['openmc*']
exclude = ['tests*']
[tool.setuptools.package-data]
"openmc.data.dose" = ["**/*.txt"]
"openmc.data.dose" = ["**/*.txt", "*.h5"]
"openmc.data" = ["*.txt", "*.DAT", "*.json", "*.h5"]
"openmc.lib" = ["libopenmc.dylib", "libopenmc.so"]

202
src/bank.cpp
View file

@ -7,6 +7,7 @@
#include "openmc/vector.h"
#include <cstdint>
#include <limits>
#include <numeric>
namespace openmc {
@ -43,6 +44,13 @@ vector<vector<double>> ifp_fission_lifetime_bank;
// used to efficiently sort the fission bank after each iteration.
vector<int64_t> progeny_per_particle;
// When shared secondary bank mode is enabled, secondaries produced during
// transport are collected in the write bank. When a secondary generation is
// complete, write is moved to read for transport, and a new empty write bank
// is created. This repeats until no secondaries remain.
SharedArray<SourceSite> shared_secondary_bank_read;
SharedArray<SourceSite> shared_secondary_bank_write;
} // namespace simulation
//==============================================================================
@ -60,6 +68,8 @@ void free_memory_bank()
simulation::ifp_source_lifetime_bank.clear();
simulation::ifp_fission_delayed_group_bank.clear();
simulation::ifp_fission_lifetime_bank.clear();
simulation::shared_secondary_bank_read.clear();
simulation::shared_secondary_bank_write.clear();
}
void init_fission_bank(int64_t max)
@ -68,13 +78,13 @@ void init_fission_bank(int64_t max)
simulation::progeny_per_particle.resize(simulation::work_per_rank);
}
// Performs an O(n) sort on the fission bank, by leveraging
// Performs an O(n) sort on a fission or secondary bank, by leveraging
// the parent_id and progeny_id fields of banked particles. See the following
// paper for more details:
// "Reproducibility and Monte Carlo Eigenvalue Calculations," F.B. Brown and
// T.M. Sutton, 1992 ANS Annual Meeting, Transactions of the American Nuclear
// Society, Volume 65, Page 235.
void sort_fission_bank()
void sort_bank(SharedArray<SourceSite>& bank, bool is_fission_bank)
{
// Ensure we don't read off the end of the array if we ran with 0 particles
if (simulation::progeny_per_particle.size() == 0) {
@ -96,44 +106,200 @@ void sort_fission_bank()
vector<vector<double>> sorted_ifp_lifetime_bank;
// If there is not enough space, allocate a temporary vector and point to it
if (simulation::fission_bank.size() >
simulation::fission_bank.capacity() / 2) {
sorted_bank_holder.resize(simulation::fission_bank.size());
if (bank.size() > bank.capacity() / 2) {
sorted_bank_holder.resize(bank.size());
sorted_bank = sorted_bank_holder.data();
} else { // otherwise, point sorted_bank to unused portion of the fission bank
sorted_bank = &simulation::fission_bank[simulation::fission_bank.size()];
sorted_bank = bank.data() + bank.size();
}
if (settings::ifp_on) {
if (settings::ifp_on && is_fission_bank) {
allocate_temporary_vector_ifp(
sorted_ifp_delayed_group_bank, sorted_ifp_lifetime_bank);
}
// Use parent and progeny indices to sort fission bank
for (int64_t i = 0; i < simulation::fission_bank.size(); i++) {
const auto& site = simulation::fission_bank[i];
int64_t offset = site.parent_id - 1 - simulation::work_index[mpi::rank];
int64_t idx = simulation::progeny_per_particle[offset] + site.progeny_id;
if (idx >= simulation::fission_bank.size()) {
// Use parent and progeny indices to sort bank
for (int64_t i = 0; i < bank.size(); i++) {
const auto& site = bank[i];
if (site.parent_id < 0 ||
site.parent_id >=
static_cast<int64_t>(simulation::progeny_per_particle.size())) {
fatal_error(fmt::format("Invalid parent_id {} for banked site (expected "
"range [0, {})).",
site.parent_id, simulation::progeny_per_particle.size()));
}
int64_t idx =
simulation::progeny_per_particle[site.parent_id] + site.progeny_id;
if (idx < 0 || idx >= bank.size()) {
fatal_error("Mismatch detected between sum of all particle progeny and "
"shared fission bank size.");
"bank size during sorting.");
}
sorted_bank[idx] = site;
if (settings::ifp_on) {
if (settings::ifp_on && is_fission_bank) {
copy_ifp_data_from_fission_banks(
i, sorted_ifp_delayed_group_bank[idx], sorted_ifp_lifetime_bank[idx]);
}
}
// Copy sorted bank into the fission bank
std::copy(sorted_bank, sorted_bank + simulation::fission_bank.size(),
simulation::fission_bank.data());
if (settings::ifp_on) {
std::copy(sorted_bank, sorted_bank + bank.size(), bank.data());
if (settings::ifp_on && is_fission_bank) {
copy_ifp_data_to_fission_banks(
sorted_ifp_delayed_group_bank.data(), sorted_ifp_lifetime_bank.data());
}
}
// This function redistributes SourceSite particles across MPI ranks to
// achieve load balancing while preserving the global ordering of particles.
//
// GUARANTEES:
// -----------
// 1. Global Order Preservation: After redistribution, each rank holds a
// contiguous slice of the original global ordering. For example, if the
// input across 3 ranks was:
// - Rank 0: IDs 0-4
// - Rank 1: IDs 5-6
// - Rank 2: IDs 7-200
// Then after redistribution (assuming ~67 particles per rank):
// - Rank 0: IDs 0-66 (contiguous)
// - Rank 1: IDs 67-133 (contiguous)
// - Rank 2: IDs 134-200 (contiguous)
// The global ordering is always preserved - no rank will ever hold
// non-contiguous ID ranges like "0-4 and 100-200".
//
// 2. Even Load Balancing: Particles are distributed as evenly as possible.
// If total % n_procs != 0, the first 'remainder' ranks each get one extra
// particle (i.e., floor division with remainder distributed to lower
// ranks). This follows the same logic as calculate_work().
//
// HOW IT WORKS:
// -------------
// The algorithm uses overlap-based redistribution:
// 1. Each rank's current data occupies a range [cumulative_before[rank],
// cumulative_before[rank+1]) in the global index space.
// 2. Each rank's target data should occupy [cumulative_target[rank],
// cumulative_target[rank+1]) in the same global index space.
// 3. For each pair of (source_rank, dest_rank), we calculate the overlap
// between what source_rank currently has and what dest_rank needs.
// 4. MPI_Alltoallv transfers exactly these overlapping regions, with
// displacements ensuring data lands at the correct position in the
// receiving buffer.
//
// EDGE CASES HANDLED:
// -------------------
// - Single rank (n_procs == 1): Returns immediately with local size, no MPI.
// - Empty total (all ranks have 0 particles): Returns 0 immediately.
// - Imbalanced input (e.g., one rank has all particles): Works correctly;
// that rank will send portions to all other ranks based on target ranges.
// - Non-divisible totals: First 'remainder' ranks get one extra particle.
int64_t synchronize_global_secondary_bank(
SharedArray<SourceSite>& shared_secondary_bank)
{
// Get current size of local bank
int64_t local_size = shared_secondary_bank.size();
if (mpi::n_procs == 1) {
return local_size;
}
#ifdef OPENMC_MPI
// Gather all sizes to all ranks
vector<int64_t> all_sizes(mpi::n_procs);
MPI_Allgather(&local_size, 1, MPI_INT64_T, all_sizes.data(), 1, MPI_INT64_T,
mpi::intracomm);
// Calculate total and check for empty case
int64_t total = 0;
for (int64_t size : all_sizes) {
total += size;
}
// If we don't have any items to distribute, return
if (total == 0) {
return total;
}
int64_t base_count = total / mpi::n_procs;
int64_t remainder = total % mpi::n_procs;
// Calculate target size for each rank
// First 'remainder' ranks get base_count + 1, rest get base_count
vector<int64_t> target_sizes(mpi::n_procs);
for (int i = 0; i < mpi::n_procs; ++i) {
target_sizes[i] = base_count + (i < remainder ? 1 : 0);
}
// Calculate send and receive counts in terms of SourceSite objects
// (not bytes)
vector<int> send_counts(mpi::n_procs, 0);
vector<int> recv_counts(mpi::n_procs, 0);
vector<int> send_displs(mpi::n_procs, 0);
vector<int> recv_displs(mpi::n_procs, 0);
// Calculate cumulative positions (starting index for each rank in the
// global array)
vector<int64_t> cumulative_before(mpi::n_procs + 1, 0);
vector<int64_t> cumulative_target(mpi::n_procs + 1, 0);
for (int i = 0; i < mpi::n_procs; ++i) {
cumulative_before[i + 1] = cumulative_before[i] + all_sizes[i];
cumulative_target[i + 1] = cumulative_target[i] + target_sizes[i];
}
// Determine send and receive amounts for each rank
int64_t my_start = cumulative_before[mpi::rank];
int64_t my_end = cumulative_before[mpi::rank + 1];
int64_t my_target_start = cumulative_target[mpi::rank];
int64_t my_target_end = cumulative_target[mpi::rank + 1];
for (int r = 0; r < mpi::n_procs; ++r) {
// Send: overlap between my current range and rank r's target range
int64_t send_overlap_start = std::max(my_start, cumulative_target[r]);
int64_t send_overlap_end = std::min(my_end, cumulative_target[r + 1]);
if (send_overlap_start < send_overlap_end) {
int64_t count = send_overlap_end - send_overlap_start;
int64_t displ = send_overlap_start - my_start;
if (count > std::numeric_limits<int>::max() ||
displ > std::numeric_limits<int>::max()) {
fatal_error("Secondary bank size exceeds MPI_Alltoallv int limit.");
}
send_counts[r] = static_cast<int>(count);
send_displs[r] = static_cast<int>(displ);
}
// Recv: overlap between rank r's current range and my target range
int64_t recv_overlap_start =
std::max(cumulative_before[r], my_target_start);
int64_t recv_overlap_end =
std::min(cumulative_before[r + 1], my_target_end);
if (recv_overlap_start < recv_overlap_end) {
int64_t count = recv_overlap_end - recv_overlap_start;
int64_t displ = recv_overlap_start - my_target_start;
if (count > std::numeric_limits<int>::max() ||
displ > std::numeric_limits<int>::max()) {
fatal_error("Secondary bank size exceeds MPI_Alltoallv int limit.");
}
recv_counts[r] = static_cast<int>(count);
recv_displs[r] = static_cast<int>(displ);
}
}
// Prepare receive buffer with target size
SharedArray<SourceSite> new_bank(target_sizes[mpi::rank]);
// Perform all-to-all redistribution using the custom MPI type
MPI_Alltoallv(shared_secondary_bank.data(), send_counts.data(),
send_displs.data(), mpi::source_site, new_bank.data(), recv_counts.data(),
recv_displs.data(), mpi::source_site, mpi::intracomm);
// Replace old bank with redistributed data
shared_secondary_bank = std::move(new_bank);
return total;
#else
return local_size;
#endif
}
//==============================================================================
// C API
//==============================================================================

10
src/boundary_condition.cpp
View file

@ -38,6 +38,9 @@ void VacuumBC::handle_particle(Particle& p, const Surface& surf) const
void ReflectiveBC::handle_particle(Particle& p, const Surface& surf) const
{
Direction u = surf.reflect(p.r(), p.u(), &p);
// normalize reflected u to ensure no floating point error leads to
// unnormalized directions
u /= u.norm();
// Handle the effects of the surface albedo on the particle's weight.
@ -53,6 +56,9 @@ void ReflectiveBC::handle_particle(Particle& p, const Surface& surf) const
void WhiteBC::handle_particle(Particle& p, const Surface& surf) const
{
Direction u = surf.diffuse_reflect(p.r(), p.u(), p.current_seed());
// normalize outgoing u to ensure no floating point error leads to
// unnormalized directions
u /= u.norm();
// Handle the effects of the surface albedo on the particle's weight.
@ -241,6 +247,10 @@ void RotationalPeriodicBC::handle_particle(
new_u[axis_1_idx_] = cos_theta * u[axis_1_idx_] - sin_theta * u[axis_2_idx_];
new_u[axis_2_idx_] = sin_theta * u[axis_1_idx_] + cos_theta * u[axis_2_idx_];
// normalize new_u to ensure no floating point error leads to unnormalized
// directions
new_u /= new_u.norm();
// Handle the effects of the surface albedo on the particle's weight.
BoundaryCondition::handle_albedo(p, surf);

94
src/cell.cpp
View file

@ -340,6 +340,63 @@ void Cell::to_hdf5(hid_t cell_group) const
close_group(group);
}
//==============================================================================
// XML parsing helpers for <cell> nodes
//==============================================================================
vector<int32_t> parse_cell_material_xml(pugi::xml_node node, int32_t cell_id)
{
vector<std::string> mats {
get_node_array<std::string>(node, "material", true)};
if (mats.empty()) {
fatal_error(fmt::format(
"An empty material element was specified for cell {}", cell_id));
}
vector<int32_t> material;
material.reserve(mats.size());
for (const auto& mat : mats) {
if (mat == "void") {
material.push_back(MATERIAL_VOID);
} else {
material.push_back(std::stoi(mat));
}
}
return material;
}
vector<double> parse_cell_temperature_xml(pugi::xml_node node, int32_t cell_id)
{
auto temperatures = get_node_array<double>(node, "temperature");
if (temperatures.empty()) {
fatal_error(fmt::format(
"An empty temperature element was specified for cell {}", cell_id));
}
for (auto T : temperatures) {
if (T < 0) {
fatal_error(fmt::format(
"Cell {} was specified with a negative temperature", cell_id));
}
}
return temperatures;
}
vector<double> parse_cell_density_xml(pugi::xml_node node, int32_t cell_id)
{
auto densities = get_node_array<double>(node, "density");
if (densities.empty()) {
fatal_error(fmt::format(
"An empty density element was specified for cell {}", cell_id));
}
for (auto rho : densities) {
if (rho <= 0) {
fatal_error(fmt::format(
"Cell {} was specified with a density less than or equal to zero",
cell_id));
}
}
return densities;
}
//==============================================================================
// CSGCell implementation
//==============================================================================
@ -390,26 +447,12 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
// universe), more than one material (distribmats), and some materials may
// be "void".
if (material_present) {
vector<std::string> mats {
get_node_array<std::string>(cell_node, "material", true)};
if (mats.size() > 0) {
material_.reserve(mats.size());
for (std::string mat : mats) {
if (mat.compare("void") == 0) {
material_.push_back(MATERIAL_VOID);
} else {
material_.push_back(std::stoi(mat));
}
}
} else {
fatal_error(fmt::format(
"An empty material element was specified for cell {}", id_));
}
material_ = parse_cell_material_xml(cell_node, id_);
}
// Read the temperature element which may be distributed like materials.
if (check_for_node(cell_node, "temperature")) {
sqrtkT_ = get_node_array<double>(cell_node, "temperature");
sqrtkT_ = parse_cell_temperature_xml(cell_node, id_);
sqrtkT_.shrink_to_fit();
// Make sure this is a material-filled cell.
@ -420,14 +463,6 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
id_));
}
// Make sure all temperatures are non-negative.
for (auto T : sqrtkT_) {
if (T < 0) {
fatal_error(fmt::format(
"Cell {} was specified with a negative temperature", id_));
}
}
// Convert to sqrt(k*T).
for (auto& T : sqrtkT_) {
T = std::sqrt(K_BOLTZMANN * T);
@ -440,7 +475,7 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
// Note: calculating the actual density multiplier is deferred until materials
// are finalized. density_mult_ contains the true density in the meantime.
if (check_for_node(cell_node, "density")) {
density_mult_ = get_node_array<double>(cell_node, "density");
density_mult_ = parse_cell_density_xml(cell_node, id_);
density_mult_.shrink_to_fit();
// Make sure this is a material-filled cell.
@ -461,15 +496,6 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
id_));
}
}
// Make sure all densities are non-negative and greater than zero.
for (auto rho : density_mult_) {
if (rho <= 0) {
fatal_error(fmt::format(
"Cell {} was specified with a density less than or equal to zero",
id_));
}
}
}
// Read the region specification.

8
src/chain.cpp
View file

@ -98,6 +98,8 @@ vector<unique_ptr<ChainNuclide>> chain_nuclides;
void read_chain_file_xml()
{
free_memory_chain();
char* chain_file_path = std::getenv("OPENMC_CHAIN_FILE");
if (!chain_file_path) {
return;
@ -120,4 +122,10 @@ void read_chain_file_xml()
}
}
void free_memory_chain()
{
data::chain_nuclides.clear();
data::chain_nuclide_map.clear();
}
} // namespace openmc

12
src/collision_track.cpp
View file

@ -200,8 +200,13 @@ void collision_track_record(Particle& particle)
return;
int cell_id = model::cells[cell_index]->id_;
const auto* nuclide_ptr = data::nuclides[particle.event_nuclide()].get();
std::string nuclide = nuclide_ptr->name_;
std::string nuclide {};
int nuclide_id = 0;
if (particle.event_nuclide() != NUCLIDE_NONE) {
const auto* nuclide_ptr = data::nuclides[particle.event_nuclide()].get();
nuclide = nuclide_ptr->name_;
nuclide_id = nuclide_ptr->particle_type().pdg_number();
}
int universe_id = model::universes[particle.lowest_coord().universe()]->id_;
double delta_E = particle.E_last() - particle.E();
int material_index = particle.material();
@ -224,8 +229,7 @@ void collision_track_record(Particle& particle)
site.event_mt = particle.event_mt();
site.delayed_group = particle.delayed_group();
site.cell_id = cell_id;
site.nuclide_id =
10000 * nuclide_ptr->Z_ + 10 * nuclide_ptr->A_ + nuclide_ptr->metastable_;
site.nuclide_id = nuclide_id;
site.material_id = material_id;
site.universe_id = universe_id;
site.n_collision = particle.n_collision();

227
src/dagmc.cpp
View file

@ -50,6 +50,10 @@ namespace openmc {
DAGUniverse::DAGUniverse(pugi::xml_node node)
{
MaterialOverrides material_overrides;
TemperatureOverrides temperature_overrides;
DensityOverrides density_overrides;
if (check_for_node(node, "id")) {
id_ = std::stoi(get_node_value(node, "id"));
} else {
@ -76,24 +80,77 @@ DAGUniverse::DAGUniverse(pugi::xml_node node)
adjust_material_ids_ = get_node_value_bool(node, "auto_mat_ids");
}
// get material assignment overloading
if (check_for_node(node, "material_overrides")) {
auto mat_node = node.child("material_overrides");
// loop over all subelements (each subelement corresponds to a material)
for (pugi::xml_node cell_node : mat_node.children("cell_override")) {
// Store assignment reference name
int32_t ref_assignment = std::stoi(get_node_value(cell_node, "id"));
// Get material assignment overrides from nested DAGMC cell elements.
if (node.child("cell")) {
for (pugi::xml_node cell_node : node.children("cell")) {
if (!check_for_node(cell_node, "id")) {
fatal_error(
"Must specify id for each DAGMC cell override in <dagmc_universe>.");
}
// Get mat name for each assignement instances
vector<int32_t> instance_mats =
get_node_array<int32_t>(cell_node, "material_ids");
int32_t cell_id = std::stoi(get_node_value(cell_node, "id"));
// Store mat name for each instances
material_overrides_.emplace(ref_assignment, instance_mats);
if (check_for_node(cell_node, "region")) {
fatal_error(fmt::format(
"DAGMC cell {} override cannot specify a region.", cell_id));
}
if (check_for_node(cell_node, "fill")) {
fatal_error(fmt::format(
"DAGMC cell {} override currently only supports material fills.",
cell_id));
}
if (check_for_node(cell_node, "universe")) {
fatal_error(fmt::format(
"DAGMC cell {} override cannot specify a universe.", cell_id));
}
if (check_for_node(cell_node, "translation") ||
check_for_node(cell_node, "rotation")) {
fatal_error(fmt::format(
"DAGMC cell {} override does not support translation or rotation.",
cell_id));
}
if (!check_for_node(cell_node, "material")) {
fatal_error(fmt::format(
"DAGMC cell {} override must specify material.", cell_id));
}
auto inserted = material_overrides.emplace(
cell_id, parse_cell_material_xml(cell_node, cell_id));
if (!inserted.second) {
fatal_error(fmt::format(
"Duplicate DAGMC cell override specified for cell {}", cell_id));
}
if (check_for_node(cell_node, "temperature")) {
temperature_overrides.emplace(
cell_id, parse_cell_temperature_xml(cell_node, cell_id));
}
if (check_for_node(cell_node, "density")) {
density_overrides.emplace(
cell_id, parse_cell_density_xml(cell_node, cell_id));
}
}
} else if (check_for_node(node, "material_overrides")) {
if (node.child("cell")) {
fatal_error("DAGMCUniverse cannot specify both <material_overrides> and "
"<cell> sub-elements. Use <cell> elements only.");
}
warning("DAGMCUniverse <material_overrides> is deprecated. Use nested "
"<cell> elements under <dagmc_universe> instead.");
for (pugi::xml_node co :
node.child("material_overrides").children("cell_override")) {
int32_t cell_id = std::stoi(get_node_value(co, "id"));
std::istringstream iss(co.child("material_ids").text().get());
vector<int32_t> mats;
for (std::string s; iss >> s;) {
mats.push_back(s == "void" ? MATERIAL_VOID : std::stoi(s));
}
material_overrides.emplace(cell_id, mats);
}
}
initialize();
initialize(material_overrides, temperature_overrides, density_overrides);
}
DAGUniverse::DAGUniverse(
@ -110,9 +167,12 @@ DAGUniverse::DAGUniverse(std::shared_ptr<moab::DagMC> dagmc_ptr,
: dagmc_instance_(dagmc_ptr), filename_(filename),
adjust_geometry_ids_(auto_geom_ids), adjust_material_ids_(auto_mat_ids)
{
MaterialOverrides material_overrides;
TemperatureOverrides temperature_overrides;
DensityOverrides density_overrides;
set_id();
init_metadata();
init_geometry();
init_geometry(material_overrides, temperature_overrides, density_overrides);
}
void DAGUniverse::set_id()
@ -130,6 +190,15 @@ void DAGUniverse::set_id()
}
void DAGUniverse::initialize()
{
MaterialOverrides material_overrides;
TemperatureOverrides temperature_overrides;
initialize(material_overrides, temperature_overrides);
}
void DAGUniverse::initialize(const MaterialOverrides& material_overrides,
const TemperatureOverrides& temperature_overrides,
const DensityOverrides& density_overrides)
{
#ifdef OPENMC_UWUW_ENABLED
// read uwuw materials from the .h5m file if present
@ -140,7 +209,7 @@ void DAGUniverse::initialize()
init_metadata();
init_geometry();
init_geometry(material_overrides, temperature_overrides, density_overrides);
}
void DAGUniverse::init_dagmc()
@ -176,7 +245,9 @@ void DAGUniverse::init_metadata()
MB_CHK_ERR_CONT(rval);
}
void DAGUniverse::init_geometry()
void DAGUniverse::init_geometry(const MaterialOverrides& material_overrides,
const TemperatureOverrides& temperature_overrides,
const DensityOverrides& density_overrides)
{
moab::ErrorCode rval;
@ -202,6 +273,9 @@ void DAGUniverse::init_geometry()
: dagmc_instance_->id_by_index(3, c->dag_index());
c->universe_ = this->id_;
c->fill_ = C_NONE; // no fill, single universe
if (dagmc_instance_->is_implicit_complement(vol_handle)) {
c->name_ = "implicit complement";
}
auto in_map = model::cell_map.find(c->id_);
if (in_map == model::cell_map.end()) {
@ -230,16 +304,68 @@ void DAGUniverse::init_geometry()
if (mat_str == "graveyard") {
graveyard = vol_handle;
}
// material void checks
if (mat_str == "void" || mat_str == "vacuum" || mat_str == "graveyard") {
if (material_overrides.count(c->id_)) {
override_assign_material(c, material_overrides);
} else if (mat_str == "void" || mat_str == "vacuum" ||
mat_str == "graveyard") {
c->material_.push_back(MATERIAL_VOID);
} else if (uses_uwuw()) {
uwuw_assign_material(vol_handle, c);
} else {
if (material_overrides_.count(c->id_)) {
override_assign_material(c);
} else if (uses_uwuw()) {
uwuw_assign_material(vol_handle, c);
} else {
legacy_assign_material(mat_str, c);
legacy_assign_material(mat_str, c);
}
if (temperature_overrides.count(c->id_)) {
if (c->material_.empty() || c->material_[0] == MATERIAL_VOID) {
fatal_error(fmt::format("DAGMC cell {} was specified with a "
"temperature but no non-void material.",
c->id_));
}
c->sqrtkT_.clear();
const auto& temp_overrides = temperature_overrides.at(c->id_);
c->sqrtkT_.reserve(temp_overrides.size());
for (auto T : temp_overrides) {
c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * T));
}
if (settings::verbosity >= 10) {
std::stringstream override_values;
for (size_t i = 0; i < temp_overrides.size(); ++i) {
if (i > 0) {
override_values << " ";
}
override_values << temp_overrides[i];
}
auto msg = fmt::format("Overriding DAGMC cell {} property "
"'temperature [K]' with value(s): {}",
c->id_, override_values.str());
write_message(msg, 10);
}
}
if (density_overrides.count(c->id_)) {
if (c->material_.empty() || c->material_[0] == MATERIAL_VOID) {
fatal_error(fmt::format("DAGMC cell {} was specified with a density "
"but no non-void material.",
c->id_));
}
// density_mult_ holds the true density until materials are finalized,
// at which point it is converted to a proper multiplier (same as CSG).
c->density_mult_ = density_overrides.at(c->id_);
if (settings::verbosity >= 10) {
const auto& dens = density_overrides.at(c->id_);
std::stringstream override_values;
for (size_t i = 0; i < dens.size(); ++i) {
if (i > 0)
override_values << " ";
override_values << dens[i];
}
write_message(fmt::format("Overriding DAGMC cell {} property "
"'density [g/cm³]' with value(s): {}",
c->id_, override_values.str()),
10);
}
}
@ -252,18 +378,21 @@ void DAGUniverse::init_geometry()
continue;
}
// assign cell temperature
const auto& mat = model::materials[model::material_map.at(c->material_[0])];
if (dagmc_instance_->has_prop(vol_handle, "temp")) {
rval = dagmc_instance_->prop_value(vol_handle, "temp", temp_value);
MB_CHK_ERR_CONT(rval);
double temp = std::stod(temp_value);
c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * temp));
} else if (mat->temperature() > 0.0) {
c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * mat->temperature()));
} else {
c->sqrtkT_.push_back(
std::sqrt(K_BOLTZMANN * settings::temperature_default));
// assign cell temperature if not explicitly overridden
if (c->sqrtkT_.empty()) {
const auto& mat =
model::materials[model::material_map.at(c->material_[0])];
if (dagmc_instance_->has_prop(vol_handle, "temp")) {
rval = dagmc_instance_->prop_value(vol_handle, "temp", temp_value);
MB_CHK_ERR_CONT(rval);
double temp = std::stod(temp_value);
c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * temp));
} else if (mat->temperature() > 0.0) {
c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * mat->temperature()));
} else {
c->sqrtkT_.push_back(
std::sqrt(K_BOLTZMANN * settings::temperature_default));
}
}
model::cells.emplace_back(std::move(c));
@ -630,7 +759,8 @@ void DAGUniverse::uwuw_assign_material(
#endif // OPENMC_UWUW_ENABLED
}
void DAGUniverse::override_assign_material(std::unique_ptr<DAGCell>& c) const
void DAGUniverse::override_assign_material(std::unique_ptr<DAGCell>& c,
const MaterialOverrides& material_overrides) const
{
// if Cell ID matches an override key, use it to override the material
// assignment else if UWUW is used, get the material assignment from the DAGMC
@ -638,17 +768,30 @@ void DAGUniverse::override_assign_material(std::unique_ptr<DAGCell>& c) const
// Notify User that an override is being applied on a DAGMCCell
write_message(fmt::format("Applying override for DAGMCCell {}", c->id_), 8);
const auto& mat_overrides = material_overrides.at(c->id_);
if (settings::verbosity >= 10) {
auto msg = fmt::format("Assigning DAGMC cell {} material(s) based on "
"override information (see input XML).",
c->id_);
std::stringstream override_values;
for (size_t i = 0; i < mat_overrides.size(); ++i) {
if (i > 0) {
override_values << " ";
}
if (mat_overrides[i] == MATERIAL_VOID) {
override_values << "void";
} else {
override_values << mat_overrides[i];
}
}
auto msg = fmt::format("Overriding DAGMC cell {} property 'material' "
"with value(s): {}",
c->id_, override_values.str());
write_message(msg, 10);
}
// Override the material assignment for each cell instance using the legacy
// assignement
for (auto mat_id : material_overrides_.at(c->id_)) {
if (model::material_map.find(mat_id) == model::material_map.end()) {
for (auto mat_id : mat_overrides) {
if (mat_id != MATERIAL_VOID &&
model::material_map.find(mat_id) == model::material_map.end()) {
fatal_error(fmt::format(
"Material with ID '{}' not found for DAGMC cell {}", mat_id, c->id_));
}

135
src/distribution.cpp
View file

@ -7,7 +7,9 @@
#include <numeric> // for accumulate
#include <stdexcept> // for runtime_error
#include <string> // for string, stod
#include <unordered_set>
#include "openmc/chain.h"
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/math_functions.h"
@ -15,6 +17,10 @@
#include "openmc/random_lcg.h"
#include "openmc/xml_interface.h"
namespace {
std::unordered_set<std::string> decay_spectrum_missing_chain_nuclides;
}
namespace openmc {
//==============================================================================
@ -546,14 +552,9 @@ double Tabular::sample_unbiased(uint64_t* seed) const
double c = prn(seed);
// Find first CDF bin which is above the sampled value
double c_i = c_[0];
int i;
std::size_t n = c_.size();
for (i = 0; i < n - 1; ++i) {
if (c <= c_[i + 1])
break;
c_i = c_[i + 1];
}
auto c_iter = std::lower_bound(c_.begin() + 1, c_.end(), c);
int i = std::distance(c_.begin(), c_iter) - 1;
double c_i = c_[i];
// Determine bounding PDF values
double x_i = x_[i];
@ -758,6 +759,8 @@ UPtrDist distribution_from_xml(pugi::xml_node node)
dist = UPtrDist {new Tabular(node)};
} else if (type == "mixture") {
dist = UPtrDist {new Mixture(node)};
} else if (type == "decay_spectrum") {
dist = UPtrDist {new DecaySpectrum(node)};
} else if (type == "muir") {
openmc::fatal_error(
"'muir' distributions are now specified using the openmc.stats.muir() "
@ -768,4 +771,120 @@ UPtrDist distribution_from_xml(pugi::xml_node node)
return dist;
}
//==============================================================================
// DecaySpectrum implementation
//==============================================================================
DecaySpectrum::DecaySpectrum(pugi::xml_node node)
{
// Read the region volume [cm^3] needed for absolute emission rate
if (!check_for_node(node, "volume"))
fatal_error("DecaySpectrum: 'volume' attribute is required.");
double volume = std::stod(get_node_value(node, "volume"));
// Read nuclide names and atom densities from XML
vector<int> nuclide_indices;
vector<double> atoms;
auto names = get_node_array<std::string>(node, "nuclides");
auto densities = get_node_array<double>(node, "parameters");
if (names.size() != densities.size()) {
fatal_error("DecaySpectrum nuclides and parameters must have the same "
"length.");
}
for (size_t i = 0; i < names.size(); ++i) {
const auto& name = names[i];
double density = densities[i];
// Look up nuclide in the depletion chain
auto it = data::chain_nuclide_map.find(name);
if (it == data::chain_nuclide_map.end()) {
if (decay_spectrum_missing_chain_nuclides.insert(name).second) {
warning("Nuclide '" + name +
"' appears in a DecaySpectrum source but is not present in "
"the depletion chain; it will be ignored.");
}
continue;
}
int nuclide_index = it->second;
const auto& chain_nuc = data::chain_nuclides[nuclide_index];
const Distribution* photon_dist = chain_nuc->photon_energy();
if (!photon_dist)
continue;
// Skip non-positive densities and warn if negative
if (density <= 0.0) {
if (density < 0.0) {
warning("Nuclide '" + name +
"' has a negative density in a DecaySpectrum source; it will "
"be ignored.");
}
continue;
}
// atoms = density [atom/b-cm] * 1e24 [b/cm^2] * volume [cm^3]
double atoms_i = density * 1.0e24 * volume;
nuclide_indices.push_back(nuclide_index);
atoms.push_back(atoms_i);
}
init(std::move(nuclide_indices), atoms);
}
void DecaySpectrum::init(
vector<int> nuclide_indices, const vector<double>& atoms)
{
if (nuclide_indices.size() != atoms.size()) {
fatal_error("DecaySpectrum nuclide index and atoms arrays must have "
"the same length.");
}
vector<double> probs;
probs.reserve(nuclide_indices.size());
for (size_t i = 0; i < nuclide_indices.size(); ++i) {
// Distribution integral is in [photons/s/atom]; multiplying by atoms gives
// the total emission rate [photons/s] for this nuclide.
const auto* dist =
data::chain_nuclides[nuclide_indices[i]]->photon_energy();
probs.push_back(atoms[i] * dist->integral());
}
nuclide_indices_ = std::move(nuclide_indices);
integral_ = std::accumulate(probs.begin(), probs.end(), 0.0);
if (nuclide_indices_.empty() || integral_ <= 0.0) {
fatal_error("DecaySpectrum source did not resolve any nuclides with decay "
"photon spectra and positive atom densities. Ensure "
"OPENMC_CHAIN_FILE is set and matches the nuclides in the "
"source definition.");
}
di_.assign(probs);
}
DecaySpectrum::Sample DecaySpectrum::sample_with_parent(uint64_t* seed) const
{
size_t idx = di_.sample(seed);
int parent_nuclide = nuclide_indices_[idx];
const auto* dist = data::chain_nuclides[parent_nuclide]->photon_energy();
auto [energy, weight] = dist->sample(seed);
return {energy, weight, parent_nuclide};
}
std::pair<double, double> DecaySpectrum::sample(uint64_t* seed) const
{
auto sample = sample_with_parent(seed);
return {sample.energy, sample.weight};
}
double DecaySpectrum::integral() const
{
return integral_;
}
double DecaySpectrum::sample_unbiased(uint64_t* seed) const
{
return sample_with_parent(seed).energy;
}
} // namespace openmc

8
src/distribution_spatial.cpp
View file

@ -262,9 +262,11 @@ std::pair<Position, double> SphericalIndependent::sample(uint64_t* seed) const
MeshSpatial::MeshSpatial(pugi::xml_node node)
{
if (get_node_value(node, "type", true, true) != "mesh") {
fatal_error(fmt::format(
"Incorrect spatial type '{}' for a MeshSpatial distribution"));
auto spatial_type = get_node_value(node, "type", true, true);
if (spatial_type != "mesh") {
fatal_error(
fmt::format("Incorrect spatial type '{}' for a MeshSpatial distribution",
spatial_type));
}
// No in-tet distributions implemented, could include distributions for the

51
src/event.cpp
View file

@ -1,6 +1,9 @@
#include "openmc/event.h"
#include "openmc/bank.h"
#include "openmc/error.h"
#include "openmc/material.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/timer.h"
@ -64,7 +67,8 @@ void process_init_events(int64_t n_particles, int64_t source_offset)
simulation::time_event_init.start();
#pragma omp parallel for schedule(runtime)
for (int64_t i = 0; i < n_particles; i++) {
initialize_history(simulation::particles[i], source_offset + i + 1);
initialize_particle_track(
simulation::particles[i], source_offset + i + 1, false);
dispatch_xs_event(i);
}
simulation::time_event_init.stop();
@ -136,7 +140,7 @@ void process_surface_crossing_events()
int64_t buffer_idx = simulation::surface_crossing_queue[i].idx;
Particle& p = simulation::particles[buffer_idx];
p.event_cross_surface();
p.event_revive_from_secondary();
p.event_check_limit_and_revive();
if (p.alive())
dispatch_xs_event(buffer_idx);
}
@ -155,7 +159,7 @@ void process_collision_events()
int64_t buffer_idx = simulation::collision_queue[i].idx;
Particle& p = simulation::particles[buffer_idx];
p.event_collide();
p.event_revive_from_secondary();
p.event_check_limit_and_revive();
if (p.alive())
dispatch_xs_event(buffer_idx);
}
@ -176,4 +180,45 @@ void process_death_events(int64_t n_particles)
simulation::time_event_death.stop();
}
void process_transport_events()
{
while (true) {
int64_t max = std::max({simulation::calculate_fuel_xs_queue.size(),
simulation::calculate_nonfuel_xs_queue.size(),
simulation::advance_particle_queue.size(),
simulation::surface_crossing_queue.size(),
simulation::collision_queue.size()});
if (max == 0) {
break;
} else if (max == simulation::calculate_fuel_xs_queue.size()) {
process_calculate_xs_events(simulation::calculate_fuel_xs_queue);
} else if (max == simulation::calculate_nonfuel_xs_queue.size()) {
process_calculate_xs_events(simulation::calculate_nonfuel_xs_queue);
} else if (max == simulation::advance_particle_queue.size()) {
process_advance_particle_events();
} else if (max == simulation::surface_crossing_queue.size()) {
process_surface_crossing_events();
} else if (max == simulation::collision_queue.size()) {
process_collision_events();
}
}
}
void process_init_secondary_events(int64_t n_particles, int64_t offset,
const SharedArray<SourceSite>& shared_secondary_bank)
{
simulation::time_event_init.start();
#pragma omp parallel for schedule(runtime)
for (int64_t i = 0; i < n_particles; i++) {
initialize_particle_track(simulation::particles[i], offset + i + 1, true);
const SourceSite& site = shared_secondary_bank[offset + i];
simulation::particles[i].event_revive_from_secondary(site);
if (simulation::particles[i].alive()) {
dispatch_xs_event(i);
}
}
simulation::time_event_init.stop();
}
} // namespace openmc

5
src/finalize.cpp
View file

@ -2,6 +2,7 @@
#include "openmc/bank.h"
#include "openmc/capi.h"
#include "openmc/chain.h"
#include "openmc/cmfd_solver.h"
#include "openmc/collision_track.h"
#include "openmc/constants.h"
@ -44,6 +45,7 @@ void free_memory()
free_memory_photon();
free_memory_settings();
free_memory_thermal();
free_memory_chain();
library_clear();
nuclides_clear();
free_memory_source();
@ -147,6 +149,8 @@ int openmc_finalize()
settings::uniform_source_sampling = false;
settings::ufs_on = false;
settings::urr_ptables_on = true;
settings::use_decay_photons = false;
settings::use_shared_secondary_bank = false;
settings::verbosity = -1;
settings::weight_cutoff = 0.25;
settings::weight_survive = 1.0;
@ -217,6 +221,7 @@ int openmc_reset()
settings::cmfd_run = false;
simulation::n_lost_particles = 0;
simulation::simulation_tracks_completed = 0;
return 0;
}

37
src/geometry.cpp
View file

@ -26,16 +26,22 @@ int n_coord_levels;
vector<int64_t> overlap_check_count;
vector<OverlapKey> overlap_keys;
std::unordered_map<OverlapKey, int, OverlapKeyHash> overlap_key_index;
} // namespace model
//==============================================================================
// Non-member functions
//==============================================================================
bool check_cell_overlap(GeometryState& p, bool error)
int check_cell_overlap(GeometryState& p, bool error)
{
int n_coord = p.n_coord();
// If no overlap found, return a nonphysical index
int overlap_index = -1;
// Loop through each coordinate level
for (int j = 0; j < n_coord; j++) {
Universe& univ = *model::universes[p.coord(j).universe()];
@ -44,21 +50,40 @@ bool check_cell_overlap(GeometryState& p, bool error)
for (auto index_cell : univ.cells_) {
Cell& c = *model::cells[index_cell];
if (c.contains(p.coord(j).r(), p.coord(j).u(), p.surface())) {
#pragma omp atomic
++model::overlap_check_count[index_cell];
if (index_cell != p.coord(j).cell()) {
if (error) {
fatal_error(
fmt::format("Overlapping cells detected: {}, {} on universe {}",
c.id_, model::cells[p.coord(j).cell()]->id_, univ.id_));
}
return true;
// With no fatal error (plotter is calling), now adds overlaps and
// ensures order does not matter when making overlap key
int cell_a = model::cells[index_cell]->id_;
int cell_b = model::cells[p.coord(j).cell()]->id_;
int a = std::min(cell_a, cell_b);
int b = std::max(cell_a, cell_b);
OverlapKey key {univ.id_, a, b};
#pragma omp critical(overlap_key_update)
{
auto it = model::overlap_key_index.find(key);
if (it != model::overlap_key_index.end()) {
overlap_index = it->second; // already exists, reuse index
} else {
int idx = int(model::overlap_keys.size());
model::overlap_keys.push_back(key);
model::overlap_key_index[key] = idx;
overlap_index = idx;
}
}
break;
}
#pragma omp atomic
++model::overlap_check_count[index_cell];
}
}
}
return false;
return overlap_index;
}
//==============================================================================

4
src/ifp.cpp
View file

@ -32,13 +32,13 @@ void ifp(const Particle& p, int64_t idx)
{
if (is_beta_effective_or_both()) {
const auto& delayed_groups =
simulation::ifp_source_delayed_group_bank[p.current_work() - 1];
simulation::ifp_source_delayed_group_bank[p.current_work()];
simulation::ifp_fission_delayed_group_bank[idx] =
_ifp(p.delayed_group(), delayed_groups);
}
if (is_generation_time_or_both()) {
const auto& lifetimes =
simulation::ifp_source_lifetime_bank[p.current_work() - 1];
simulation::ifp_source_lifetime_bank[p.current_work()];
simulation::ifp_fission_lifetime_bank[idx] = _ifp(p.lifetime(), lifetimes);
}
}

78
src/initialize.cpp
View file

@ -161,7 +161,7 @@ void initialize_mpi(MPI_Comm intracomm)
// Create bank datatype
SourceSite b;
MPI_Aint disp[11];
MPI_Aint disp[14];
MPI_Get_address(&b.r, &disp[0]);
MPI_Get_address(&b.u, &disp[1]);
MPI_Get_address(&b.E, &disp[2]);
@ -173,14 +173,35 @@ void initialize_mpi(MPI_Comm intracomm)
MPI_Get_address(&b.parent_nuclide, &disp[8]);
MPI_Get_address(&b.parent_id, &disp[9]);
MPI_Get_address(&b.progeny_id, &disp[10]);
for (int i = 10; i >= 0; --i) {
MPI_Get_address(&b.wgt_born, &disp[11]);
MPI_Get_address(&b.wgt_ww_born, &disp[12]);
MPI_Get_address(&b.n_split, &disp[13]);
for (int i = 13; i >= 0; --i) {
disp[i] -= disp[0];
}
int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1};
MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE,
MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_LONG, MPI_LONG};
MPI_Type_create_struct(11, blocks, disp, types, &mpi::source_site);
// Block counts for each field
int blocks[] = {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
// Types for each field
MPI_Datatype types[] = {
MPI_DOUBLE, // r (3 doubles)
MPI_DOUBLE, // u (3 doubles)
MPI_DOUBLE, // E
MPI_DOUBLE, // time
MPI_DOUBLE, // wgt
MPI_INT, // delayed_group
MPI_INT, // surf_id
MPI_INT, // particle (enum)
MPI_INT, // parent_nuclide
MPI_INT64_T, // parent_id
MPI_INT64_T, // progeny_id
MPI_DOUBLE, // wgt_born
MPI_DOUBLE, // wgt_ww_born
MPI_INT64_T // n_split
};
MPI_Type_create_struct(14, blocks, disp, types, &mpi::source_site);
MPI_Type_commit(&mpi::source_site);
CollisionTrackSite bc;
@ -300,6 +321,11 @@ int parse_command_line(int argc, char* argv[])
settings::run_mode = RunMode::VOLUME;
} else if (arg == "-s" || arg == "--threads") {
// Read number of threads
if (i + 1 >= argc) {
std::string msg {"Number of threads not specified."};
strcpy(openmc_err_msg, msg.c_str());
return OPENMC_E_INVALID_ARGUMENT;
}
i += 1;
#ifdef _OPENMP
@ -361,6 +387,28 @@ int parse_command_line(int argc, char* argv[])
return 0;
}
// TODO: Pulse-height tallies require per-history scoring across the full
// particle tree (parent + all descendants). The shared secondary bank
// transports each secondary as an independent Particle, breaking this
// assumption. A proper fix would defer pulse-height scoring: save
// (root_source_id, cell, pht_storage) per particle, then aggregate by
// root_source_id after all secondary generations complete before scoring
// into the histogram. For now, disable shared secondary when pulse-height
// tallies are present.
static void check_pulse_height_compatibility()
{
if (settings::use_shared_secondary_bank) {
for (const auto& t : model::tallies) {
if (t->type_ == TallyType::PULSE_HEIGHT) {
settings::use_shared_secondary_bank = false;
warning("Pulse-height tallies are not yet compatible with the shared "
"secondary bank. Disabling shared secondary bank.");
break;
}
}
}
}
bool read_model_xml()
{
std::string model_filename = settings::path_input;
@ -405,6 +453,10 @@ bool read_model_xml()
write_message(
fmt::format("Reading model XML file '{}' ...", model_filename), 5);
// Read chain data before settings so DecaySpectrum source distributions can
// resolve nuclides while sources are constructed.
read_chain_file_xml();
read_settings_xml(settings_root);
// If other XML files are present, display warning
@ -420,9 +472,6 @@ bool read_model_xml()
}
}
// Read data from chain file
read_chain_file_xml();
// Read materials and cross sections
if (!check_for_node(root, "materials")) {
fatal_error(fmt::format(
@ -454,6 +503,8 @@ bool read_model_xml()
if (check_for_node(root, "tallies"))
read_tallies_xml(root.child("tallies"));
check_pulse_height_compatibility();
// Initialize distribcell_filters
prepare_distribcell();
@ -475,14 +526,15 @@ bool read_model_xml()
void read_separate_xml_files()
{
// Read chain data before settings so DecaySpectrum source distributions can
// resolve nuclides while sources are constructed.
read_chain_file_xml();
read_settings_xml();
if (settings::run_mode != RunMode::PLOTTING) {
read_cross_sections_xml();
}
// Read data from chain file
read_chain_file_xml();
read_materials_xml();
read_geometry_xml();
@ -498,6 +550,8 @@ void read_separate_xml_files()
read_tallies_xml();
check_pulse_height_compatibility();
// Initialize distribcell_filters
prepare_distribcell();

33
src/lattice.cpp
View file

@ -10,6 +10,7 @@
#include "openmc/geometry.h"
#include "openmc/geometry_aux.h"
#include "openmc/hdf5_interface.h"
#include "openmc/math_functions.h"
#include "openmc/string_utils.h"
#include "openmc/vector.h"
#include "openmc/xml_interface.h"
@ -260,25 +261,33 @@ std::pair<double, array<int, 3>> RectLattice::distance(
// Determine the oncoming edge.
double x0 {copysign(0.5 * pitch_[0], u.x)};
double y0 {copysign(0.5 * pitch_[1], u.y)};
double z0;
double d = std::min(
u.x != 0.0 ? (x0 - x) / u.x : INFTY, u.y != 0.0 ? (y0 - y) / u.y : INFTY);
// Evaluate the distance to each oncoming edge independently. Comparing these
// distances directly (rather than reconstructing the crossing position)
// avoids the floating-point cancellation that occurs for large pitches.
double dx = u.x != 0.0 ? (x0 - x) / u.x : INFTY;
double dy = u.y != 0.0 ? (y0 - y) / u.y : INFTY;
double dz = INFTY;
if (is_3d_) {
z0 = copysign(0.5 * pitch_[2], u.z);
d = std::min(d, u.z != 0.0 ? (z0 - z) / u.z : INFTY);
double z0 {copysign(0.5 * pitch_[2], u.z)};
dz = u.z != 0.0 ? (z0 - z) / u.z : INFTY;
}
// Determine which lattice boundaries are being crossed
// The distance to the nearest lattice boundary is the smallest axial
// distance.
double d = std::min({dx, dy, dz});
// Determine which lattice boundaries are being crossed. The axis attaining
// the minimum is exactly equal to d, so at least one translation is always
// set for a finite crossing; a near-equal second axis indicates a corner
// crossing.
array<int, 3> lattice_trans = {0, 0, 0};
if (u.x != 0.0 && std::abs(x + u.x * d - x0) < FP_PRECISION)
if (isclose(d, dx, FP_COINCIDENT, FP_PRECISION))
lattice_trans[0] = copysign(1, u.x);
if (u.y != 0.0 && std::abs(y + u.y * d - y0) < FP_PRECISION)
if (isclose(d, dy, FP_COINCIDENT, FP_PRECISION))
lattice_trans[1] = copysign(1, u.y);
if (is_3d_) {
if (u.z != 0.0 && std::abs(z + u.z * d - z0) < FP_PRECISION)
lattice_trans[2] = copysign(1, u.z);
}
if (is_3d_ && isclose(d, dz, FP_COINCIDENT, FP_PRECISION))
lattice_trans[2] = copysign(1, u.z);
return {d, lattice_trans};
}

50
src/math_functions.cpp
View file

@ -1,5 +1,7 @@
#include "openmc/math_functions.h"
#include <limits> // for numeric_limits
#include "openmc/external/Faddeeva.hh"
#include "openmc/constants.h"
@ -944,6 +946,46 @@ double log1prel(double x)
}
}
double cyl_bessel_j(int n, double x)
{
// Handle negative arguments via the parity relation
// J_n(-x) = (-1)^n J_n(x); std::cyl_bessel_j has a domain error for x < 0.
double sign = 1.0;
if (x < 0.0) {
x = -x;
if (n % 2 == 1)
sign = -1.0;
}
#if defined(__cpp_lib_math_special_functions) && \
__cpp_lib_math_special_functions >= 201603L
return sign * std::cyl_bessel_j(static_cast<double>(n), x);
#else
// Ascending power series (e.g., Abramowitz & Stegun eq. 9.1.10):
// J_n(x) = sum_{m=0}^inf (-1)^m / (m! (m+n)!) * (x/2)^(2m+n)
// The term ratio is -(x/2)^2 / (m*(m+n)), so for |x| <= 2 the series
// converges to machine precision within ~20 terms.
double half_x = 0.5 * x;
// First term: (x/2)^n / n!
double term = 1.0;
for (int k = 1; k <= n; ++k) {
term *= half_x / k;
}
double sum = term;
double neg_half_x_sq = -half_x * half_x;
for (int m = 1; m <= 50; ++m) {
term *= neg_half_x_sq / (m * (m + n));
sum += term;
if (std::abs(term) <=
std::numeric_limits<double>::epsilon() * std::abs(sum))
break;
}
return sign * sum;
#endif
}
// Helper function to get index and interpolation function on an incident energy
// grid
void get_energy_index(
@ -959,4 +1001,12 @@ void get_energy_index(
}
}
// Return true if two floating-point values are approximately equal within a
// combined relative and absolute tolerance.
bool isclose(double a, double b, double rel_tol, double abs_tol)
{
return std::abs(a - b) <=
std::max(rel_tol * std::max(std::abs(a), std::abs(b)), abs_tol);
}
} // namespace openmc

26
src/mesh.cpp
View file

@ -6,7 +6,8 @@
#define _USE_MATH_DEFINES // to make M_PI declared in Intel and MSVC compilers
#include <cmath> // for ceil
#include <cstddef> // for size_t
#include <numeric> // for accumulate
#include <limits>
#include <numeric> // for accumulate
#include <string>
#ifdef _MSC_VER
@ -1080,13 +1081,26 @@ int StructuredMesh::get_bin(Position r) const
int StructuredMesh::n_bins() const
{
return std::accumulate(
shape_.begin(), shape_.begin() + n_dimension_, 1, std::multiplies<>());
// Bin indices are stored as 32-bit ints in the tally system.
int64_t n = 1;
for (int i = 0; i < n_dimension_; ++i)
n *= shape_[i];
if (n > std::numeric_limits<int>::max()) {
fatal_error(fmt::format(
"Mesh {} has too many bins ({}) for 32-bit tally indexing", id_, n));
}
return static_cast<int>(n);
}
int StructuredMesh::n_surface_bins() const
{
return 4 * n_dimension_ * n_bins();
// Surface bin indices are stored as 32-bit ints in the tally system.
int64_t n = static_cast<int64_t>(n_bins()) * 4 * n_dimension_;
if (n > std::numeric_limits<int>::max()) {
fatal_error(fmt::format(
"Mesh {} has too many surface bins ({}) for tally indexing", id_, n));
}
return static_cast<int>(n);
}
tensor::Tensor<double> StructuredMesh::count_sites(
@ -3670,7 +3684,7 @@ Position LibMesh::sample_element(int32_t bin, uint64_t* seed) const
// Get tet vertex coordinates from LibMesh
std::array<Position, 4> tet_verts;
for (int i = 0; i < elem.n_nodes(); i++) {
auto node_ref = elem.node_ref(i);
const auto& node_ref = elem.node_ref(i);
tet_verts[i] = {node_ref(0), node_ref(1), node_ref(2)};
}
// Samples position within tet using Barycentric coordinates
@ -3700,7 +3714,7 @@ int LibMesh::n_vertices() const
Position LibMesh::vertex(int vertex_id) const
{
const auto node_ref = m_->node_ref(vertex_id);
const auto& node_ref = m_->node_ref(vertex_id);
if (length_multiplier_ > 0.0) {
return length_multiplier_ * Position(node_ref(0), node_ref(1), node_ref(2));
} else {

17
src/output.cpp
View file

@ -501,6 +501,14 @@ void print_runtime()
show_rate("Calculation Rate (inactive)", speed_inactive);
}
show_rate("Calculation Rate (active)", speed_active);
// Display track rate when weight windows are enabled
if (settings::weight_windows_on) {
double speed_tracks =
simulation::simulation_tracks_completed / time_active.elapsed();
fmt::print(
" {:<33} = {:.6} tracks/second\n", "Track Rate (active)", speed_tracks);
}
}
//==============================================================================
@ -613,8 +621,15 @@ void write_tallies()
if (model::tallies.empty())
return;
// Tag tallies.out written during the forward solve of an adjoint run
const char* forward =
(FlatSourceDomain::solve_ == RandomRaySolve::FORWARD_FOR_ADJOINT)
? "forward."
: "";
// Set filename for tallies_out
std::string filename = fmt::format("{}tallies.out", settings::path_output);
std::string filename =
fmt::format("{}tallies.{}out", settings::path_output, forward);
// Open the tallies.out file.
std::ofstream tallies_out;

206
src/particle.cpp
View file

@ -94,7 +94,7 @@ bool Particle::create_secondary(
// Increment number of secondaries created (for ParticleProductionFilter)
n_secondaries()++;
auto& bank = secondary_bank().emplace_back();
SourceSite bank;
bank.particle = type;
bank.wgt = wgt;
bank.r = r();
@ -102,12 +102,21 @@ bool Particle::create_secondary(
bank.E = settings::run_CE ? E : g();
bank.time = time();
bank_second_E() += bank.E;
bank.parent_id = current_work();
if (settings::use_shared_secondary_bank) {
bank.progeny_id = n_progeny()++;
}
bank.wgt_born = wgt_born();
bank.wgt_ww_born = wgt_ww_born();
bank.n_split = n_split();
local_secondary_bank().emplace_back(bank);
return true;
}
void Particle::split(double wgt)
{
auto& bank = secondary_bank().emplace_back();
SourceSite bank;
bank.particle = type();
bank.wgt = wgt;
bank.r = r();
@ -122,6 +131,16 @@ void Particle::split(double wgt)
int surf_id = model::surfaces[surface_index()]->id_;
bank.surf_id = (surface() > 0) ? surf_id : -surf_id;
}
bank.wgt_born = wgt_born();
bank.wgt_ww_born = wgt_ww_born();
bank.n_split = n_split();
bank.parent_id = current_work();
if (settings::use_shared_secondary_bank) {
bank.progeny_id = n_progeny()++;
}
local_secondary_bank().emplace_back(bank);
}
void Particle::from_source(const SourceSite* src)
@ -168,6 +187,10 @@ void Particle::from_source(const SourceSite* src)
int index_plus_one = model::surface_map[std::abs(src->surf_id)] + 1;
surface() = (src->surf_id > 0) ? index_plus_one : -index_plus_one;
}
wgt_born() = src->wgt_born;
wgt_ww_born() = src->wgt_ww_born;
n_split() = src->n_split;
}
void Particle::event_calculate_xs()
@ -315,13 +338,14 @@ void Particle::event_cross_surface()
boundary().lattice_translation()[2] != 0) {
// Particle crosses lattice boundary
int i_lattice = coord(boundary().coord_level() - 1).lattice();
bool verbose = settings::verbosity >= 10 || trace();
cross_lattice(*this, boundary(), verbose);
event() = TallyEvent::LATTICE;
// Score cell to cell partial currents
if (!model::active_surface_tallies.empty()) {
auto& lat {*model::lattices[lowest_coord().lattice()]};
auto& lat {*model::lattices[i_lattice]};
bool is_valid;
Direction normal =
lat.get_normal(boundary().lattice_translation(), is_valid);
@ -449,61 +473,72 @@ void Particle::event_collide()
#endif
}
void Particle::event_revive_from_secondary()
void Particle::event_revive_from_secondary(const SourceSite& site)
{
// Write final position for the previous track (skip if this is a freshly
// constructed particle with no prior track, e.g., Phase 2 of shared
// secondary transport)
if (write_track() && n_event() > 0) {
write_particle_track(*this);
}
from_source(&site);
n_event() = 0;
if (!settings::use_shared_secondary_bank) {
n_tracks()++;
}
bank_second_E() = 0.0;
// Subtract secondary particle energy from interim pulse-height results.
// In shared secondary mode, this subtraction was already done on the parent
// particle during create_secondary(), so skip it here.
if (!settings::use_shared_secondary_bank &&
!model::active_pulse_height_tallies.empty() && this->type().is_photon()) {
// Since the birth cell of the particle has not been set we
// have to determine it before the energy of the secondary particle can be
// removed from the pulse-height of this cell.
if (lowest_coord().cell() == C_NONE) {
bool verbose = settings::verbosity >= 10 || trace();
if (!exhaustive_find_cell(*this, verbose)) {
mark_as_lost("Could not find the cell containing particle " +
std::to_string(id()));
return;
}
// Set birth cell attribute
if (cell_born() == C_NONE)
cell_born() = lowest_coord().cell();
// Initialize last cells from current cell
for (int j = 0; j < n_coord(); ++j) {
cell_last(j) = coord(j).cell();
}
n_coord_last() = n_coord();
}
pht_secondary_particles();
}
// Enter new particle in particle track file
if (write_track())
add_particle_track(*this);
}
void Particle::event_check_limit_and_revive()
{
// If particle has too many events, display warning and kill it
++n_event();
n_event()++;
if (n_event() == settings::max_particle_events) {
warning("Particle " + std::to_string(id()) +
" underwent maximum number of events.");
wgt() = 0.0;
}
// Check for secondary particles if this particle is dead
if (!alive()) {
// Write final position for this particle
if (write_track()) {
write_particle_track(*this);
}
// If no secondary particles, break out of event loop
if (secondary_bank().empty())
return;
from_source(&secondary_bank().back());
secondary_bank().pop_back();
n_event() = 0;
bank_second_E() = 0.0;
// Subtract secondary particle energy from interim pulse-height results
if (!model::active_pulse_height_tallies.empty() &&
this->type().is_photon()) {
// Since the birth cell of the particle has not been set we
// have to determine it before the energy of the secondary particle can be
// removed from the pulse-height of this cell.
if (lowest_coord().cell() == C_NONE) {
bool verbose = settings::verbosity >= 10 || trace();
if (!exhaustive_find_cell(*this, verbose)) {
mark_as_lost("Could not find the cell containing particle " +
std::to_string(id()));
return;
}
// Set birth cell attribute
if (cell_born() == C_NONE)
cell_born() = lowest_coord().cell();
// Initialize last cells from current cell
for (int j = 0; j < n_coord(); ++j) {
cell_last(j) = coord(j).cell();
}
n_coord_last() = n_coord();
}
pht_secondary_particles();
}
// Enter new particle in particle track file
if (write_track())
add_particle_track(*this);
// In non-shared-secondary mode, revive from local secondary bank
if (!alive() && !settings::use_shared_secondary_bank &&
!local_secondary_bank().empty()) {
SourceSite& site = local_secondary_bank().back();
event_revive_from_secondary(site);
local_secondary_bank().pop_back();
}
}
@ -515,18 +550,34 @@ void Particle::event_death()
// Finish particle track output.
if (write_track()) {
write_particle_track(*this);
finalize_particle_track(*this);
}
// Contribute tally reduction variables to global accumulator
// Contribute tally reduction variables to global accumulator
const auto k_absorption = keff_tally_absorption();
const auto k_collision = keff_tally_collision();
const auto k_tracklength = keff_tally_tracklength();
const auto leakage = keff_tally_leakage();
if (settings::run_mode == RunMode::EIGENVALUE) {
if (k_absorption != 0.0) {
#pragma omp atomic
global_tally_absorption += keff_tally_absorption();
global_tally_absorption += k_absorption;
}
if (k_collision != 0.0) {
#pragma omp atomic
global_tally_collision += keff_tally_collision();
global_tally_collision += k_collision;
}
if (k_tracklength != 0.0) {
#pragma omp atomic
global_tally_tracklength += keff_tally_tracklength();
global_tally_tracklength += k_tracklength;
}
}
if (leakage != 0.0) {
#pragma omp atomic
global_tally_leakage += keff_tally_leakage();
global_tally_leakage += leakage;
}
// Reset particle tallies once accumulated
keff_tally_absorption() = 0.0;
@ -538,11 +589,17 @@ void Particle::event_death()
score_pulse_height_tally(*this, model::active_pulse_height_tallies);
}
// Accumulate track count for this particle history
if (!settings::use_shared_secondary_bank) {
#pragma omp atomic
simulation::simulation_tracks_completed += n_tracks();
}
// Record the number of progeny created by this particle.
// This data will be used to efficiently sort the fission bank.
if (settings::run_mode == RunMode::EIGENVALUE) {
int64_t offset = id() - 1 - simulation::work_index[mpi::rank];
simulation::progeny_per_particle[offset] = n_progeny();
if (settings::run_mode == RunMode::EIGENVALUE ||
settings::use_shared_secondary_bank) {
simulation::progeny_per_particle[current_work()] = n_progeny();
}
}
@ -863,28 +920,27 @@ void Particle::write_restart() const
write_dataset(file_id, "id", id());
write_dataset(file_id, "type", type().pdg_number());
// Get source site data for the particle that got lost
int64_t i = current_work();
SourceSite site;
if (settings::run_mode == RunMode::EIGENVALUE) {
// take source data from primary bank for eigenvalue simulation
write_dataset(file_id, "weight", simulation::source_bank[i - 1].wgt);
write_dataset(file_id, "energy", simulation::source_bank[i - 1].E);
write_dataset(file_id, "xyz", simulation::source_bank[i - 1].r);
write_dataset(file_id, "uvw", simulation::source_bank[i - 1].u);
write_dataset(file_id, "time", simulation::source_bank[i - 1].time);
site = simulation::source_bank[i];
} else if (settings::run_mode == RunMode::FIXED_SOURCE &&
settings::use_shared_secondary_bank &&
i < simulation::shared_secondary_bank_read.size()) {
site = simulation::shared_secondary_bank_read[i];
} else if (settings::run_mode == RunMode::FIXED_SOURCE) {
// re-sample using rng random number seed used to generate source particle
int64_t id = (simulation::total_gen + overall_generation() - 1) *
settings::n_particles +
simulation::work_index[mpi::rank] + i;
// Re-sample using the same seed used to generate the source particle.
// current_work() is 0-indexed, compute_particle_id expects 1-indexed.
int64_t id = compute_transport_seed(compute_particle_id(i + 1));
uint64_t seed = init_seed(id, STREAM_SOURCE);
// re-sample source site
auto site = sample_external_source(&seed);
write_dataset(file_id, "weight", site.wgt);
write_dataset(file_id, "energy", site.E);
write_dataset(file_id, "xyz", site.r);
write_dataset(file_id, "uvw", site.u);
write_dataset(file_id, "time", site.time);
site = sample_external_source(&seed);
}
write_dataset(file_id, "weight", site.wgt);
write_dataset(file_id, "energy", site.E);
write_dataset(file_id, "xyz", site.r);
write_dataset(file_id, "uvw", site.u);
write_dataset(file_id, "time", site.time);
// Close file
file_close(file_id);

23
src/particle_restart.cpp
View file

@ -1,6 +1,7 @@
#include "openmc/particle_restart.h"
#include "openmc/array.h"
#include "openmc/bank.h"
#include "openmc/constants.h"
#include "openmc/hdf5_interface.h"
#include "openmc/mgxs_interface.h"
@ -106,20 +107,16 @@ void run_particle_restart()
// Set all tallies to 0 for now (just tracking errors)
model::tallies.clear();
// Compute random number seed
int64_t particle_seed;
switch (previous_run_mode) {
case RunMode::EIGENVALUE:
case RunMode::FIXED_SOURCE:
particle_seed = (simulation::total_gen + overall_generation() - 1) *
settings::n_particles +
p.id();
break;
default:
throw std::runtime_error {
"Unexpected run mode: " +
std::to_string(static_cast<int>(previous_run_mode))};
// Allocate progeny_per_particle if needed for shared secondary mode
// (event_death() writes to this array). Set current_work to 0 since we
// only have one particle being restarted.
if (settings::use_shared_secondary_bank) {
p.current_work() = 0;
simulation::progeny_per_particle.resize(1, 0);
}
// Compute random number seed
int64_t particle_seed = compute_transport_seed(p.id());
init_particle_seeds(particle_seed, p.seeds());
// Force calculation of cross-sections by setting last energy to zero

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