Merge branch 'develop' into batchwise_pr

.github/workflows/ci.yml

@@ -35,9 +35,6 @@ jobs:
         vectfit: [n]
 
         include:
-          - python-version: "3.7"
-            omp: n
-            mpi: n
           - python-version: "3.8"
             omp: n
             mpi: n
@@ -47,6 +44,9 @@ jobs:
           - python-version: "3.11"
             omp: n
             mpi: n
+          - python-version: "3.12"
+            omp: n
+            mpi: n
           - dagmc: y
             python-version: "3.10"
             mpi: y
@@ -111,16 +111,26 @@ jobs:
         run: |
           sudo apt -y update
           sudo apt install -y libpng-dev \
-                              libmpich-dev \
                               libnetcdf-dev \
                               libpnetcdf-dev \
                               libhdf5-serial-dev \
-                              libhdf5-mpich-dev \
                               libeigen3-dev
+
+      - name: Optional apt dependencies for MPI
+        shell: bash
+        if: ${{ matrix.mpi == 'y' }}
+        run: |
+          sudo apt install -y libhdf5-mpich-dev \
+                              libmpich-dev
           sudo update-alternatives --set mpi /usr/bin/mpicc.mpich
           sudo update-alternatives --set mpirun /usr/bin/mpirun.mpich
           sudo update-alternatives --set mpi-x86_64-linux-gnu /usr/include/x86_64-linux-gnu/mpich
 
+      - name: Optional apt dependencies for vectfit
+        shell: bash
+        if: ${{ matrix.vectfit == 'y' }}
+        run: sudo apt install -y libblas-dev liblapack-dev
+
       - name: install
         shell: bash
         run: |

CMakeLists.txt

@@ -39,6 +39,7 @@ option(OPENMC_USE_LIBMESH     "Enable support for libMesh unstructured mesh tall
 option(OPENMC_USE_MPI         "Enable MPI"                                           OFF)
 option(OPENMC_USE_MCPL        "Enable MCPL"                                          OFF)
 option(OPENMC_USE_NCRYSTAL    "Enable support for NCrystal scattering"               OFF)
+option(OPENMC_USE_UWUW        "Enable UWUW"                                          OFF)
 
 # Warnings for deprecated options
 foreach(OLD_OPT IN ITEMS "openmp" "profile" "coverage" "dagmc" "libmesh")
@@ -80,6 +81,14 @@ if(NOT CMAKE_BUILD_TYPE)
     set(CMAKE_BUILD_TYPE RelWithDebInfo CACHE STRING "Choose the type of build" FORCE)
 endif()
 
+#===============================================================================
+# OpenMP for shared-memory parallelism (and GPU support some day!)
+#===============================================================================
+
+if(OPENMC_USE_OPENMP)
+  find_package(OpenMP REQUIRED)
+endif()
+
 #===============================================================================
 # MPI for distributed-memory parallelism
 #===============================================================================
@@ -117,10 +126,15 @@ endif()
 if(OPENMC_USE_DAGMC)
   find_package(DAGMC REQUIRED PATH_SUFFIXES lib/cmake)
   if (${DAGMC_VERSION} VERSION_LESS 3.2.0)
-    message(FATAL_ERROR "Discovered DAGMC Version: ${DAGMC_VERSION}. \
-    Please update DAGMC to version 3.2.0 or greater.")
+    message(FATAL_ERROR "Discovered DAGMC Version: ${DAGMC_VERSION}."
+      "Please update DAGMC to version 3.2.0 or greater.")
   endif()
   message(STATUS "Found DAGMC: ${DAGMC_DIR} (version ${DAGMC_VERSION})")
+
+  # Check if UWUW is needed and available
+  if(OPENMC_USE_UWUW AND NOT DAGMC_BUILD_UWUW)
+    message(FATAL_ERROR "UWUW is enabled but DAGMC was not configured with UWUW.")
+  endif()
 endif()
 
 #===============================================================================
@@ -192,13 +206,6 @@ endif()
 # Skip for Visual Studio which has its own configurations through GUI
 if(NOT MSVC)
 
-if(OPENMC_USE_OPENMP)
-  find_package(OpenMP REQUIRED)
-  # In CMake 3.9+, can use the OpenMP::OpenMP_CXX imported target
-  list(APPEND cxxflags ${OpenMP_CXX_FLAGS})
-  list(APPEND ldflags ${OpenMP_CXX_FLAGS})
-endif()
-
 set(CMAKE_POSITION_INDEPENDENT_CODE ON)
 
 if(OPENMC_ENABLE_PROFILE)
@@ -380,6 +387,9 @@ list(APPEND libopenmc_SOURCES
   src/progress_bar.cpp
   src/random_dist.cpp
   src/random_lcg.cpp
+  src/random_ray/random_ray_simulation.cpp
+  src/random_ray/random_ray.cpp
+  src/random_ray/flat_source_domain.cpp
   src/reaction.cpp
   src/reaction_product.cpp
   src/scattdata.cpp
@@ -411,6 +421,7 @@ list(APPEND libopenmc_SOURCES
   src/tallies/filter_material.cpp
   src/tallies/filter_materialfrom.cpp
   src/tallies/filter_mesh.cpp
+  src/tallies/filter_meshborn.cpp
   src/tallies/filter_meshsurface.cpp
   src/tallies/filter_mu.cpp
   src/tallies/filter_particle.cpp
@@ -505,7 +516,7 @@ endif()
 
 if(OPENMC_USE_DAGMC)
   target_compile_definitions(libopenmc PRIVATE DAGMC)
-  target_link_libraries(libopenmc dagmc-shared uwuw-shared)
+  target_link_libraries(libopenmc dagmc-shared)
 endif()
 
 if(OPENMC_USE_LIBMESH)
@@ -518,6 +529,10 @@ if (PNG_FOUND)
   target_link_libraries(libopenmc PNG::PNG)
 endif()
 
+if (OPENMC_USE_OPENMP)
+  target_link_libraries(libopenmc OpenMP::OpenMP_CXX)
+endif()
+
 if (OPENMC_USE_MPI)
   target_link_libraries(libopenmc MPI::MPI_CXX)
 endif()
@@ -538,6 +553,11 @@ if(OPENMC_USE_NCRYSTAL)
   target_link_libraries(libopenmc NCrystal::NCrystal)
 endif()
 
+if (OPENMC_USE_UWUW)
+  target_compile_definitions(libopenmc PRIVATE UWUW)
+  target_link_libraries(libopenmc uwuw-shared)
+endif()
+
 #===============================================================================
 # Log build info that this executable can report later
 #===============================================================================

CODEOWNERS

@@ -5,9 +5,9 @@ openmc/data/ @paulromano
 openmc/lib/ @paulromano
 
 # Depletion
-openmc/deplete/ @drewejohnson
-tests/regression_tests/deplete/ @drewejohnson
-tests/unit_tests/test_deplete_*.py @drewejohnson
+openmc/deplete/ @paulromano
+tests/regression_tests/deplete/ @paulromano
+tests/unit_tests/test_deplete_*.py @paulromano
 
 # MG-related functionality
 openmc/mgxs_library.py @nelsonag
@@ -26,6 +26,12 @@ src/dagmc.cpp @pshriwise
 tests/regression_tests/dagmc/ @pshriwise
 tests/unit_tests/dagmc/ @pshriwise
 
+# Weight windows
+openmc/weight_windows.py @pshriwise
+openmc/lib/weight_windows.py @pshriwise
+src/weight_windows.py @pshriwise
+tests/unit_tests/weightwindows/ @pshriwise
+
 # Photon transport
 openmc/data/BREMX.DAT @amandalund
 openmc/data/compton_profiles.h5 @amandalund
@@ -49,3 +55,12 @@ openmc/data/resonance_covariance.py @icmeyer
 
 # Docker
 Dockerfile @shimwell
+
+# Random ray
+src/random_ray/ @jtramm
+
+# NCrystal interface
+src/ncrystal_interface.cpp @marquezj
+
+# MCPL interface
+src/mcpl_interface.cpp @ebknudsen

Dockerfile

@@ -24,7 +24,7 @@ ARG compile_cores=1
 ARG build_dagmc=off
 ARG build_libmesh=off
 
-FROM debian:bullseye-slim AS dependencies
+FROM debian:bookworm-slim AS dependencies
 
 ARG compile_cores
 ARG build_dagmc
@@ -34,21 +34,21 @@ ARG build_libmesh
 ENV HOME=/root
 
 # Embree variables
-ENV EMBREE_TAG='v3.12.2'
+ENV EMBREE_TAG='v4.3.1'
 ENV EMBREE_REPO='https://github.com/embree/embree'
 ENV EMBREE_INSTALL_DIR=$HOME/EMBREE/
 
 # MOAB variables
-ENV MOAB_TAG='5.3.0'
+ENV MOAB_TAG='5.5.1'
 ENV MOAB_REPO='https://bitbucket.org/fathomteam/moab/'
 
 # Double-Down variables
-ENV DD_TAG='v1.0.0'
+ENV DD_TAG='v1.1.0'
 ENV DD_REPO='https://github.com/pshriwise/double-down'
 ENV DD_INSTALL_DIR=$HOME/Double_down
 
 # DAGMC variables
-ENV DAGMC_BRANCH='v3.2.1'
+ENV DAGMC_BRANCH='v3.2.3'
 ENV DAGMC_REPO='https://github.com/svalinn/DAGMC'
 ENV DAGMC_INSTALL_DIR=$HOME/DAGMC/
 
@@ -71,9 +71,13 @@ RUN apt-get update -y && \
     apt-get install -y \
         python3-pip python-is-python3 wget git build-essential cmake \
         mpich libmpich-dev libhdf5-serial-dev libhdf5-mpich-dev \
-        libpng-dev && \
+        libpng-dev python3-venv && \
     apt-get autoremove
 
+# create virtual enviroment to avoid externally managed environment error
+RUN python3 -m venv openmc_venv
+ENV PATH=/openmc_venv/bin:$PATH
+
 # Update system-provided pip
 RUN pip install --upgrade pip
 

cmake/OpenMCConfig.cmake.in

@@ -31,6 +31,14 @@ if(@OPENMC_USE_MPI@)
   find_package(MPI REQUIRED)
 endif()
 
+if(@OPENMC_USE_OPENMP@)
+  find_package(OpenMP REQUIRED)
+endif()
+
 if(@OPENMC_USE_MCPL@)
   find_package(MCPL REQUIRED)
 endif()
+
+if(@OPENMC_USE_UWUW@)
+  find_package(UWUW REQUIRED)
+endif()

docs/source/capi/index.rst

@@ -420,6 +420,16 @@ Functions
    :return: Return status (negative if an error occurred)
    :rtype: int
 
+.. c:function:: int openmc_mesh_filter_get_mesh(int32_t index, int32_t* index_mesh)
+
+   Get the mesh for a mesh filter
+
+   :param int32_t index: Index in the filters array
+   :param index_mesh: Index in the meshes array
+   :type index_mesh: int32_t*
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
 .. c:function:: int openmc_mesh_filter_set_mesh(int32_t index, int32_t index_mesh)
 
    Set the mesh for a mesh filter
@@ -429,6 +439,98 @@ Functions
    :return: Return status (negative if an error occurred)
    :rtype: int
 
+.. c:function:: int openmc_mesh_filter_get_translation(int32_t index, double translation[3])
+
+   Get the 3-D translation coordinates for a mesh filter
+
+   :param int32_t index: Index in the filters array
+   :param double[3] translation: 3-D translation coordinates
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_mesh_filter_set_translation(int32_t index, double translation[3])
+
+   Set the 3-D translation coordinates for a mesh filter
+
+   :param int32_t index: Index in the filters array
+   :param double[3] translation: 3-D translation coordinates
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_meshborn_filter_get_mesh(int32_t index, int32_t* index_mesh)
+
+   Get the mesh for a meshborn filter
+
+   :param int32_t index: Index in the filters array
+   :param index_mesh: Index in the meshes array
+   :type index_mesh: int32_t*
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_meshborn_filter_set_mesh(int32_t index, int32_t index_mesh)
+
+   Set the mesh for a meshborn filter
+
+   :param int32_t index: Index in the filters array
+   :param int32_t index_mesh: Index in the meshes array
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_meshborn_filter_get_translation(int32_t index, double translation[3])
+
+   Get the 3-D translation coordinates for a meshborn filter
+
+   :param int32_t index: Index in the filters array
+   :param double[3] translation: 3-D translation coordinates
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_meshborn_filter_set_translation(int32_t index, double translation[3])
+
+   Set the 3-D translation coordinates for a meshborn filter
+
+   :param int32_t index: Index in the filters array
+   :param double[3] translation: 3-D translation coordinates
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh)
+
+   Get the mesh for a mesh surface filter
+
+   :param int32_t index: Index in the filters array
+   :param index_mesh: Index in the meshes array
+   :type index_mesh: int32_t*
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh)
+
+   Set the mesh for a mesh surface filter
+
+   :param int32_t index: Index in the filters array
+   :param int32_t index_mesh: Index in the meshes array
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_meshsurface_filter_get_translation(int32_t index, double translation[3])
+
+   Get the 3-D translation coordinates for a mesh surface filter
+
+   :param int32_t index: Index in the filters array
+   :param double[3] translation: 3-D translation coordinates
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_meshsurface_filter_set_translation(int32_t index, double translation[3])
+
+   Set the 3-D translation coordinates for a mesh surface filter
+
+   :param int32_t index: Index in the filters array
+   :param double[3] translation: 3-D translation coordinates
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
 .. c:function:: int openmc_next_batch()
 
    Simulate next batch of particles. Must be called after openmc_simulation_init().

docs/source/devguide/styleguide.rst

@@ -146,7 +146,7 @@ Style for Python code should follow PEP8_.
 
 Docstrings for functions and methods should follow numpydoc_ style.
 
-Python code should work with Python 3.7+.
+Python code should work with Python 3.8+.
 
 Use of third-party Python packages should be limited to numpy_, scipy_,
 matplotlib_, pandas_, and h5py_. Use of other third-party packages must be

docs/source/io_formats/settings.rst

@@ -252,9 +252,9 @@ to false.
 
   *Default*: true
 
-----------------------------------------
+-------------------------------------
 ``<max_particles_in_flight>`` Element
-----------------------------------------
+-------------------------------------
 
 This element indicates the number of neutrons to run in flight concurrently
 when using event-based parallelism. A higher value uses more memory, but
@@ -262,9 +262,17 @@ may be more efficient computationally.
 
   *Default*: 100000
 
----------------------------
+---------------------------------
+``<max_particle_events>`` Element
+---------------------------------
+
+This element indicates the maximum number of events a particle can undergo.
+
+  *Default*: 1000000
+
+-----------------------
 ``<max_order>`` Element
----------------------------
+-----------------------
 
 The ``<max_order>`` element allows the user to set a maximum scattering order
 to apply to every nuclide/material in the problem.  That is, if the data
@@ -276,11 +284,11 @@ then, OpenMC will only use up to the :math:`P_1` data.
   .. note:: This element is not used in the continuous-energy
     :ref:`energy_mode`.
 
----------------------------
+------------------------
 ``<max_splits>`` Element
----------------------------
+------------------------
 
-The ``<max_splits>`` element indicates the number of times a particle can split during a history. 
+The ``<max_splits>`` element indicates the number of times a particle can split during a history.
 
   *Default*: 1000
 
@@ -405,6 +413,32 @@ or sub-elements and can be set to either "false" or "true".
 
   .. note:: This element is not used in the multi-group :ref:`energy_mode`.
 
+------------------------
+``<random_ray>`` Element
+------------------------
+
+The ``<random_ray>`` element enables random ray mode and contains a number of
+settings relevant to the solver. Tips for selecting these parameters can be
+found in the :ref:`random ray user guide <random_ray>`.
+
+  :distance_inactive:
+    The inactive ray length (dead zone length) in [cm].
+
+    *Default*: None
+
+  :distance_active:
+    The active ray length in [cm].
+
+    *Default*: None
+
+  :source:
+    Specifies the starting ray distribution, and follows the format for
+    :ref:`source_element`. It must be uniform in space and angle and cover the
+    full domain. It does not represent a physical neutron or photon source -- it
+    is only used to sample integrating ray starting locations and directions.
+
+    *Default*: None
+
 ----------------------------------
 ``<resonance_scattering>`` Element
 ----------------------------------
@@ -500,8 +534,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``, or
+    ``mesh``. The type of the source will be determined by this attribute if it
+    is present.
 
   :particle:
     The source particle type, either ``neutron`` or ``photon``.
@@ -682,6 +717,39 @@ 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.
 
+  :constraints:
+    This sub-element indicates the presence of constraints on sampled source
+    sites (see :ref:`usersguide_source_constraints` for details). It may have
+    the following sub-elements:
+
+    :domain_ids:
+      The unique IDs of domains for which source sites must be within.
+
+      *Default*: None
+
+    :domain_type:
+      The type of each domain for source rejection ("cell", "material", or
+      "universe").
+
+      *Default*: None
+
+    :fissionable:
+      A boolean indicating whether source sites must be sampled within a
+      material that is fissionable in order to be accepted.
+
+    :time_bounds:
+      A pair of times in [s] indicating the lower and upper bound for a time
+      interval that source particles must be within.
+
+    :energy_bounds:
+      A pair of energies in [eV] indicating the lower and upper bound for an
+      energy interval that source particles must be within.
+
+    :rejection_strategy:
+      Either "resample", indicating that source sites should be resampled when
+      one is rejected, or "kill", indicating that a rejected source site is
+      assigned zero weight.
+
 .. _univariate:
 
 Univariate Probability Distributions

docs/source/io_formats/statepoint.rst

@@ -96,6 +96,8 @@ The current version of the statepoint file format is 18.1.
               - **library** (*char[]*) -- Mesh library used to represent the
                                           mesh ("moab" or "libmesh").
               - **length_multiplier** (*double*) Scaling factor applied to the mesh.
+              - **options** (*char[]*) -- Special options that control spatial
+                                          search data structures used.
               - **volumes** (*double[]*) -- Volume of each mesh cell.
               - **vertices** (*double[]*) -- x, y, z values of the mesh vertices.
               - **connectivity** (*int[]*) -- Connectivity array for the mesh

docs/source/io_formats/tallies.rst

@@ -100,6 +100,18 @@ The ``<tally>`` element accepts the following sub-elements:
 
      *Default*: None
 
+   :ignore_zeros:
+     Whether to allow zero tally bins to be ignored when assessing the
+     convergece of the precision trigger. If True, only nonzero tally scores
+     will be compared to the trigger's threshold.
+
+     .. note:: The ``ignore_zeros`` option can cause the tally trigger to fire
+               prematurely if there are no hits in any bins at the first
+               evalulation. It is the user's responsibility to specify enough
+               particles per batch to get a nonzero score in at least one bin.
+     
+     *Default*: False
+
    :scores:
      The score(s) in this tally to which the trigger should be applied.
 
@@ -364,6 +376,10 @@ attributes/sub-elements:
     The mesh library used to represent an unstructured mesh. This can be either
     "moab" or "libmesh". (For unstructured mesh only.)
 
+  :options:
+    Special options that control spatial search data structures used. (For
+    unstructured mesh using MOAB only)
+
   :filename:
     The name of the mesh file to be loaded at runtime. (For unstructured mesh
     only.)

docs/source/methods/index.rst

@@ -20,3 +20,4 @@ Theory and Methodology
     energy_deposition
     parallelization
     cmfd
+    random_ray

docs/source/methods/random_ray.rst

@@ -0,0 +1,804 @@
+.. _methods_random_ray:
+
+==========
+Random Ray
+==========
+
+.. _methods_random_ray_intro:
+
+-------------------
+What is Random Ray?
+-------------------
+
+`Random ray <Tramm-2017a>`_ is a stochastic transport method, closely related to
+the deterministic Method of Characteristics (MOC) [Askew-1972]_. Rather than
+each ray representing a single neutron as in Monte Carlo, it represents a
+characteristic line through the simulation geometry upon which the transport
+equation can be written as an ordinary differential equation that can be solved
+analytically (although with discretization required in energy, making it a
+multigroup method). The behavior of the governing transport equation can be
+approximated by solving along many characteristic tracks (rays) through the
+system. Unlike particles in Monte Carlo, rays in random ray or MOC are not
+affected by the material characteristics of the simulated problem---rays are
+selected so as to explore the full simulation problem with a statistically equal
+distribution in space and angle.
+
+.. raw:: html
+
+    <iframe width="560" height="315" src="https://www.youtube.com/embed/pHQq3FE4PDo?si=kPm9ngMBr95wLRGC" title="YouTube video player" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture; web-share" allowfullscreen></iframe>
+
+The above animation is an example of the random ray integration process at work,
+showing a series of random rays being sampled and transported through the
+geometry. In the following sections, we will discuss how the random ray solver
+works.
+
+----------------------------------------------
+Why is a Random Ray Solver Included in OpenMC?
+----------------------------------------------
+
+* One area that Monte Carlo struggles with is maintaining numerical efficiency
+  in regions of low physical particle flux. Random ray, on the other hand, has
+  approximately even variance throughout the entire global simulation domain,
+  such that areas with low neutron flux are no less well known that areas of
+  high neutron flux. Absent weight windows in MC, random ray can be several
+  orders of magnitude faster than multigroup Monte Carlo in classes of problems
+  where areas with low physical neutron flux need to be resolved. While MC
+  uncertainty can be greatly improved with variance reduction techniques, they
+  add some user complexity, and weight windows can often be expensive to
+  generate via MC transport alone (e.g., via the `MAGIC method
+  <https://doi.org/10.1016/j.fusengdes.2011.01.059>`_). The random ray solver
+  may be used in future versions of OpenMC as a fast way to generate weight
+  windows for subsequent usage by the MC solver in OpenMC.
+
+* In practical implementation terms, random ray is mechanically very similar to
+  how Monte Carlo works, in terms of the process of ray tracing on constructive
+  solid geometry (CSG) and handling stochastic convergence, etc. In the original
+  1972 paper by Askew that introduces MOC (which random ray is a variant of), he
+  stated:
+
+    .. epigraph::
+
+        "One of the features of the method proposed [MoC] is that ... the
+        tracking process needed to perform this operation is common to the
+        proposed method ... and to Monte Carlo methods. Thus a single tracking
+        routine capable of recognizing a geometric arrangement could be utilized
+        to service all types of solution, choice being made depending which was
+        more appropriate to the problem size and required accuracy."
+
+        -- Askew [Askew-1972]_
+
+    This prediction holds up---the additional requirements needed in OpenMC to
+    handle random ray transport turned out to be fairly small.
+
+* It amortizes the code complexity in OpenMC for representing multigroup cross
+  sections. There is a significant amount of interface code, documentation, and
+  complexity in allowing OpenMC to generate and use multigroup XS data in its
+  MGMC mode. Random ray allows the same multigroup data to be used, making full
+  reuse of these existing capabilities.
+
+-------------------------------
+Random Ray Numerical Derivation
+-------------------------------
+
+In this section, we will derive the numerical basis for the random ray solver
+mode in OpenMC. The derivation of random ray is also discussed in several papers
+(`1 <Tramm-2017a>`_, `2 <Tramm-2017b>`_, `3 <Tramm-2018>`_), and some of those
+derivations are reproduced here verbatim. Several extensions are also made to
+add clarity, particularly on the topic of OpenMC's treatment of cell volumes in
+the random ray solver.
+
+~~~~~~~~~~~~~~~~~~~~~~~~~
+Method of Characteristics
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The Boltzmann neutron transport equation is a partial differential equation
+(PDE) that describes the angular flux within a system. It is a balance equation,
+with the streaming and absorption terms typically appearing on the left hand
+side, which are balanced by the scattering source and fission source terms on
+the right hand side.
+
+.. math::
+    :label: transport
+
+    \begin{align*}
+    \mathbf{\Omega} \cdot \mathbf{\nabla} \psi(\mathbf{r},\mathbf{\Omega},E) & + \Sigma_t(\mathbf{r},E) \psi(\mathbf{r},\mathbf{\Omega},E) = \\
+    & \int_0^\infty d E^\prime \int_{4\pi} d \Omega^{\prime} \Sigma_s(\mathbf{r},\mathbf{\Omega}^\prime \rightarrow \mathbf{\Omega}, E^\prime \rightarrow E) \psi(\mathbf{r},\mathbf{\Omega}^\prime, E^\prime) \\
+    & + \frac{\chi(\mathbf{r}, E)}{4\pi k_{eff}} \int_0^\infty dE^\prime \nu \Sigma_f(\mathbf{r},E^\prime) \int_{4\pi}d \Omega^\prime \psi(\mathbf{r},\mathbf{\Omega}^\prime,E^\prime)
+    \end{align*}
+
+In Equation :eq:`transport`, :math:`\psi` is the angular neutron flux. This
+parameter represents the total distance traveled by all neutrons in a particular
+direction inside of a control volume per second, and is often given in units of
+:math:`1/(\text{cm}^{2} \text{s})`. As OpenMC does not support time dependence
+in the random ray solver mode, we consider the steady state equation, where the
+units of flux become :math:`1/\text{cm}^{2}`. The angular direction unit vector,
+:math:`\mathbf{\Omega}`, represents the direction of travel for the neutron. The
+spatial position vector, :math:`\mathbf{r}`,  represents the location within the
+simulation. The neutron energy, :math:`E`, or speed in continuous space, is
+often given in units of electron volts. The total macroscopic neutron cross
+section is :math:`\Sigma_t`. This value represents the total probability of
+interaction between a neutron traveling at a certain speed (i.e., neutron energy
+:math:`E`) and a target nucleus (i.e., the material through which the neutron is
+traveling) per unit path length, typically given in units of
+:math:`1/\text{cm}`. Macroscopic cross section data is a combination of
+empirical data and quantum mechanical modeling employed in order to generate an
+evaluation represented either in pointwise form or resonance parameters for each
+target isotope of interest in a material, as well as the density of the
+material, and is provided as input to a simulation. The scattering neutron cross
+section, :math:`\Sigma_s`, is similar to the total cross section but only
+measures scattering interactions between the neutron and the target nucleus, and
+depends on the change in angle and energy the neutron experiences as a result of
+the interaction. Several additional reactions like (n,2n) and (n,3n) are
+included in the scattering transfer cross section. The fission neutron cross
+section, :math:`\Sigma_f`, is also similar to the total cross section but only
+measures the fission interaction between a neutron and a target nucleus. The
+energy spectrum for neutrons born from fission, :math:`\chi`, represents a known
+distribution of outgoing neutron energies based on the material that fissioned,
+which is taken as input data to a computation. The average number of neutrons
+born per fission is :math:`\nu`. The eigenvalue of the equation,
+:math:`k_{eff}`, represents the effective neutron multiplication factor. If the
+right hand side of Equation :eq:`transport` is condensed into a single term,
+represented by the total neutron source term :math:`Q(\mathbf{r}, \mathbf{\Omega},E)`,
+the form given in Equation :eq:`transport_simple` is reached.
+
+.. math::
+    :label: transport_simple
+
+    \overbrace{\mathbf{\Omega} \cdot \mathbf{\nabla} \psi(\mathbf{r},\mathbf{\Omega},E)}^{\text{streaming term}} + \overbrace{\Sigma_t(\mathbf{r},E) \psi(\mathbf{r},\mathbf{\Omega},E)}^{\text{absorption term}} = \overbrace{Q(\mathbf{r}, \mathbf{\Omega},E)}^{\text{total neutron source term}}
+
+Fundamentally, MOC works by solving Equation :eq:`transport_simple` along a
+single characteristic line, thus altering the full spatial and angular scope of
+the transport equation into something that holds true only for a particular
+linear path (or track) through the reactor. These tracks are linear for neutral
+particles that are not subject to field effects. With our transport equation in
+hand, we will now derive the solution along a track. To accomplish this, we
+parameterize :math:`\mathbf{r}` with respect to some reference location
+:math:`\mathbf{r}_0` such that :math:`\mathbf{r} = \mathbf{r}_0 + s\mathbf{\Omega}`. In this
+manner, Equation :eq:`transport_simple` can be rewritten for a specific segment
+length :math:`s` at a specific angle :math:`\mathbf{\Omega}` through a constant
+cross section region of the reactor geometry as in Equation :eq:`char_long`.
+
+.. math::
+    :label: char_long
+
+    \mathbf{\Omega} \cdot \mathbf{\nabla} \psi(\mathbf{r}_0 + s\mathbf{\Omega},\mathbf{\Omega},E) + \Sigma_t(\mathbf{r}_0 + s\mathbf{\Omega},E) \psi(\mathbf{r}_0 + s\mathbf{\Omega},\mathbf{\Omega},E) = Q(\mathbf{r}_0 + s\mathbf{\Omega}, \mathbf{\Omega},E)
+
+As this equation holds along a one dimensional path, we can assume the
+dependence of :math:`s` on :math:`\mathbf{r}_0` and :math:`\mathbf{\Omega}` such that
+:math:`\mathbf{r}_0 + s\mathbf{\Omega}` simplifies to :math:`s`. When the differential
+operator is also applied to the angular flux :math:`\psi`, we arrive at the
+characteristic form of the Boltzmann Neutron Transport Equation given in
+Equation :eq:`char`.
+
+.. math::
+    :label: char
+
+    \frac{d}{ds} \psi(s,\mathbf{\Omega},E) + \Sigma_t(s,E) \psi(s,\mathbf{\Omega},E) = Q(s, \mathbf{\Omega},E)
+
+An analytical solution to this characteristic equation can be achieved with the
+use of an integrating factor:
+
+.. math::
+    :label: int_factor
+
+    e^{ \int_0^s ds' \Sigma_t (s', E)}
+
+to arrive at the final form of the characteristic equation shown in Equation
+:eq:`full_char`.
+
+.. math::
+    :label: full_char
+
+    \psi(s,\mathbf{\Omega},E) = \psi(\mathbf{r}_0,\mathbf{\Omega},E) e^{-\int_0^s ds^\prime \Sigma_t(s^\prime,E)} + \int_0^s ds^{\prime\prime} Q(s^{\prime\prime},\mathbf{\Omega}, E) e^{-\int_{s^{\prime\prime}}^s ds^\prime \Sigma_t(s^\prime,E)}
+
+With this characteristic form of the transport equation, we now have an
+analytical solution along a linear path through any constant cross section
+region of a system. While the solution only holds along a linear track, no
+discretizations have yet been made.
+
+Similar to many other solution approaches to the Boltzmann neutron transport
+equation, the MOC approach also uses a "multigroup" approximation in order to
+discretize the continuous energy spectrum of neutrons traveling through the
+system into fixed set of energy groups :math:`G`, where each group :math:`g \in
+G` has its own specific cross section parameters. This makes the difficult
+non-linear continuous energy dependence much more manageable as group wise cross
+section data can be precomputed and fed into a simulation as input data. The
+computation of multigroup cross section data is not a trivial task and can
+introduce errors in the simulation. However, this is an active field of research
+common to all multigroup methods, and there are numerous generation methods
+available that are capable of reducing the biases introduced by the multigroup
+approximation. Commonly used methods include the subgroup self-shielding method
+and use of fast (unconverged) Monte Carlo simulations to produce cross section
+estimates. It is important to note that Monte Carlo methods are capable of
+treating the energy variable of the neutron continuously, meaning that they do
+not need to make this approximation and are therefore not subject to any
+multigroup errors.
+
+Following the multigroup discretization, another assumption made is that a large
+and complex problem can be broken up into small constant cross section regions,
+and that these regions have group dependent, flat, isotropic sources (fission
+and scattering), :math:`Q_g`. Anisotropic as well as higher order sources are
+also possible with MOC-based methods but are not used yet in OpenMC for
+simplicity. With these key assumptions, the multigroup MOC form of the neutron
+transport equation can be written as in Equation :eq:`moc_final`.
+
+.. math::
+    :label: moc_final
+
+    \psi_g(s, \mathbf{\Omega}) = \psi_g(\mathbf{r_0}, \mathbf{\Omega}) e^{-\int_0^s ds^\prime \Sigma_{t_g}(s^\prime)} + \int_0^s ds^{\prime\prime} Q_g(s^{\prime\prime},\mathbf{\Omega}) e^{-\int_{s^{\prime\prime}}^s ds^\prime \Sigma_{t_g}(s^\prime)}
+
+The CSG definition of the system is used to create spatially defined source
+regions (each region being denoted as :math:`i`). These neutron source regions
+are often approximated as being constant
+(flat) in source intensity but can also be defined using a higher order source
+(linear, quadratic, etc.) that allows for fewer source regions to be required to
+achieve a specified solution fidelity. In OpenMC, the approximation of a
+spatially constant isotropic fission and scattering source :math:`Q_{i,g}` in
+cell :math:`i` leads
+to simple exponential attenuation along an individual characteristic of length
+:math:`s` given by Equation :eq:`fsr_attenuation`.
+
+.. math::
+    :label: fsr_attenuation
+
+    \psi_g(s) = \psi_g(0) e^{-\Sigma_{t,i,g} s} + \frac{Q_{i,g}}{\Sigma_{t,i,g}} \left( 1 - e^{-\Sigma_{t,i,g} s} \right)
+
+For convenience, we can also write this equation in terms of the incoming and
+outgoing angular flux (:math:`\psi_g^{in}` and :math:`\psi_g^{out}`), and
+consider a specific tracklength for a particular ray :math:`r` crossing cell
+:math:`i` as :math:`\ell_r`, as in:
+
+.. math::
+    :label: fsr_attenuation_in_out
+
+    \psi_g^{out} = \psi_g^{in} e^{-\Sigma_{t,i,g} \ell_r} + \frac{Q_{i,g}}{\Sigma_{t,i,g}} \left( 1 - e^{-\Sigma_{t,i,g} \ell_r} \right) .
+
+We can then define the average angular flux of a single ray passing through the
+cell as:
+
+.. math::
+    :label: average
+
+    \overline{\psi}_{r,i,g} = \frac{1}{\ell_r} \int_0^{\ell_r} \psi_{g}(s)ds .
+
+We can then substitute in Equation :eq:`fsr_attenuation` and solve, resulting
+in:
+
+.. math::
+    :label: average_solved
+
+    \overline{\psi}_{r,i,g} = \frac{Q_{i,g}}{\Sigma_{t,i,g}} - \frac{\psi_{r,g}^{out} - \psi_{r,g}^{in}}{\ell_r \Sigma_{t,i,g}} .
+
+By rearranging Equation :eq:`fsr_attenuation_in_out`, we can then define
+:math:`\Delta \psi_{r,g}` as the change in angular flux for ray :math:`r`
+passing through region :math:`i` as:
+
+.. math::
+    :label: delta_psi
+
+    \Delta \psi_{r,g} = \psi_{r,g}^{in} - \psi_{r,g}^{out} = \left(\psi_{r,g}^{in} - \frac{Q_{i,g}}{\Sigma_{t,i,g}} \right) \left( 1 - e^{-\Sigma_{t,i,g} \ell_r} \right) .
+
+Equation :eq:`delta_psi` is a useful expression as it is easily computed with
+the known inputs for a ray crossing through the region.
+
+By substituting :eq:`delta_psi` into :eq:`average_solved`, we can arrive at a
+final expression for the average angular flux for a ray crossing a region as:
+
+.. math::
+    :label: average_psi_final
+
+    \overline{\psi}_{r,i,g} = \frac{Q_{i,g}}{\Sigma_{t,i,g}} + \frac{\Delta \psi_{r,g}}{\ell_r \Sigma_{t,i,g}}
+
+~~~~~~~~~~~
+Random Rays
+~~~~~~~~~~~
+
+In the previous subsection, the governing characteristic equation along a 1D
+line through the system was written, such that an analytical solution for the
+ODE can be computed. If enough characteristic tracks (ODEs) are solved, then the
+behavior of the governing PDE can be numerically approximated. In traditional
+deterministic MOC, the selection of tracks is chosen deterministically, where
+azimuthal and polar quadratures are defined along with even track spacing in
+three dimensions. This is the point at which random ray diverges from
+deterministic MOC numerically. In the random ray method, rays are randomly
+sampled from a uniform distribution in space and angle and tracked along a
+predefined distance through the geometry before terminating. **Importantly,
+different rays are sampled each power iteration, leading to a fully stochastic
+convergence process.** This results in a need to utilize both inactive and
+active batches as in the Monte Carlo method.
+
+While Monte Carlo implicitly converges the scattering source fully within each
+iteration, random ray (and MOC) solvers are not typically written to fully
+converge the scattering source within a single iteration. Rather, both the
+fission and scattering sources are updated each power iteration, thus requiring
+enough outer iterations to reach a stationary distribution in both the fission
+source and scattering source. So, even in a low dominance ratio problem like a
+2D pincell, several hundred inactive batches may still be required with random
+ray to allow the scattering source to fully develop, as neutrons undergoing
+hundreds of scatters may constitute a non-trivial contribution to the fission
+source. We note that use of a two-level second iteration scheme is sometimes
+used by some MOC or random ray solvers so as to fully converge the scattering
+source with many inner iterations before updating the fission source in the
+outer iteration. It is typically more efficient to use the single level
+iteration scheme, as there is little reason to spend so much work converging the
+scattering source if the fission source is not yet converged.
+
+Overall, the difference in how random ray and Monte Carlo converge the
+scattering source means that in practice, random ray typically requires more
+inactive iterations than are required in Monte Carlo. While a Monte Carlo
+simulation may need 100 inactive iterations to reach a stationary source
+distribution for many problems, a random ray solve will likely require 1,000
+iterations or more. Source convergence metrics (e.g., Shannon entropy) are thus
+recommended when performing random ray simulations to ascertain when the source
+has fully developed.
+
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+Converting Angular Flux to Scalar Flux
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Thus far in our derivation, we have been able to write analytical equations that
+solve for the change in angular flux of a ray crossing a flat source region
+(Equation :eq:`delta_psi`) as well as the ray's average angular flux through
+that region (Equation :eq:`average_psi_final`). To determine the source for the
+next power iteration, we need to assemble our estimates of angular fluxes from
+all the sampled rays into scalar fluxes within each FSR.
+
+We can define the scalar flux in region :math:`i` as:
+
+.. math::
+    :label: integral
+
+    \phi_i = \frac{\int_{V_i} \int_{4\pi} \psi(r, \Omega) d\Omega d\mathbf{r}}{\int_{V_i} d\mathbf{r}} .
+
+The integral in the numerator:
+
+.. math::
+    :label: numerator
+
+    \int_{V_i} \int_{4\pi} \psi(r, \Omega) d\Omega d\mathbf{r} .
+
+is not known analytically, but with random ray, we are going the numerically
+approximate it by discretizing over a finite number of tracks (with a finite
+number of locations and angles) crossing the domain. We can then use the
+characteristic method to determine the total angular flux along that line.
+
+Conceptually, this can be thought of as taking a volume-weighted sum of angular
+fluxes for all :math:`N_i` rays that happen to pass through cell :math:`i` that
+iteration. When written in discretized form (with the discretization happening
+in terms of individual ray segments :math:`r` that pass through region
+:math:`i`), we arrive at:
+
+.. math::
+    :label: discretized
+
+    \phi_{i,g} = \frac{\int_{V_i} \int_{4\pi} \psi(r, \Omega) d\Omega d\mathbf{r}}{\int_{V_i} d\mathbf{r}} = \overline{\overline{\psi}}_{i,g} \approx \frac{\sum\limits_{r=1}^{N_i} \ell_r w_r \overline{\psi}_{r,i,g}}{\sum\limits_{r=1}^{N_i} \ell_r w_r} .
+
+Here we introduce the term :math:`w_r`, which represents the "weight" of the ray
+(its 2D area), such that the volume that a ray is responsible for can be
+determined by multiplying its length :math:`\ell` by its weight :math:`w`. As
+the scalar flux vector is a shape function only, we are actually free to
+multiply all ray weights :math:`w` by any constant such that the overall shape
+is still maintained, even if the magnitude of the shape function changes. Thus,
+we can simply set :math:`w_r` to be unity for all rays, such that:
+
+.. math::
+    :label: weights
+
+    \text{Volume of cell } i = V_i \approx \sum\limits_{r=1}^{N_i} \ell_r w_r = \sum\limits_{r=1}^{N_i} \ell_r .
+
+We can then rewrite our discretized equation as:
+
+.. math::
+    :label: discretized_2
+
+    \phi_{i,g} \approx \frac{\sum\limits_{r=1}^{N_i} \ell_r w_r \overline{\psi}_{r,i,g}}{\sum\limits_{r=1}^{N_i} \ell_r w_r} = \frac{\sum\limits_{r=1}^{N_i} \ell_r \overline{\psi}_{r,i,g}}{\sum\limits_{r=1}^{N_i} \ell_r} .
+
+Thus, the scalar flux can be inferred if we know the volume weighted sum of the
+average angular fluxes that pass through the cell. Substituting
+:eq:`average_psi_final` into :eq:`discretized_2`, we arrive at:
+
+.. math::
+    :label: scalar_full
+
+    \phi_{i,g} = \frac{\int_{V_i} \int_{4\pi} \psi(r, \Omega) d\Omega d\mathbf{r}}{\int_{V_i} d\mathbf{r}} = \overline{\overline{\psi}}_{i,g} = \frac{\sum\limits_{r=1}^{N_i} \ell_r \overline{\psi}_{r,i,g}}{\sum\limits_{r=1}^{N_i} \ell_r} = \frac{\sum\limits_{r=1}^{N_i} \ell_r \frac{Q_{i,g}}{\Sigma_{t,i,g}} + \frac{\Delta \psi_{r,g}}{\ell_r \Sigma_{t,i,g}}}{\sum\limits_{r=1}^{N_i} \ell_r},
+
+which when partially simplified becomes:
+
+.. math::
+    :label: scalar_four_vols
+
+    \phi =  \frac{Q_{i,g} \sum\limits_{r=1}^{N_i} \ell_r}{\Sigma_{t,i,g} \sum\limits_{r=1}^{N_i} \ell_r} + \frac{\sum\limits_{r=1}^{N_i} \ell_r \frac{\Delta \psi_i}{\ell_r}}{\Sigma_{t,i,g} \sum\limits_{r=1}^{N_i} \ell_r} .
+
+Note that there are now four (seemingly identical) volume terms in this equation.
+
+~~~~~~~~~~~~~~
+Volume Dilemma
+~~~~~~~~~~~~~~
+
+At first glance, Equation :eq:`scalar_four_vols` appears ripe for cancellation
+of terms. Mathematically, such cancellation allows us to arrive at the following
+"naive" estimator for the scalar flux:
+
+.. math::
+    :label: phi_naive
+
+    \phi_{i,g}^{naive} = \frac{Q_{i,g} }{\Sigma_{t,i,g}} + \frac{\sum\limits_{r=1}^{N_i} \Delta \psi_{r,g}}{\Sigma_{t,i,g} \sum\limits_{r=1}^{N_i} \ell_r} .
+
+This derivation appears mathematically sound at first glance but unfortunately
+raises a serious issue as discussed in more depth by `Tramm et al.
+<Tramm-2020>`_ and `Cosgrove and Tramm <Cosgrove-2023>`_. Namely, the second
+term:
+
+.. math::
+    :label: ratio_estimator
+
+     \frac{\sum\limits_{r=1}^{N_i} \Delta \psi_{r,g}}{\Sigma_{t,i,g} \sum\limits_{r=1}^{N_i} \ell_r}
+
+features stochastic variables (the sums over random ray lengths and angular
+fluxes) in both the numerator and denominator, making it a stochastic ratio
+estimator, which is inherently biased. In practice, usage of the naive estimator
+does result in a biased, but "consistent"  estimator (i.e., it is biased, but
+the bias tends towards zero as the sample size increases). Experimentally, the
+right answer can be obtained with this estimator, though a very fine ray density
+is required to eliminate the bias.
+
+How might we solve the biased ratio estimator problem? While there is no obvious
+way to alter the numerator term (which arises from the characteristic
+integration approach itself), there is potentially more flexibility in how we
+treat the stochastic term in the denominator, :math:`\sum\limits_{r=1}^{N_i}
+\ell_r` . From Equation :eq:`weights` we know that this term can be directly
+inferred from the volume of the problem, which does not actually change between
+iterations. Thus, an alternative treatment for this "volume" term in the
+denominator is to replace the actual stochastically sampled total track length
+with the expected value of the total track length. For instance, if the true
+volume of the FSR is known (as is the total volume of the full simulation domain
+and the total tracklength used for integration that iteration), then we know the
+true expected value of the tracklength in that FSR. That is, if a FSR accounts
+for 2% of the overall volume of a simulation domain, then we know that the
+expected value of tracklength in that FSR will be 2% of the total tracklength
+for all rays that iteration. This is a key insight, as it allows us to the
+replace the actual tracklength that was accumulated inside that FSR each
+iteration with the expected value.
+
+If we know the analytical volumes, then those can be used to directly compute
+the expected value of the tracklength in each cell. However, as the analytical
+volumes are not typically known in OpenMC due to the usage of user-defined
+constructive solid geometry, we need to source this quantity from elsewhere. An
+obvious choice is to simply accumulate the total tracklength through each FSR
+across all iterations (batches) and to use that sum to compute the expected
+average length per iteration, as:
+
+.. math::
+    :label: sim_estimator
+
+       \sum\limits^{}_{i} \ell_i \approx \frac{\sum\limits^{B}_{b}\sum\limits^{N_i}_{r} \ell_{b,r} }{B}
+
+where :math:`b` is a single batch in :math:`B` total batches simulated so far.
+
+In this manner, the expected value of the tracklength will become more refined
+as iterations continue, until after many iterations the variance of the
+denominator term becomes trivial compared to the numerator term, essentially
+eliminating the presence of the stochastic ratio estimator. A "simulation
+averaged" estimator is therefore:
+
+.. math::
+    :label: phi_sim
+
+    \phi_{i,g}^{simulation} = \frac{Q_{i,g} }{\Sigma_{t,i,g}} + \frac{\sum\limits_{r=1}^{N_i} \Delta \psi_{r,g}}{\Sigma_{t,i,g} \frac{\sum\limits^{B}_{b}\sum\limits^{N_i}_{r} \ell_{b,r} }{B}}
+
+In practical terms, the "simulation averaged" estimator is virtually
+indistinguishable numerically from use of the true analytical volume to estimate
+this term. Note also that the term "simulation averaged" refers only to the
+volume/length treatment, the scalar flux estimate itself is computed fully again
+each iteration.
+
+There are some drawbacks to this method. Recall, this denominator volume term
+originally stemmed from taking a volume weighted integral of the angular flux,
+in which case the denominator served as a normalization term for the numerator
+integral in Equation :eq:`integral`. Essentially, we have now used a different
+term for the volume in the numerator as compared to the normalizing volume in
+the denominator. The inevitable mismatch (due to noise) between these two
+quantities results in a significant increase in variance. Notably, the same
+problem occurs if using a tracklength estimate based on the analytical volume,
+as again the numerator integral and the normalizing denominator integral no
+longer match on a per-iteration basis.
+
+In practice, the simulation averaged method does completely remove the bias,
+though at the cost of a notable increase in variance. Empirical testing reveals
+that on most problems, the simulation averaged estimator does win out overall in
+numerical performance, as a much coarser quadrature can be used resulting in
+faster runtimes overall. Thus, OpenMC uses the simulation averaged estimator in
+its random ray mode.
+
+~~~~~~~~~~~~~~~
+Power Iteration
+~~~~~~~~~~~~~~~
+
+Given a starting source term, we now have a way of computing an estimate of the
+scalar flux in each cell by way of transporting rays randomly through the
+domain, recording the change in angular flux for the rays into each cell as they
+make their traversals, and summing these contributions up as in Equation
+:eq:`phi_sim`. How then do we turn this into an iterative process such that we
+improve the estimate of the source and scalar flux over many iterations, given
+that our initial starting source will just be a guess?
+
+The source :math:`Q^{n}` for iteration :math:`n` can be inferred
+from the scalar flux from the previous iteration :math:`n-1` as:
+
+.. math::
+    :label: source_update
+
+    Q^{n}(i, g) = \frac{\chi}{k^{n-1}_{eff}} \nu \Sigma_f(i, g) \phi^{n-1}(g) + \sum\limits^{G}_{g'} \Sigma_{s}(i,g,g') \phi^{n-1}(g')
+
+where :math:`Q^{n}(i, g)` is the total source (fission + scattering) in region
+:math:`i` and energy group :math:`g`. Notably, the in-scattering source in group
+:math:`g` must be computed by summing over the contributions from all groups
+:math:`g' \in G`.
+
+In a similar manner, the eigenvalue for iteration :math:`n` can be computed as:
+
+.. math::
+    :label: eigenvalue_update
+
+    k^{n}_{eff} = k^{n-1}_{eff} \frac{F^n}{F^{n-1}},
+
+where the total spatial- and energy-integrated fission rate :math:`F^n` in
+iteration :math:`n` can be computed as:
+
+.. math::
+    :label: fission_source
+
+    F^n = \sum\limits^{M}_{i} \left( V_i \sum\limits^{G}_{g} \nu \Sigma_f(i, g) \phi^{n}(g) \right)
+
+where :math:`M` is the total number of FSRs in the simulation. Similarly, the
+total spatial- and energy-integrated fission rate :math:`F^{n-1}` in iteration
+:math:`n-1` can be computed as:
+
+.. math::
+    :label: fission_source_prev
+
+    F^{n-1} = \sum\limits^{M}_{i} \left( V_i \sum\limits^{G}_{g} \nu \Sigma_f(i, g) \phi^{n-1}(g) \right)
+
+Notably, the volume term :math:`V_i` appears in the eigenvalue update equation.
+The same logic applies to the treatment of this term as was discussed earlier.
+In OpenMC, we use the "simulation averaged" volume derived from summing over all
+ray tracklength contributions to a FSR over all iterations and dividing by the
+total integration tracklength to date. Thus, Equation :eq:`fission_source`
+becomes:
+
+.. math::
+    :label: fission_source_volumed
+
+    F^n = \sum\limits^{M}_{i} \left( \frac{\sum\limits^{B}_{b}\sum\limits^{N_i}_{r} \ell_{b,r} }{B} \sum\limits^{G}_{g} \nu \Sigma_f(i, g) \phi^{n}(g) \right)
+
+and a similar substitution can be made to update Equation
+:eq:`fission_source_prev` . In OpenMC, the most up-to-date version of the volume
+estimate is used, such that the total fission source from the previous iteration
+(:math:`n-1`) is also recomputed each iteration.
+
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+Ray Starting Conditions and Inactive Length
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Another key area of divergence between deterministic MOC and random ray is the
+starting conditions for rays. In deterministic MOC, the angular flux spectrum
+for rays are stored at any reflective or periodic boundaries so as to provide a
+starting condition for the next iteration. As there are many tracks, storage of
+angular fluxes can become costly in terms of memory consumption unless there are
+only vacuum boundaries present.
+
+In random ray, as the starting locations of rays are sampled anew each
+iteration, the initial angular flux spectrum for the ray is unknown. While a
+guess can be made by taking the isotropic source from the FSR the ray was
+sampled in, direct usage of this quantity would result in significant bias and
+error being imparted on the simulation.
+
+Thus, an `on-the-fly approximation method <Tramm-2017a>`_ was developed (known
+as the "dead zone"), where the first several mean free paths of a ray are
+considered to be "inactive" or "read only". In this sense, the angular flux is
+solved for using the MOC equation, but the ray does not "tally" any scalar flux
+back to the FSRs that it travels through. After several mean free paths have
+been traversed, the ray's angular flux spectrum typically becomes dominated by
+the accumulated source terms from the cells it has traveled through, while the
+(incorrect) starting conditions have been attenuated away. In the animation in
+the :ref:`introductory section on this page <methods_random_ray_intro>`, the
+yellow portion of the ray lengths is the dead zone. As can be seen in this
+animation, the tallied :math:`\sum\limits_{r=1}^{N_i} \Delta \psi_{r,g}` term
+that is plotted is not affected by the ray when the ray is within its inactive
+length. Only when the ray enters its active mode does the ray contribute to the
+:math:`\sum\limits_{r=1}^{N_i} \Delta \psi_{r,g}` sum for the iteration.
+
+~~~~~~~~~~~~~~~~~~~~~
+Ray Ending Conditions
+~~~~~~~~~~~~~~~~~~~~~
+
+To ensure that a uniform density of rays is integrated in space and angle
+throughout the simulation domain, after exiting the initial inactive "dead zone"
+portion of the ray, the rays are run for a user-specified distance. Typically, a
+choice of at least several times the length of the inactive "dead zone" is made
+so as to amortize the cost of the dead zone. For example, if a dead zone of 30
+cm is selected, then an active length of 300 cm might be selected so that the
+cost of the dead zone is at most 10% of the overall runtime.
+
+--------------------
+Simplified Algorithm
+--------------------
+
+A simplified set of functions that execute a single random ray power iteration
+are given below. Not all global variables are defined in this illustrative
+example, but the high level components of the algorithm are shown. A number of
+significant simplifications are made for clarity---for example, no inactive
+"dead zone" length is shown, geometry operations are abstracted, no parallelism
+(or thread safety) is expressed, a naive exponential treatment is used, and rays
+are not halted at their exact termination distances, among other subtleties.
+Nonetheless, the below algorithms may be useful for gaining intuition on the
+basic components of the random ray process. Rather than expressing the algorithm
+in abstract pseudocode, C++ is used to make the control flow easier to
+understand.
+
+The first block below shows the logic for a single power iteration (batch):
+
+.. code-block:: C++
+
+    double power_iteration(double k_eff) {
+
+        // Update source term (scattering + fission)
+        update_neutron_source(k_eff);
+
+        // Reset scalar fluxes to zero
+        fill<float>(global::scalar_flux_new, 0.0f);
+
+        // Transport sweep over all random rays for the iteration
+        for (int i = 0; i < nrays; i++) {
+            RandomRay ray;
+            initialize_ray(ray);
+            transport_single_ray(ray);
+        }
+
+        // Normalize scalar flux and update volumes
+        normalize_scalar_flux_and_volumes();
+
+        // Add source to scalar flux, compute number of FSR hits
+        add_source_to_scalar_flux();
+
+        // Compute k-eff using updated scalar flux
+        k_eff = compute_k_eff(k_eff);
+
+        // Set phi_old = phi_new
+        global::scalar_flux_old.swap(global::scalar_flux_new);
+
+        return k_eff;
+    }
+
+The second function shows the logic for transporting a single ray within the
+transport loop:
+
+.. code-block:: C++
+
+    void transport_single_ray(RandomRay& ray) {
+
+        // Reset distance to zero
+        double distance = 0.0;
+
+        // Continue transport of ray until active length is reached
+        while (distance < user_setting::active_length) {
+            // Ray trace to find distance to next surface (i.e., segment length)
+            double s = distance_to_nearest_boundary(ray);
+
+            // Attenuate flux (and accumulate source/attenuate) on segment
+            attenuate_flux(ray, s);
+
+            // Advance particle to next surface
+            ray.location = ray.location + s * ray.direction;
+
+            // Move ray across the surface
+            cross_surface(ray);
+
+            // Add segment length "s" to total distance traveled
+            distance += s;
+        }
+    }
+
+The final function below shows the logic for solving for the characteristic MOC
+equation (and accumulating the scalar flux contribution of the ray into the
+scalar flux value for the FSR).
+
+.. code-block:: C++
+
+    void attenuate_flux(RandomRay& ray, double s) {
+
+        // Determine which flat source region (FSR) the ray is currently in
+        int fsr = get_fsr_id(ray.location);
+
+        // Determine material type
+        int material = get_material_type(fsr);
+
+        // MOC incoming flux attenuation + source contribution/attenuation equation
+         for (int e = 0; e < global::n_energy_groups; e++) {
+            float sigma_t = global::macro_xs[material].total;
+            float tau = sigma_t * s;
+            float delta_psi = (ray.angular_flux[e] - global::source[fsr][e] / sigma_t) * (1 - exp(-tau));
+            ray.angular_flux_[e] -= delta_psi;
+            global::scalar_flux_new[fsr][e] += delta_psi;
+        }
+
+        // Record total tracklength in this FSR (to compute volume)
+        global::volume[fsr] += s;
+    }
+
+------------------------
+How are Tallies Handled?
+------------------------
+
+Most tallies, filters, and scores that you would expect to work with a
+multigroup solver like random ray should work. For example, you can define 3D
+mesh tallies with energy filters and flux, fission, and nu-fission scores, etc.
+There are some restrictions though. For starters, it is assumed that all filter
+mesh boundaries will conform to physical surface boundaries (or lattice
+boundaries) in the simulation geometry. It is acceptable for multiple cells
+(FSRs) to be contained within a filter mesh cell (e.g., pincell-level or
+assembly-level tallies should work), but it is currently left as undefined
+behavior if a single simulation cell is able to score to multiple filter mesh
+cells. In the future, the capability to fully support mesh tallies may be added
+to OpenMC, but for now this restriction needs to be respected.
+
+---------------------------
+Fundamental Sources of Bias
+---------------------------
+
+Compared to continuous energy Monte Carlo simulations, the known sources of bias
+in random ray particle transport are:
+
+    - **Multigroup Energy Discretization:** The multigroup treatment of flux and
+      cross sections incurs a significant bias, as a reaction rate (:math:`R_g =
+      V \phi_g \Sigma_g`) for an energy group :math:`g` can only be conserved
+      for a given choice of multigroup cross section :math:`\Sigma_g` if the
+      flux (:math:`\phi_g`) is known a priori. If the flux was already known,
+      then there would be no point to the simulation, resulting in a fundamental
+      need for approximating this quantity. There are numerous methods for
+      generating relatively accurate multigroup cross section libraries that can
+      each be applied to a narrow design area reliably, although there are
+      always limitations and/or complexities that arise with a multigroup energy
+      treatment. This is by far the most significant source of simulation bias
+      between Monte Carlo and random ray for most problems. While the other
+      areas typically have solutions that are highly effective at mitigating
+      bias, error stemming from multigroup energy discretization is much harder
+      to remedy.
+    - **Flat Source Approximation:**. In OpenMC, a "flat" (0th order) source
+      approximation is made, wherein the scattering and fission sources within a
+      cell are assumed to be spatially uniform. As the source in reality is a
+      continuous function, this leads to bias, although the bias can be reduced
+      to acceptable levels if the flat source regions are sufficiently small.
+      The bias can also be mitigated by assuming a higher-order source (e.g.,
+      linear or quadratic), although OpenMC does not yet have this capability.
+      In practical terms, this source of bias can become very large if cells are
+      large (with dimensions beyond that of a typical particle mean free path),
+      but the subdivision of cells can often reduce this bias to trivial levels.
+    - **Anisotropic Source Approximation:** In OpenMC, the source is not only
+      assumed to be flat but also isotropic, leading to bias. It is possible for
+      MOC (and likely random ray) to treat anisotropy explicitly, but this is
+      not currently supported in OpenMC. This source of bias is not significant
+      for some problems, but becomes more problematic for others. Even in the
+      absence of explicit treatment of anistropy, use of transport-corrected
+      multigroup cross sections can often mitigate this bias, particularly for
+      light water reactor simulation problems.
+    - **Angular Flux Initial Conditions:** Each time a ray is sampled, its
+      starting angular flux is unknown, so a guess must be made (typically the
+      source term for the cell it starts in). Usage of an adequate inactive ray
+      length (dead zone) mitigates this error. As the starting guess is
+      attenuated at a rate of :math:`\exp(-\Sigma_t \ell)`, this bias can driven
+      below machine precision in a low cost manner on many problems.
+
+.. _Tramm-2017a: https://doi.org/10.1016/j.jcp.2017.04.038
+.. _Tramm-2017b: https://doi.org/10.1016/j.anucene.2017.10.015
+.. _Tramm-2018: https://dspace.mit.edu/handle/1721.1/119038
+.. _Tramm-2020: https://doi.org/10.1051/EPJCONF/202124703021
+.. _Cosgrove-2023: https://doi.org/10.1080/00295639.2023.2270618
+
+.. only:: html
+
+   .. rubric:: References
+
+.. [Askew-1972] Askew, “A Characteristics Formulation of the Neutron Transport
+    Equation in Complicated Geometries.” Technical Report AAEW-M 1108, UK Atomic
+    Energy Establishment (1972).

docs/source/pythonapi/base.rst

@@ -128,6 +128,7 @@ Constructing Tallies
    openmc.CollisionFilter
    openmc.SurfaceFilter
    openmc.MeshFilter
+   openmc.MeshBornFilter
    openmc.MeshSurfaceFilter
    openmc.EnergyFilter
    openmc.EnergyoutFilter

docs/source/pythonapi/capi.rst

@@ -16,7 +16,6 @@ Functions
    current_batch
    export_properties
    export_weight_windows
-   import_weight_windows
    finalize
    find_cell
    find_material
@@ -25,6 +24,7 @@ Functions
    hard_reset
    id_map
    import_properties
+   import_weight_windows
    init
    is_statepoint_batch
    iter_batches
@@ -40,9 +40,10 @@ Functions
    run
    run_in_memory
    sample_external_source
-   simulation_init
    simulation_finalize
+   simulation_init
    source_bank
+   statepoint_load
    statepoint_write
 
 Classes
@@ -53,15 +54,86 @@ Classes
    :nosignatures:
    :template: myclass.rst
 
+   AzimuthalFilter
    Cell
+   CellFilter
+   CellInstanceFilter
+   CellbornFilter
+   CellfromFilter
+   CollisionFilter
    CylindricalMesh
+   DelayedGroupFilter
+   DistribcellFilter
    EnergyFilter
-   MaterialFilter
+   EnergyFunctionFilter
+   EnergyoutFilter
+   Filter
+   LegendreFilter
    Material
+   MaterialFilter
+   MaterialFromFilter
+   Mesh
    MeshFilter
+   MeshBornFilter
    MeshSurfaceFilter
+   MuFilter
    Nuclide
+   ParticleFilter
+   PolarFilter
    RectilinearMesh
    RegularMesh
+   SpatialLegendreFilter
+   SphericalHarmonicsFilter
    SphericalMesh
+   SurfaceFilter
    Tally
+   UniverseFilter
+   UnstructuredMesh
+   WeightWindows
+   ZernikeFilter
+   ZernikeRadialFilter
+
+Data
+----
+
+.. data:: cells
+
+   Mapping of cell ID to :class:`openmc.lib.Cell` instances.
+
+   :type: dict
+
+.. data:: filters
+
+   Mapping of filter ID to :class:`openmc.lib.Filter` instances.
+
+   :type: dict
+
+.. data:: materials
+
+   Mapping of material ID to :class:`openmc.lib.Material` instances.
+
+   :type: dict
+
+.. data:: meshes
+
+   Mapping of mesh ID to :class:`openmc.lib.Mesh` instances.
+
+   :type: dict
+
+.. data:: nuclides
+
+   Mapping of nuclide name to :class:`openmc.lib.Nuclide` instances.
+
+   :type: dict
+
+.. data:: tallies
+
+   Mapping of tally ID to :class:`openmc.lib.Tally` instances.
+
+   :type: dict
+
+.. data:: weight_windows
+
+   Mapping of weight window ID to :class:`openmc.lib.WeightWindows` instances.
+
+   :type: dict

docs/source/quickinstall.rst

@@ -107,31 +107,54 @@ can be used to access the installed packages.
 .. _Spack: https://spack.readthedocs.io/en/latest/
 .. _setup guide: https://spack.readthedocs.io/en/latest/getting_started.html
 
---------------------------------
-Installing from Source on Ubuntu
---------------------------------
+-------------------------------
+Manually Installing from Source
+-------------------------------
 
-To build OpenMC from source, several :ref:`prerequisites <prerequisites>` are
-needed. If you are using Ubuntu or higher, all prerequisites can be installed
-directly from the package manager:
+Obtaining prerequisites on Ubuntu
+---------------------------------
+
+When building OpenMC from source, all :ref:`prerequisites <prerequisites>` can
+be installed using the package manager:
 
 .. code-block:: sh
 
     sudo apt install g++ cmake libhdf5-dev libpng-dev
 
-After the packages have been installed, follow the instructions below for
-building and installing OpenMC from source.
+After the packages have been installed, follow the instructions to build from
+source below.
 
--------------------------------------------
-Installing from Source on Linux or Mac OS X
--------------------------------------------
+Obtaining prerequisites on macOS
+--------------------------------
+
+For an OpenMC build with multithreading enabled, a package manager like
+`Homebrew <https://brew.sh>`_ should first be installed. Then, the following
+packages should be installed, for example in Homebrew via:
+
+.. code-block:: sh
+
+   brew install llvm cmake xtensor hdf5 python libomp libpng
+
+The compiler provided by the above LLVM package should be used in place of the
+one provisioned by XCode, which does not support the multithreading library used
+by OpenMC. Consequently, the C++ compiler should explicitly be set before
+proceeding:
+
+.. code-block:: sh
+
+   export CXX=/opt/homebrew/opt/llvm/bin/clang++
+
+After the packages have been installed, follow the instructions to build from
+source below.
+
+Building Source on Linux or macOS
+---------------------------------
 
 All OpenMC source code is hosted on `GitHub
 <https://github.com/openmc-dev/openmc>`_. If you have `git
-<https://git-scm.com>`_, the `gcc <https://gcc.gnu.org/>`_ compiler suite,
-`CMake <https://cmake.org>`_, and `HDF5
-<https://www.hdfgroup.org/solutions/hdf5/>`_ installed, you can download and
-install OpenMC be entering the following commands in a terminal:
+<https://git-scm.com>`_, a modern C++ compiler, `CMake <https://cmake.org>`_,
+and `HDF5 <https://www.hdfgroup.org/solutions/hdf5/>`_ installed, you can
+download and install OpenMC by entering the following commands in a terminal:
 
 .. code-block:: sh
 
@@ -151,14 +174,14 @@ should specify an installation directory where you have write access, e.g.
     cmake -DCMAKE_INSTALL_PREFIX=$HOME/.local ..
 
 The :mod:`openmc` Python package must be installed separately. The easiest way
-to install it is using `pip <https://pip.pypa.io/en/stable/>`_, which is
-included by default in Python 3.4+. From the root directory of the OpenMC
-distribution/repository, run:
+to install it is using `pip <https://pip.pypa.io/en/stable/>`_.
+From the root directory of the OpenMC repository, run:
 
 .. code-block:: sh
 
     python -m pip install .
 
-If you want to build a parallel version of OpenMC (using OpenMP or MPI),
-directions can be found in the :ref:`detailed installation instructions
+By default, OpenMC will be built with multithreading support. To build
+distributed-memory parallel versions of OpenMC using MPI or to configure other
+options, directions can be found in the :ref:`detailed installation instructions
 <usersguide_build>`.

docs/source/usersguide/data.rst

@@ -279,6 +279,8 @@ The `official ENDF/B-VII.1 HDF5 library
 multipole library, so if you are using this library, the windowed multipole data
 will already be available to you.
 
+.. _create_mgxs:
+
 -------------------------
 Multigroup Cross Sections
 -------------------------

docs/source/usersguide/index.rst

@@ -25,4 +25,6 @@ essential aspects of using OpenMC to perform simulations.
     processing
     parallel
     volume
+    random_ray
     troubleshoot
+    

docs/source/usersguide/install.rst

@@ -464,11 +464,11 @@ can typically be set for a single command, i.e.
 
 .. _compile_linux:
 
-Compiling on Linux and Mac OS X
--------------------------------
+Compiling on Linux and macOS
+----------------------------
 
-To compile OpenMC on Linux or Max OS X, run the following commands from within
-the root directory of the source code:
+To compile OpenMC on Linux or macOS, run the following commands from within the
+root directory of the source code:
 
 .. code-block:: sh
 
@@ -540,7 +540,7 @@ to install the Python package in :ref:`"editable" mode <devguide_editable>`.
 Prerequisites
 -------------
 
-The Python API works with Python 3.7+. In addition to Python itself, the API
+The Python API works with Python 3.8+. In addition to Python itself, the API
 relies on a number of third-party packages. All prerequisites can be installed
 using Conda_ (recommended), pip_, or through the package manager in most Linux
 distributions.

docs/source/usersguide/random_ray.rst

@@ -0,0 +1,480 @@
+.. _random_ray:
+
+=================
+Random Ray Solver
+=================
+
+In general, the random ray solver mode uses most of the same settings and
+:ref:`run strategies <usersguide_particles>` as the standard Monte Carlo solver
+mode. For instance, random ray solves are also split up into :ref:`inactive and
+active batches <usersguide_batches>`. However, there are a couple of settings
+that are unique to the random ray solver and a few areas that the random ray
+run strategy differs, both of which will be described in this section.
+
+------------------------
+Enabling Random Ray Mode
+------------------------
+
+To utilize the random ray solver, the :attr:`~openmc.Settings.random_ray`
+dictionary must be present in the :class:`openmc.Settings` Python class. There
+are a number of additional settings that must be specified within this
+dictionary that will be discussed below. Additionally, the "multi-group" energy
+mode must be specified.
+
+-------
+Batches
+-------
+
+In Monte Carlo simulations, inactive batches are used to let the fission source
+develop into a stationary distribution before active batches are performed that
+actually accumulate statistics. While this is true of random ray as well, in the
+random ray mode the inactive batches are also used to let the scattering source
+develop. Monte Carlo fully represents the scattering source within each
+iteration (by its nature of fully simulating particles from birth to death
+through any number of physical scattering events), whereas the scattering source
+in random ray can only represent as many scattering events as batches have been
+completed. For example, by iteration 10 in random ray, the scattering source
+only captures the behavior of neutrons through their 10th scattering event.
+Thus, while inactive batches are only required in an eigenvalue solve in Monte
+Carlo, **inactive batches are required for both eigenvalue and fixed source
+solves in random ray mode** due to this additional need to converge the
+scattering source.
+
+The additional burden of converging the scattering source generally results in a
+higher requirement for the number of inactive batches---often by an order of
+magnitude or more. For instance, it may be reasonable to only use 50 inactive
+batches for a light water reactor simulation with Monte Carlo, but random ray
+might require 500 or more inactive batches. Similar to Monte Carlo,
+:ref:`Shannon entropy <usersguide_entropy>` can be used to gauge whether the
+combined scattering and fission source has fully developed.
+
+Similar to Monte Carlo, active batches are used in the random ray solver mode to
+accumulate and converge statistics on unknown quantities (i.e., the random ray
+sources, scalar fluxes, as well as any user-specified tallies).
+
+The batch parameters are set in the same manner as with the regular Monte Carlo
+solver::
+
+    settings = openmc.Settings()
+    settings.energy_mode = "multi-group"
+    settings.batches = 1200
+    settings.inactive = 600
+
+-------------------------------
+Inactive Ray Length (Dead Zone)
+-------------------------------
+
+A major issue with random ray is that the starting angular flux distribution for
+each sampled ray is unknown. Thus, an on-the-fly method is used to build a high
+quality approximation of the angular flux of the ray each iteration. This is
+accomplished by running the ray through an inactive length (also known as a dead
+zone length), where the ray is moved through the geometry and its angular flux
+is solved for via the normal :ref:`MOC <methods_random_ray_intro>` equation, but
+no information is written back to the system. Thus, the ray is run in a "read
+only" mode for the set inactive length. This parameter can be adjusted, in units
+of cm, as::
+
+    settings.random_ray['distance_inactive'] = 40.0
+
+After several mean free paths are traversed, the angular flux spectrum of the
+ray becomes dominated by the in-scattering and fission source components that it
+picked up when travelling through the geometry, while its original (incorrect)
+starting angular flux is attenuated toward zero. Thus, longer selections of
+inactive ray length will asymptotically approach the true angular flux.
+
+In practice, 10 mean free paths are sufficient (with light water reactors often
+requiring only about 10--50 cm of inactive ray length for the error to become
+undetectable). However, we caution that certain models with large quantities of
+void regions (even if just limited to a few streaming channels) may require
+significantly longer inactive ray lengths to ensure that the angular flux is
+accurate before the conclusion of the inactive ray length. Additionally,
+problems where a sensitive estimate of the uncollided flux is required (e.g.,
+the detector response to fast neutrons is required, and the detected is located
+far away from the source in a moderator region) may require the user to specify
+an inactive length that is derived from the pyhsical geometry of the simulation
+problem rather than its material properties. For instance, consider a detector
+placed 30 cm outside of a reactor core, with a moderator region separating the
+detector from the core. In this case, rays sampled in the moderator region and
+heading toward the detector will begin life with a highly scattered thermal
+spectrum and will have an inaccurate fast spectrum. If the dead zone length is
+only 20 cm, we might imagine such rays writing to the detector tally within
+their active lengths, despite their innaccurate estimate of the uncollided fast
+angular flux. Thus, an inactive length of 100--200 cm would ensure that any such
+rays would still be within their inactive regions, and only rays that have
+actually traversed through the core (and thus have an accurate representation of
+the core's emitted fast flux) will score to the detector region while in their
+active phase.
+
+
+------------------------------------
+Active Ray Length and Number of Rays
+------------------------------------
+
+Once the inactive length of the ray has completed, the active region of the ray
+begins. The ray is now run in regular mode, where changes in angular flux as it
+traverses through each flat source region are written back to the system, so as
+to contribute to the estimate for the iteration scalar flux (which is used to
+compute the source for the next iteration). The active ray length can be
+adjusted, in units of [cm], as::
+
+    settings.random_ray['distance_active'] = 400.0
+
+Assuming that a sufficient inactive ray length is used so that the starting
+angular flux is highly accurate, any selection of active length greater than
+zero is theoretically acceptable. However, in order to adequately sample the
+full integration domain, a selection of a very short track length would require
+a very high number of rays to be selected. Due to the static costs per ray of
+computing the starting angular flux in the dead zone, typically very short ray
+lengths are undesireable. Thus, to amortize the per-ray cost of the inactive
+region of the ray, it is desirable to select a very long inactive ray length.
+For example, if the inactive length is set to 20 cm, a 200 cm active ray length
+ensures that only about 10% of the overall simulation runtime is spent in the
+inactive ray phase integration, making the dead zone a relatively inexpensive
+way of estimating the angular flux.
+
+Thus, to fully amortize the cost of the dead zone integration, one might ask why
+not simply run a single ray per iteration with an extremely long active length?
+While this is also theoretically possible, this results in two issues. The first
+problem is that each ray only represents a single angular sample. As we want to
+sample the angular phase space of the simulation with similar fidelity to the
+spatial phase space, we naturally want a lot of angles. This means in practice,
+we want to balance the need to amortize the cost of the inactive region of the
+ray with the need to sample lots of angles. The second problem is that
+parallelism in OpenMC is expressed in terms of rays, with each being processed
+by an independent MPI rank and/or OpenMP thread, thus we want to ensure each
+thread has many rays to process.
+
+In practical terms, the best strategy is typically to set an active ray length
+that is about 10 times that of the inactive ray length. This is often the right
+balance between ensuring not too much time is spent in the dead zone, while
+still adequately sampling the angular phase space. However, as discussed in the
+previous section, some types of simulation may demand that additional thought be
+applied to this parameter. For instance, in the same example where we have a
+detector region far outside a reactor core, we want to make sure that there is
+enough active ray length that rays exiting the core can reach the detector
+region. For example, if the detector were to be 30 cm outside of the core, then
+we would need to ensure that at least a few hundred cm of active length were
+used so as to ensure even rays with indirect angles will be able to reach the
+target region.
+
+The number of rays each iteration can be set by reusing the normal Monte Carlo
+particle count selection parameter, as::
+
+    settings.particles = 2000
+
+-----------
+Ray Density
+-----------
+
+In the preceding sections, it was argued that for most use cases, the inactive
+length for a ray can be determined by taking a multiple of the mean free path
+for the limiting energy group. The active ray length could then be set by taking
+a multiple of the inactive length. With these parameters set, how many rays per
+iteration should be run?
+
+There are three basic settings that control the density of the stochastic
+quadrature being used to integrate the domain each iteration. These three
+variables are:
+
+- The number of rays (in OpenMC settings parlance, "particles")
+- The inactive distance per ray
+- The active distance per ray
+
+While the inactive and active ray lengths can usually be chosen by simply
+examining the geometry, tallies, and cross section data, one has much more
+flexibility in the choice of the number of rays to run. Consider a few
+scenarios:
+
+- If a choice of zero rays is made, then no information is gained by the system
+  after each batch.
+- If a choice of rays close to zero is made, then some information is gained
+  after each batch, but many source regions may not have been visited that
+  iteration, which is not ideal numerically and can result in instability.
+  Empirically, we have found that the simulation can remain stable and produce
+  accurate results even when on average 20% or more of the cells have zero rays
+  passing through them each iteration. However, besides the cost of transporting
+  rays, a new neutron source must be computed based on the scalar flux at each
+  iteration. This cost is dictated only by the number of source regions and
+  energy groups---it is independent of the number of rays. Thus, in practical
+  terms, if too few rays are run, then the simulation runtime becomes dominated
+  by the fixed cost of source updates, making it inefficient overall given that
+  a huge number of active batches will likely be required to converge statistics
+  to acceptable levels. Additionally, if many cells are missed each iteration,
+  then the fission and scattering sources may not develop very quickly,
+  resulting in a need for far more inactive batches than might otherwise be
+  required.
+- If a choice of running a very large number of rays is made such that you
+  guarantee that all cells are hit each iteration, this avoids any issues with
+  numerical instability. As even more rays are run, this reduces the number of
+  active batches that must be used to converge statistics and therefore
+  minimizes the fixed per-iteration source update costs. While this seems
+  advantageous, it has the same practical downside as with Monte Carlo---namely,
+  that the inactive batches tend to be overly well integrated, resulting in a
+  lot of wasted time. This issue is actually much more serious than in Monte
+  Carlo (where typically only tens of inactive batches are needed), as random
+  ray often requires hundreds or even thousands of inactive batches. Thus,
+  minimizing the cost of the source updates in the active phase needs to be
+  balanced against the increased cost of the inactive phase of the simulation.
+- If a choice of rays is made such that relatively few (e.g., around 0.1%) of
+  cells are missed each iteration, the cost of the inactive batches of the
+  simulation is minimized. In this "goldilocks" regime, there is very little
+  chance of numerical instability, and enough information is gained by each cell
+  to progress the fission and scattering sources forward at their maximum rate.
+  However, the inactive batches can proceed with minimal cost. While this will
+  result in the active phase of the simulation requiring more batches (and
+  correspondingly higher source update costs), the added cost is typically far
+  less than the savings by making the inactive phase much cheaper.
+
+To help you set this parameter, OpenMC will report the average flat source
+region miss rate at the end of the simulation. Additionally, OpenMC will alert
+you if very high miss rates are detected, indicating that more rays and/or a
+longer active ray length might improve numerical performance. Thus, a "guess and
+check" approach to this parameter is recommended, where a very low guess is
+made, a few iterations are performed, and then the simulation is restarted with
+a larger value until the "low ray density" messages go away.
+
+.. note::
+    In summary, the user should select an inactive length corresponding to many
+    times the mean free path of a particle, generally O(10--100) cm, to ensure accuracy of
+    the starting angular flux. The active length should be 10× the inactive
+    length to amortize its cost. The number of rays should be enough so that
+    nearly all :ref:`FSRs <subdivision_fsr>` are hit at least once each power iteration (the hit fraction
+    is reported by OpenMC for empirical user adjustment).
+
+.. warning::
+    For simulations where long range uncollided flux estimates need to be
+    accurately resolved (e.g., shielding, detector response, and problems with
+    significant void areas), make sure that selections for inactive and active
+    ray lengths are sufficiently long to allow for transport to occur between
+    source and target regions of interest.
+
+----------
+Ray Source
+----------
+
+Random ray requires that the ray source be uniform in space and isotropic in
+angle. To facilitate sampling, the user must specify a single random ray source
+for sampling rays in both eigenvalue and fixed source solver modes. The random
+ray integration source should be of type :class:`openmc.IndependentSource`, and
+is specified as part of the :attr:`openmc.Settings.random_ray` dictionary. Note
+that the source must not be limited to only fissionable regions. Additionally,
+the source box must cover the entire simulation domain. In the case of a
+simulation domain that is not box shaped, a box source should still be used to
+bound the domain but with the source limited to rejection sampling the actual
+simulation universe (which can be specified via the ``domains`` field of the
+:class:`openmc.IndependentSource` Python class). Similar to Monte Carlo sources,
+for two-dimensional problems (e.g., a 2D pincell) it is desirable to make the
+source bounded near the origin of the infinite dimension. An example of an
+acceptable ray source for a two-dimensional 2x2 lattice would look like:
+
+::
+
+    pitch = 1.26
+    lower_left  = (-pitch, -pitch, -pitch)
+    upper_right = ( pitch,  pitch,  pitch)
+    uniform_dist = openmc.stats.Box(lower_left, upper_right)
+    settings.random_ray['ray_source'] = openmc.IndependentSource(space=uniform_dist)
+
+.. note::
+    The random ray source is not related to the underlying particle flux or
+    source distribution of the simulation problem. It is akin to the selection
+    of an integration quadrature. Thus, in fixed source mode, the ray source
+    still needs to be provided and still needs to be uniform in space and angle
+    throughout the simulation domain. In fixed source mode, the user will
+    provide physical particle fixed sources in addition to the random ray
+    source.
+
+.. _subdivision_fsr:
+
+----------------------------------
+Subdivision of Flat Source Regions
+----------------------------------
+
+While the scattering and fission sources in Monte Carlo
+are treated continuously, they are assumed to be invariant (flat) within a
+MOC or random ray flat source region (FSR). This introduces bias into the
+simulation, which can be remedied by reducing the physical size of the FSR
+to dimensions below that of typical mean free paths of particles.
+
+In OpenMC, this subdivision currently must be done manually. The level of
+subdivision needed will be dependent on the fidelity the user requires. For
+typical light water reactor analysis, consider the following example subdivision
+of a two-dimensional 2x2 reflective pincell lattice:
+
+.. figure:: ../_images/2x2_materials.jpeg
+    :class: with-border
+    :width: 400
+
+    Material definition for an asymmetrical 2x2 lattice (1.26 cm pitch)
+
+.. figure:: ../_images/2x2_fsrs.jpeg
+    :class: with-border
+    :width: 400
+
+    FSR decomposition for an asymmetrical 2x2 lattice (1.26 cm pitch)
+
+In the future, automated subdivision of FSRs via mesh overlay may be supported.
+
+-------
+Tallies
+-------
+
+Most tallies, filters, and scores that you would expect to work with a
+multigroup solver like random ray are supported. For example, you can define 3D
+mesh tallies with energy filters and flux, fission, and nu-fission scores, etc.
+There are some restrictions though. For starters, it is assumed that all filter
+mesh boundaries will conform to physical surface boundaries (or lattice
+boundaries) in the simulation geometry. It is acceptable for multiple cells
+(FSRs) to be contained within a mesh element (e.g., pincell-level or
+assembly-level tallies should work), but it is currently left as undefined
+behavior if a single simulation cell is contained in multiple mesh elements.
+
+Supported scores:
+    - flux
+    - total
+    - fission
+    - nu-fission
+    - events
+
+Supported Estimators:
+    - tracklength
+
+Supported Filters:
+    - cell
+    - cell instance
+    - distribcell
+    - energy
+    - material
+    - mesh
+    - universe
+
+Note that there is no difference between the analog, tracklength, and collision
+estimators in random ray mode as individual particles are not being simulated.
+Tracklength-style tally estimation is inherent to the random ray method.
+
+--------
+Plotting
+--------
+
+Visualization of geometry is handled in the same way as normal with OpenMC (see
+:ref:`plotting guide <usersguide_plots>` for more details). That is, ``openmc
+--plot`` is handled without any modifications, as the random ray solver uses the
+same geometry definition as in Monte Carlo.
+
+In addition to OpenMC's standard geometry plotting mode, the random ray solver
+also features an additional method of data visualization. If a ``plots.xml``
+file is present, any voxel plots that are defined will be output at the end of a
+random ray simulation. Rather than being stored in HDF5 file format, the random
+ray plotting will generate ``.vtk`` files that can be directly read and plotted
+with `Paraview <https://www.paraview.org/>`_.
+
+In fixed source Monte Carlo (MC) simulations, by default the only thing global
+tally provided is the leakage fraction. In a k-eigenvalue MC simulation, by
+default global tallies are collected for the eigenvalue and leakage fraction.
+Spatial flux information must be manually requested, and often fine-grained
+spatial meshes are considered costly/unnecessary, so it is impractical in MC
+mode to plot spatial flux or power info by default. Conversely, in random ray,
+the solver functions by estimating the multigroup source and flux spectrums in
+every fine-grained FSR each iteration. Thus, for random ray, in both fixed
+source and eigenvalue simulations, the simulation always finishes with a well
+converged flux estimate for all areas. As such, it is much more common in random
+ray, MOC, and other deterministic codes to provide spatial flux information by
+default. In the future, all FSR data will be made available in the statepoint
+file, which facilitates plotting and manipulation through the Python API; at
+present, statepoint support is not available.
+
+Only voxel plots will be used to generate output; other plot types present in
+the ``plots.xml`` file will be ignored. The following fields will be written to
+the VTK structured grid file:
+
+    - material
+    - FSR index
+    - flux spectrum (for each energy group)
+    - total fission source (integrated across all energy groups)
+
+------------------------------------------
+Inputting Multigroup Cross Sections (MGXS)
+------------------------------------------
+
+Multigroup cross sections for use with OpenMC's random ray solver are input the
+same way as with OpenMC's traditional multigroup Monte Carlo mode. There is more
+information on generating multigroup cross sections via OpenMC in the
+:ref:`multigroup materials <create_mgxs>` user guide. You may also wish to
+use an existing multigroup library. An example of using OpenMC's Python
+interface to generate a correctly formatted ``mgxs.h5`` input file is given
+in the `OpenMC Jupyter notebook collection
+<https://nbviewer.org/github/openmc-dev/openmc-notebooks/blob/main/mg-mode-part-i.ipynb>`_.
+
+.. note::
+    Currently only isotropic and isothermal multigroup cross sections are
+    supported in random ray mode. To represent multiple material temperatures,
+    separate materials can be defined each with a separate multigroup dataset
+    corresponding to a given temperature.
+
+---------------------------------------
+Putting it All Together: Example Inputs
+---------------------------------------
+
+An example of a settings definition for random ray is given below::
+
+    # Geometry and MGXS material definition of 2x2 lattice (not shown)
+    pitch = 1.26
+    group_edges = [1e-5, 0.0635, 10.0, 1.0e2, 1.0e3, 0.5e6, 1.0e6, 20.0e6]
+    ...
+
+    # Instantiate a settings object for a random ray solve
+    settings = openmc.Settings()
+    settings.energy_mode = "multi-group"
+    settings.batches = 1200
+    settings.inactive = 600
+    settings.particles = 2000
+
+    settings.random_ray['distance_inactive'] = 40.0
+    settings.random_ray['distance_active'] = 400.0
+
+    # Create an initial uniform spatial source distribution for sampling rays
+    lower_left  = (-pitch, -pitch, -pitch)
+    upper_right = ( pitch,  pitch,  pitch)
+    uniform_dist = openmc.stats.Box(lower_left, upper_right)
+    settings.random_ray['ray_source'] = openmc.IndependentSource(space=uniform_dist)
+
+    settings.export_to_xml()
+
+    # Define tallies
+
+    # Create a mesh filter
+    mesh = openmc.RegularMesh()
+    mesh.dimension = (2, 2)
+    mesh.lower_left = (-pitch/2, -pitch/2)
+    mesh.upper_right = (pitch/2, pitch/2)
+    mesh_filter = openmc.MeshFilter(mesh)
+
+    # Create a multigroup energy filter
+    energy_filter = openmc.EnergyFilter(group_edges)
+
+    # Create tally using our two filters and add scores
+    tally = openmc.Tally()
+    tally.filters = [mesh_filter, energy_filter]
+    tally.scores = ['flux', 'fission', 'nu-fission']
+
+    # Instantiate a Tallies collection and export to XML
+    tallies = openmc.Tallies([tally])
+    tallies.export_to_xml()
+
+    # Create voxel plot
+    plot = openmc.Plot()
+    plot.origin = [0, 0, 0]
+    plot.width = [2*pitch, 2*pitch, 1]
+    plot.pixels = [1000, 1000, 1]
+    plot.type = 'voxel'
+
+    # Instantiate a Plots collection and export to XML
+    plots = openmc.Plots([plot])
+    plots.export_to_xml()
+
+All other inputs (e.g., geometry, materials) will be unchanged from a typical
+Monte Carlo run (see the :ref:`geometry <usersguide_geometry>` and
+:ref:`multigroup materials <create_mgxs>` user guides for more information).
+
+There is also a complete example of a pincell available in the
+``openmc/examples/pincell_random_ray`` folder.

docs/source/usersguide/settings.rst

@@ -420,6 +420,54 @@ string, which gets passed down to the ``openmc_create_source()`` function::
 
   settings.source = openmc.CompiledSource('libsource.so', '3.5e6')
 
+.. _usersguide_source_constraints:
+
+Source Constraints
+------------------
+
+All source classes in OpenMC have the ability to apply a set of "constraints"
+that limit which sampled source sites are actually used for transport. The most
+common use case is to sample source sites over some simple spatial distribution
+(e.g., uniform over a box) and then only accept those that appear in a given
+cell or material. This can be done with a domain constraint, which can be
+specified as follows::
+
+  source_cell = openmc.Cell(...)
+  ...
+
+  spatial_dist = openmc.stats.Box((-10., -10., -10.), (10., 10., 10.))
+  source = openmc.IndependentSource(
+      space=spatial_dist,
+      constraints={'domains': [source_cell]}
+  )
+
+For k-eigenvalue problems, a convenient constraint is available that limits
+source sites to those sampled in a fissionable material::
+
+  source = openmc.IndependentSource(
+      space=spatial_dist, constraints={'fissionable': True}
+  )
+
+Constraints can also be placed on a range of energies or times::
+
+  # Only use source sites between 500 keV and 1 MeV and with times under 1 sec
+  source = openmc.FileSource(
+      'source.h5',
+      constraints={'energy_bounds': [500.0e3, 1.0e6], 'time_bounds': [0.0, 1.0]}
+  )
+
+Normally, when a source site is rejected, a new one will be resampled until one
+is found that meets the constraints. However, the rejection strategy can be
+changed so that a rejected site will just not be simulated by specifying::
+
+  source = openmc.IndependentSource(
+      space=spatial_dist,
+      constraints={'domains': [cell], 'rejection_strategy': 'kill'}
+  )
+
+In this case, the actual number of particles simulated may be less than what you
+specified in :attr:`Settings.particles`.
+
 .. _usersguide_entropy:
 
 ---------------
@@ -628,3 +676,37 @@ instance, whereas the :meth:`openmc.Track.filter` method returns a new
           .. code-block:: sh
 
             openmc-track-combine tracks_p*.h5 --out tracks.h5
+
+-----------------------
+Restarting a Simulation
+-----------------------
+
+OpenMC can be run in a mode where it reads in a statepoint file and continues a
+simulation from the ending point of the statepoint file. A restart simulation
+can be performed by passing the path to the statepoint file to the OpenMC
+executable:
+
+.. code-block:: sh
+
+    openmc -r statepoint.100.h5
+
+From the Python API, the `restart_file` argument provides the same behavior:
+
+.. code-block:: python
+
+    openmc.run(restart_file='statepoint.100.h5')
+
+or if using the :class:`~openmc.Model` class:
+
+.. code-block:: python
+
+    model.run(restart_file='statepoint.100.h5')
+
+The restart simulation will execute until the number of batches specified in the
+:class:`~openmc.Settings` object on a model (or in the :ref:`settings XML file
+<io_settings>`) is satisfied. Note that if the number of batches in the
+statepoint file is the same as that specified in the settings object (i.e., if
+the inputs were not modified before the restart run), no particles will be
+transported and OpenMC will exit immediately.
+
+.. note:: A statepoint file must match the input model to be successfully used in a restart simulation.

examples/assembly/assembly.py

@@ -112,11 +112,10 @@ def assembly_model():
     model.settings.batches = 150
     model.settings.inactive = 50
     model.settings.particles = 1000
-    model.settings.source = openmc.IndependentSource(space=openmc.stats.Box(
-        (-pitch/2, -pitch/2, -1),
-        (pitch/2, pitch/2, 1),
-        only_fissionable=True
-    ))
+    model.settings.source = openmc.IndependentSource(
+        space=openmc.stats.Box((-pitch/2, -pitch/2, -1), (pitch/2, pitch/2, 1)),
+        constraints={'fissionable': True}
+    )
 
     # NOTE: We never actually created a Materials object. When you export/run
     # using the Model object, if no materials were assigned it will look through

examples/custom_source/CMakeLists.txt

@@ -1,4 +1,4 @@
-cmake_minimum_required(VERSION 3.3 FATAL_ERROR)
+cmake_minimum_required(VERSION 3.10 FATAL_ERROR)
 project(openmc_sources CXX)
 add_library(source SHARED source_ring.cpp)
 find_package(OpenMC REQUIRED)

examples/parameterized_custom_source/CMakeLists.txt

@@ -1,4 +1,4 @@
-cmake_minimum_required(VERSION 3.3 FATAL_ERROR)
+cmake_minimum_required(VERSION 3.10 FATAL_ERROR)
 project(openmc_sources CXX)
 add_library(parameterized_source SHARED parameterized_source_ring.cpp)
 find_package(OpenMC REQUIRED)

examples/pincell_depletion/restart_depletion.py

@@ -26,8 +26,9 @@ settings.particles = 10000
 
 # Create an initial uniform spatial source distribution over fissionable zones
 bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1]
-uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
-settings.source = openmc.IndependentSource(space=uniform_dist)
+uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:])
+settings.source = openmc.IndependentSource(
+    space=uniform_dist, constraints={'fissionable': True})
 
 entropy_mesh = openmc.RegularMesh()
 entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50]

examples/pincell_depletion/run_depletion.py

@@ -71,8 +71,9 @@ settings.particles = 1000
 
 # Create an initial uniform spatial source distribution over fissionable zones
 bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1]
-uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
-settings.source = openmc.IndependentSource(space=uniform_dist)
+uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:])
+settings.source = openmc.IndependentSource(
+    space=uniform_dist, constraints={'fissionable': True})
 
 entropy_mesh = openmc.RegularMesh()
 entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50]

examples/pincell_random_ray/build_xml.py

@@ -0,0 +1,203 @@
+import numpy as np
+import openmc
+import openmc.mgxs
+
+###############################################################################
+# Create multigroup data
+
+# Instantiate the energy group data
+group_edges = [1e-5, 0.0635, 10.0, 1.0e2, 1.0e3, 0.5e6, 1.0e6, 20.0e6]
+groups = openmc.mgxs.EnergyGroups(group_edges)
+
+# Instantiate the 7-group (C5G7) cross section data
+uo2_xsdata = openmc.XSdata('UO2', groups)
+uo2_xsdata.order = 0
+uo2_xsdata.set_total(
+    [0.1779492, 0.3298048, 0.4803882, 0.5543674, 0.3118013, 0.3951678,
+     0.5644058])
+uo2_xsdata.set_absorption([8.0248e-03, 3.7174e-03, 2.6769e-02, 9.6236e-02,
+                           3.0020e-02, 1.1126e-01, 2.8278e-01])
+scatter_matrix = np.array(
+    [[[0.1275370, 0.0423780, 0.0000094, 0.0000000, 0.0000000, 0.0000000, 0.0000000],
+      [0.0000000, 0.3244560, 0.0016314, 0.0000000, 0.0000000, 0.0000000, 0.0000000],
+      [0.0000000, 0.0000000, 0.4509400, 0.0026792, 0.0000000, 0.0000000, 0.0000000],
+      [0.0000000, 0.0000000, 0.0000000, 0.4525650, 0.0055664, 0.0000000, 0.0000000],
+      [0.0000000, 0.0000000, 0.0000000, 0.0001253, 0.2714010, 0.0102550, 0.0000000],
+      [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0012968, 0.2658020, 0.0168090],
+      [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0085458, 0.2730800]]])
+scatter_matrix = np.rollaxis(scatter_matrix, 0, 3)
+uo2_xsdata.set_scatter_matrix(scatter_matrix)
+uo2_xsdata.set_fission([7.21206e-03, 8.19301e-04, 6.45320e-03,
+                        1.85648e-02, 1.78084e-02, 8.30348e-02,
+                        2.16004e-01])
+uo2_xsdata.set_nu_fission([2.005998e-02, 2.027303e-03, 1.570599e-02,
+                           4.518301e-02, 4.334208e-02, 2.020901e-01,
+                           5.257105e-01])
+uo2_xsdata.set_chi([5.8791e-01, 4.1176e-01, 3.3906e-04, 1.1761e-07, 0.0000e+00,
+                    0.0000e+00, 0.0000e+00])
+
+h2o_xsdata = openmc.XSdata('LWTR', groups)
+h2o_xsdata.order = 0
+h2o_xsdata.set_total([0.15920605, 0.412969593, 0.59030986, 0.58435,
+                      0.718, 1.2544497, 2.650379])
+h2o_xsdata.set_absorption([6.0105e-04, 1.5793e-05, 3.3716e-04,
+                           1.9406e-03, 5.7416e-03, 1.5001e-02,
+                           3.7239e-02])
+scatter_matrix = np.array(
+    [[[0.0444777, 0.1134000, 0.0007235, 0.0000037, 0.0000001, 0.0000000, 0.0000000],
+      [0.0000000, 0.2823340, 0.1299400, 0.0006234, 0.0000480, 0.0000074, 0.0000010],
+      [0.0000000, 0.0000000, 0.3452560, 0.2245700, 0.0169990, 0.0026443, 0.0005034],
+      [0.0000000, 0.0000000, 0.0000000, 0.0910284, 0.4155100, 0.0637320, 0.0121390],
+      [0.0000000, 0.0000000, 0.0000000, 0.0000714, 0.1391380, 0.5118200, 0.0612290],
+      [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0022157, 0.6999130, 0.5373200],
+      [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.1324400, 2.4807000]]])
+scatter_matrix = np.rollaxis(scatter_matrix, 0, 3)
+h2o_xsdata.set_scatter_matrix(scatter_matrix)
+
+mg_cross_sections_file = openmc.MGXSLibrary(groups)
+mg_cross_sections_file.add_xsdatas([uo2_xsdata, h2o_xsdata])
+mg_cross_sections_file.export_to_hdf5()
+
+###############################################################################
+# Create materials for the problem
+
+# Instantiate some Materials and register the appropriate macroscopic data
+uo2 = openmc.Material(name='UO2 fuel')
+uo2.set_density('macro', 1.0)
+uo2.add_macroscopic('UO2')
+
+water = openmc.Material(name='Water')
+water.set_density('macro', 1.0)
+water.add_macroscopic('LWTR')
+
+# Instantiate a Materials collection and export to XML
+materials_file = openmc.Materials([uo2, water])
+materials_file.cross_sections = "mgxs.h5"
+materials_file.export_to_xml()
+
+###############################################################################
+# Define problem geometry
+
+# The geometry we will define a simplified pincell with fuel radius 0.54 cm
+# surrounded by moderator (same as in the multigroup example).
+# In random ray, we typically want several radial regions and azimuthal
+# sectors in both the fuel and moderator areas of the pincell. This is
+# due to the flat source approximation requiring that source regions are
+# small compared to the typical mean free path of a neutron. Below we
+# sudivide the basic pincell into 8 aziumthal sectors (pizza slices) and
+# 5 concentric rings in both the fuel and moderator.
+
+# TODO: When available in OpenMC, use cylindrical lattice instead to
+# simplify definition and improve runtime performance.
+
+pincell_base = openmc.Universe()
+
+# These are the subdivided radii (creating 5 concentric regions in the
+# fuel and moderator)
+ring_radii = [0.241, 0.341, 0.418, 0.482, 0.54, 0.572, 0.612, 0.694, 0.786]
+fills = [uo2, uo2, uo2, uo2, uo2, water, water, water, water, water]
+
+# We then create cells representing the bounded rings, with special
+# treatment for both the innermost and outermost cells
+cells = []
+for r in range(10):
+    cell = []
+    if r == 0:
+        outer_bound = openmc.ZCylinder(r=ring_radii[r])
+        cell = openmc.Cell(fill=fills[r], region=-outer_bound)
+    elif r == 9:
+        inner_bound = openmc.ZCylinder(r=ring_radii[r-1])
+        cell = openmc.Cell(fill=fills[r], region=+inner_bound)
+    else:
+        inner_bound = openmc.ZCylinder(r=ring_radii[r-1])
+        outer_bound = openmc.ZCylinder(r=ring_radii[r])
+        cell = openmc.Cell(fill=fills[r], region=+inner_bound & -outer_bound)
+    pincell_base.add_cell(cell)
+
+# We then generate 8 planes to bound 8 azimuthal sectors
+azimuthal_planes = []
+for i in range(8):
+    angle = 2 * i * openmc.pi / 8
+    normal_vector = (-openmc.sin(angle), openmc.cos(angle), 0)
+    azimuthal_planes.append(openmc.Plane(a=normal_vector[0], b=normal_vector[1], c=normal_vector[2], d=0))
+
+# Create a cell for each azimuthal sector using the pincell base class
+azimuthal_cells = []
+for i in range(8):
+    azimuthal_cell = openmc.Cell(name=f'azimuthal_cell_{i}')
+    azimuthal_cell.fill = pincell_base
+    azimuthal_cell.region = +azimuthal_planes[i] & -azimuthal_planes[(i+1) % 8]
+    azimuthal_cells.append(azimuthal_cell)
+
+# Create the (subdivided) geometry with the azimuthal universes
+pincell = openmc.Universe(cells=azimuthal_cells)
+
+# Create a region represented as the inside of a rectangular prism
+pitch = 1.26
+box = openmc.model.RectangularPrism(pitch, pitch, boundary_type='reflective')
+pincell_bounded = openmc.Cell(fill=pincell, region=-box, name='pincell')
+
+# Create a geometry (specifying merge surfaces option to remove
+# all the redundant cylinder/plane surfaces) and export to XML
+geometry = openmc.Geometry([pincell_bounded], merge_surfaces=True)
+geometry.export_to_xml()
+
+###############################################################################
+# Define problem settings
+
+# Instantiate a Settings object, set all runtime parameters, and export to XML
+settings = openmc.Settings()
+settings.energy_mode = "multi-group"
+settings.batches = 600
+settings.inactive = 300
+settings.particles = 50
+
+# Create an initial uniform spatial source distribution for sampling rays.
+# Note that this must be uniform in space and angle.
+lower_left = (-pitch/2, -pitch/2, -1)
+upper_right = (pitch/2, pitch/2, 1)
+uniform_dist = openmc.stats.Box(lower_left, upper_right)
+settings.random_ray['ray_source'] = openmc.IndependentSource(space=uniform_dist)
+settings.random_ray['distance_inactive'] = 40.0
+settings.random_ray['distance_active'] = 400.0
+
+settings.export_to_xml()
+
+###############################################################################
+# Define tallies
+
+# Create a mesh that will be used for tallying
+mesh = openmc.RegularMesh()
+mesh.dimension = (2, 2)
+mesh.lower_left = (-pitch/2, -pitch/2)
+mesh.upper_right = (pitch/2, pitch/2)
+
+# Create a mesh filter that can be used in a tally
+mesh_filter = openmc.MeshFilter(mesh)
+
+# Let's also create a filter to measure each group
+# indepdendently
+energy_filter = openmc.EnergyFilter(group_edges)
+
+# Now use the mesh filter in a tally and indicate what scores are desired
+tally = openmc.Tally(name="Mesh and Energy tally")
+tally.filters = [mesh_filter, energy_filter]
+tally.scores = ['flux', 'fission', 'nu-fission']
+
+# Instantiate a Tallies collection and export to XML
+tallies = openmc.Tallies([tally])
+tallies.export_to_xml()
+
+###############################################################################
+#                   Exporting to OpenMC plots.xml file
+###############################################################################
+
+plot = openmc.Plot()
+plot.origin = [0, 0, 0]
+plot.width = [pitch, pitch, pitch]
+plot.pixels = [1000, 1000, 1]
+plot.type = 'voxel'
+
+# Instantiate a Plots collection and export to XML
+plots = openmc.Plots([plot])
+plots.export_to_xml()

include/openmc/boundary_condition.h

@@ -10,6 +10,7 @@ namespace openmc {
 
 // Forward declare some types used in function arguments.
 class Particle;
+class RandomRay;
 class Surface;
 
 //==============================================================================

include/openmc/capi.h

@@ -150,6 +150,7 @@ int openmc_sphharm_filter_get_cosine(int32_t index, char cosine[]);
 int openmc_sphharm_filter_set_order(int32_t index, int order);
 int openmc_sphharm_filter_set_cosine(int32_t index, const char cosine[]);
 int openmc_statepoint_write(const char* filename, bool* write_source);
+int openmc_statepoint_load(const char* filename);
 int openmc_tally_allocate(int32_t index, const char* type);
 int openmc_tally_get_active(int32_t index, bool* active);
 int openmc_tally_get_estimator(int32_t index, int* estimator);

include/openmc/cell.h

@@ -249,6 +249,42 @@ public:
   std::unordered_map<int32_t, vector<int32_t>> get_contained_cells(
     int32_t instance = 0, Position* hint = nullptr) const;
 
+  //! Determine the material index corresponding to a specific cell instance,
+  //! taking into account presence of distribcell material
+  //! \param[in] instance of the cell
+  //! \return material index
+  int32_t material(int32_t instance) const
+  {
+    // If distributed materials are used, then each instance has its own
+    // material definition. If distributed materials are not used, then
+    // all instances used the same material stored at material_[0]. The
+    // presence of distributed materials is inferred from the size of
+    // the material_ vector being greater than one.
+    if (material_.size() > 1) {
+      return material_[instance];
+    } else {
+      return material_[0];
+    }
+  }
+
+  //! Determine the temperature index corresponding to a specific cell instance,
+  //! taking into account presence of distribcell temperature
+  //! \param[in] instance of the cell
+  //! \return temperature index
+  double sqrtkT(int32_t instance) const
+  {
+    // If distributed materials are used, then each instance has its own
+    // temperature definition. If distributed materials are not used, then
+    // all instances used the same temperature stored at sqrtkT_[0]. The
+    // presence of distributed materials is inferred from the size of
+    // the sqrtkT_ vector being greater than one.
+    if (sqrtkT_.size() > 1) {
+      return sqrtkT_[instance];
+    } else {
+      return sqrtkT_[0];
+    }
+  }
+
 protected:
   //! Determine the path to this cell instance in the geometry hierarchy
   //! \param[in] instance of the cell to find parent cells for

include/openmc/constants.h

@@ -340,6 +340,8 @@ enum class RunMode {
   VOLUME
 };
 
+enum class SolverType { MONTE_CARLO, RANDOM_RAY };
+
 //==============================================================================
 // Geometry Constants
 

include/openmc/dagmc.h

@@ -3,7 +3,8 @@
 
 namespace openmc {
 extern "C" const bool DAGMC_ENABLED;
-}
+extern "C" const bool UWUW_ENABLED;
+} // namespace openmc
 
 // always include the XML interface header
 #include "openmc/xml_interface.h"
@@ -120,6 +121,12 @@ public:
   void write_uwuw_materials_xml(
     const std::string& outfile = "uwuw_materials.xml") const;
 
+  //! Assign a material to a cell from uwuw material library
+  //! \param[in] vol_handle The DAGMC material assignment string
+  //! \param[in] c The OpenMC cell to which the material is assigned
+  void uwuw_assign_material(
+    moab::EntityHandle vol_handle, std::unique_ptr<DAGCell>& c) const;
+
   //! Assign a material to a cell based
   //! \param[in] mat_string The DAGMC material assignment string
   //! \param[in] c The OpenMC cell to which the material is assigned

include/openmc/distribution_spatial.h

@@ -116,6 +116,11 @@ public:
   //! \return Sampled element index and position within that element
   std::pair<int32_t, Position> sample_mesh(uint64_t* seed) const;
 
+  //! Sample a mesh element
+  //! \param seed Pseudorandom number seed pointer
+  //! \return Sampled element index
+  int32_t sample_element_index(uint64_t* seed) const;
+
   //! For unstructured meshes, ensure that elements are all linear tetrahedra
   void check_element_types() const;
 

include/openmc/mesh.h

@@ -659,11 +659,6 @@ protected:
   //! Set the length multiplier to apply to each point in the mesh
   void set_length_multiplier(const double length_multiplier);
 
-  // Data members
-  double length_multiplier_ {
-    1.0}; //!< Constant multiplication factor to apply to mesh coordinates
-  bool specified_length_multiplier_ {false};
-
   //! Sample barycentric coordinates given a seed and the vertex positions and
   //! return the sampled position
   //
@@ -672,6 +667,11 @@ protected:
   //! \return Sampled position within the tetrahedron
   Position sample_tet(std::array<Position, 4> coords, uint64_t* seed) const;
 
+  // Data members
+  double length_multiplier_ {
+    -1.0};              //!< Multiplicative factor applied to mesh coordinates
+  std::string options_; //!< Options for search data structures
+
 private:
   //! Setup method for the mesh. Builds data structures,
   //! sets up element mapping, creates bounding boxes, etc.

include/openmc/mgxs.h

@@ -29,7 +29,6 @@ private:
   int num_delayed_groups; // number of delayed neutron groups
   vector<XsData> xs;      // Cross section data
   // MGXS Incoming Flux Angular grid information
-  bool is_isotropic; // used to skip search for angle indices if isotropic
   int n_pol;
   int n_azi;
   vector<double> polar;
@@ -85,6 +84,8 @@ public:
   std::string name; // name of dataset, e.g., UO2
   double awr;       // atomic weight ratio
   bool fissionable; // Is this fissionable
+  bool is_isotropic {
+    true}; // used to skip search for angle indices if isotropic
 
   Mgxs() = default;
 

include/openmc/output.h

@@ -57,6 +57,8 @@ void print_results();
 
 void write_tallies();
 
+void show_time(const char* label, double secs, int indent_level = 0);
+
 } // namespace openmc
 #endif // OPENMC_OUTPUT_H
 
@@ -80,4 +82,4 @@ struct formatter<std::array<T, 2>> {
   }
 };
 
-} // namespace fmt
+} // namespace fmt

include/openmc/particle_data.h

@@ -27,9 +27,6 @@ constexpr int MAX_DELAYED_GROUPS {8};
 
 constexpr double CACHE_INVALID {-1.0};
 
-// Maximum number of collisions/crossings
-constexpr int MAX_EVENTS {1000000};
-
 //==========================================================================
 // Aliases and type definitions
 
@@ -262,6 +259,10 @@ public:
   int& cell_last(int i) { return cell_last_[i]; }
   const int& cell_last(int i) const { return cell_last_[i]; }
 
+  // Coordinates at birth
+  Position& r_born() { return r_born_; }
+  const Position& r_born() const { return r_born_; }
+
   // Coordinates of last collision or reflective/periodic surface
   // crossing for current tallies
   Position& r_last_current() { return r_last_current_; }
@@ -323,6 +324,7 @@ private:
   int n_coord_last_ {1};  //!< number of current coordinates
   vector<int> cell_last_; //!< coordinates for all levels
 
+  Position r_born_;         //!< coordinates at birth
   Position r_last_current_; //!< coordinates of the last collision or
                             //!< reflective/periodic surface crossing for
                             //!< current tallies

include/openmc/plot.h

@@ -100,6 +100,7 @@ public:
   virtual void print_info() const = 0;
 
   const std::string& path_plot() const { return path_plot_; }
+  std::string& path_plot() { return path_plot_; }
   int id() const { return id_; }
   int level() const { return level_; }
 

include/openmc/random_dist.h

@@ -48,9 +48,9 @@ extern "C" double maxwell_spectrum(double T, uint64_t* seed);
 extern "C" double watt_spectrum(double a, double b, uint64_t* seed);
 
 //==============================================================================
-//! Samples an energy from the Gaussian energy-dependent fission distribution.
+//! Samples an energy from the Gaussian distribution.
 //!
-//! Samples from a Normal distribution with a given mean and standard deviation
+//! Samples from a normal distribution with a given mean and standard deviation
 //! The PDF is defined as s(x) = (1/2*sigma*sqrt(2) * e-((mu-x)/2*sigma)^2
 //! Its sampled according to
 //! http://www-pdg.lbl.gov/2009/reviews/rpp2009-rev-monte-carlo-techniques.pdf

include/openmc/random_ray/flat_source_domain.h

@@ -0,0 +1,141 @@
+#ifndef OPENMC_RANDOM_RAY_FLAT_SOURCE_DOMAIN_H
+#define OPENMC_RANDOM_RAY_FLAT_SOURCE_DOMAIN_H
+
+#include "openmc/openmp_interface.h"
+#include "openmc/position.h"
+
+namespace openmc {
+
+/*
+ * The FlatSourceDomain class encompasses data and methods for storing
+ * scalar flux and source region for all flat source regions in a
+ * random ray simulation domain.
+ */
+
+class FlatSourceDomain {
+public:
+  //----------------------------------------------------------------------------
+  // Helper Structs
+
+  // A mapping object that is used to map between a specific random ray
+  // source region and an OpenMC native tally bin that it should score to
+  // every iteration.
+  struct TallyTask {
+    int tally_idx;
+    int filter_idx;
+    int score_idx;
+    int score_type;
+    TallyTask(int tally_idx, int filter_idx, int score_idx, int score_type)
+      : tally_idx(tally_idx), filter_idx(filter_idx), score_idx(score_idx),
+        score_type(score_type)
+    {}
+  };
+
+  //----------------------------------------------------------------------------
+  // Constructors
+  FlatSourceDomain();
+
+  //----------------------------------------------------------------------------
+  // Methods
+  void update_neutron_source(double k_eff);
+  double compute_k_eff(double k_eff_old) const;
+  void normalize_scalar_flux_and_volumes(
+    double total_active_distance_per_iteration);
+  int64_t add_source_to_scalar_flux();
+  void batch_reset();
+  void convert_source_regions_to_tallies();
+  void random_ray_tally() const;
+  void accumulate_iteration_flux();
+  void output_to_vtk() const;
+  void all_reduce_replicated_source_regions();
+
+  //----------------------------------------------------------------------------
+  // Public Data members
+
+  bool mapped_all_tallies_ {false}; // If all source regions have been visited
+
+  int64_t n_source_regions_ {0}; // Total number of source regions in the model
+
+  // 1D array representing source region starting offset for each OpenMC Cell
+  // in model::cells
+  vector<int64_t> source_region_offsets_;
+
+  // 1D arrays representing values for all source regions
+  vector<OpenMPMutex> lock_;
+  vector<int> was_hit_;
+  vector<double> volume_;
+  vector<int> position_recorded_;
+  vector<Position> position_;
+
+  // 2D arrays stored in 1D representing values for all source regions x energy
+  // groups
+  vector<float> scalar_flux_old_;
+  vector<float> scalar_flux_new_;
+  vector<float> source_;
+
+  //----------------------------------------------------------------------------
+  // Private data members
+private:
+  int negroups_;                  // Number of energy groups in simulation
+  int64_t n_source_elements_ {0}; // Total number of source regions in the model
+                                  // times the number of energy groups
+
+  // 2D array representing values for all source regions x energy groups x tally
+  // tasks
+  vector<vector<TallyTask>> tally_task_;
+
+  // 1D arrays representing values for all source regions
+  vector<int> material_;
+  vector<double> volume_t_;
+
+  // 2D arrays stored in 1D representing values for all source regions x energy
+  // groups
+  vector<float> scalar_flux_final_;
+
+}; // class FlatSourceDomain
+
+//============================================================================
+//! Non-member functions
+//============================================================================
+
+// Returns the inputted value in big endian byte ordering. If the system is
+// little endian, the byte ordering is flipped. If the system is big endian,
+// the inputted value is returned as is. This function is necessary as
+// .vtk binary files use big endian byte ordering.
+template<typename T>
+T convert_to_big_endian(T in)
+{
+  // 4 byte integer
+  uint32_t test = 1;
+
+  // 1 byte pointer to first byte of test integer
+  uint8_t* ptr = reinterpret_cast<uint8_t*>(&test);
+
+  // If the first byte of test is 0, then the system is big endian. In this
+  // case, we don't have to do anything as .vtk files are big endian
+  if (*ptr == 0)
+    return in;
+
+  // Otherwise, the system is in little endian, so we need to flip the
+  // endianness
+  uint8_t* orig = reinterpret_cast<uint8_t*>(&in);
+  uint8_t swapper[sizeof(T)];
+  for (int i = 0; i < sizeof(T); i++) {
+    swapper[i] = orig[sizeof(T) - i - 1];
+  }
+  T out = *reinterpret_cast<T*>(&swapper);
+  return out;
+}
+
+template<typename T>
+void parallel_fill(vector<T>& arr, T value)
+{
+#pragma omp parallel for schedule(static)
+  for (int i = 0; i < arr.size(); i++) {
+    arr[i] = value;
+  }
+}
+
+} // namespace openmc
+
+#endif // OPENMC_RANDOM_RAY_FLAT_SOURCE_DOMAIN_H

include/openmc/random_ray/random_ray.h

@@ -0,0 +1,55 @@
+#ifndef OPENMC_RANDOM_RAY_H
+#define OPENMC_RANDOM_RAY_H
+
+#include "openmc/memory.h"
+#include "openmc/particle.h"
+#include "openmc/random_ray/flat_source_domain.h"
+#include "openmc/source.h"
+
+namespace openmc {
+
+/*
+ * The RandomRay class encompasses data and methods for transporting random rays
+ * through the model. It is a small extension of the Particle class.
+ */
+
+// TODO: Inherit from GeometryState instead of Particle
+class RandomRay : public Particle {
+public:
+  //----------------------------------------------------------------------------
+  // Constructors
+  RandomRay();
+  RandomRay(uint64_t ray_id, FlatSourceDomain* domain);
+
+  //----------------------------------------------------------------------------
+  // Methods
+  void event_advance_ray();
+  void attenuate_flux(double distance, bool is_active);
+  void initialize_ray(uint64_t ray_id, FlatSourceDomain* domain);
+  uint64_t transport_history_based_single_ray();
+
+  //----------------------------------------------------------------------------
+  // Static data members
+  static double distance_inactive_;      // Inactive (dead zone) ray length
+  static double distance_active_;        // Active ray length
+  static unique_ptr<Source> ray_source_; // Starting source for ray sampling
+
+  //----------------------------------------------------------------------------
+  // Public data members
+  vector<float> angular_flux_;
+
+private:
+  //----------------------------------------------------------------------------
+  // Private data members
+  vector<float> delta_psi_;
+  int negroups_;
+  FlatSourceDomain* domain_ {nullptr}; // pointer to domain that has flat source
+                                       // data needed for ray transport
+  double distance_travelled_ {0};
+  bool is_active_ {false};
+  bool is_alive_ {true};
+}; // class RandomRay
+
+} // namespace openmc
+
+#endif // OPENMC_RANDOM_RAY_H

include/openmc/random_ray/random_ray_simulation.h

@@ -0,0 +1,59 @@
+#ifndef OPENMC_RANDOM_RAY_SIMULATION_H
+#define OPENMC_RANDOM_RAY_SIMULATION_H
+
+#include "openmc/random_ray/flat_source_domain.h"
+
+namespace openmc {
+
+/*
+ * The RandomRaySimulation class encompasses data and methods for running a
+ * random ray simulation.
+ */
+
+class RandomRaySimulation {
+public:
+  //----------------------------------------------------------------------------
+  // Constructors
+  RandomRaySimulation();
+
+  //----------------------------------------------------------------------------
+  // Methods
+  void simulate();
+  void reduce_simulation_statistics();
+  void output_simulation_results() const;
+  void instability_check(
+    int64_t n_hits, double k_eff, double& avg_miss_rate) const;
+  void print_results_random_ray(uint64_t total_geometric_intersections,
+    double avg_miss_rate, int negroups, int64_t n_source_regions) const;
+
+  //----------------------------------------------------------------------------
+  // Data members
+private:
+  // Contains all flat source region data
+  FlatSourceDomain domain_;
+
+  // Random ray eigenvalue
+  double k_eff_ {1.0};
+
+  // Tracks the average FSR miss rate for analysis and reporting
+  double avg_miss_rate_ {0.0};
+
+  // Tracks the total number of geometric intersections by all rays for
+  // reporting
+  uint64_t total_geometric_intersections_ {0};
+
+  // Number of energy groups
+  int negroups_;
+
+}; // class RandomRaySimulation
+
+//============================================================================
+//! Non-member functions
+//============================================================================
+
+void openmc_run_random_ray();
+void validate_random_ray_inputs();
+
+} // namespace openmc
+
+#endif // OPENMC_RANDOM_RAY_SIMULATION_H

include/openmc/settings.h

@@ -93,8 +93,8 @@ extern "C" int32_t gen_per_batch; //!< number of generations per batch
 extern "C" int64_t n_particles;   //!< number of particles per generation
 
 extern int64_t
-  max_particles_in_flight; //!< Max num. event-based particles in flight
-
+  max_particles_in_flight;      //!< Max num. event-based particles in flight
+extern int max_particle_events; //!< Maximum number of particle events
 extern ElectronTreatment
   electron_treatment; //!< how to treat secondary electrons
 extern array<double, 4>
@@ -112,8 +112,9 @@ extern ResScatMethod res_scat_method; //!< resonance upscattering method
 extern double res_scat_energy_min; //!< Min energy in [eV] for res. upscattering
 extern double res_scat_energy_max; //!< Max energy in [eV] for res. upscattering
 extern vector<std::string>
-  res_scat_nuclides;     //!< Nuclides using res. upscattering treatment
-extern RunMode run_mode; //!< Run mode (eigenvalue, fixed src, etc.)
+  res_scat_nuclides;           //!< Nuclides using res. upscattering treatment
+extern RunMode run_mode;       //!< Run mode (eigenvalue, fixed src, etc.)
+extern SolverType solver_type; //!< Solver Type (Monte Carlo or Random Ray)
 extern std::unordered_set<int>
   sourcepoint_batch; //!< Batches when source should be written
 extern std::unordered_set<int>
@@ -139,6 +140,7 @@ extern int trigger_batch_interval; //!< Batch interval for triggers
 extern "C" int verbosity;          //!< How verbose to make output
 extern double weight_cutoff;       //!< Weight cutoff for Russian roulette
 extern double weight_survive;      //!< Survival weight after Russian roulette
+
 } // namespace settings
 
 //==============================================================================

include/openmc/source.h

@@ -4,6 +4,7 @@
 #ifndef OPENMC_SOURCE_H
 #define OPENMC_SOURCE_H
 
+#include <limits>
 #include <unordered_set>
 
 #include "pugixml.hpp"
@@ -39,19 +40,72 @@ extern vector<unique_ptr<Source>> external_sources;
 
 //==============================================================================
 //! Abstract source interface
+//
+//! The Source class provides the interface that must be implemented by derived
+//! classes, namely the sample() method that returns a sampled source site. From
+//! this base class, source rejection is handled within the
+//! sample_with_constraints() method. However, note that some classes directly
+//! check for constraints for efficiency reasons (like IndependentSource), in
+//! which case the constraints_applied() method indicates that constraints
+//! should not be checked a second time from the base class.
 //==============================================================================
 
 class Source {
 public:
+  // Constructors, destructors
+  Source() = default;
+  explicit Source(pugi::xml_node node);
   virtual ~Source() = default;
 
-  // Methods that must be implemented
+  // Methods that can be overridden
+  virtual double strength() const { return strength_; }
+
+  //! Sample a source site and apply constraints
+  //
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Sampled site
+  SourceSite sample_with_constraints(uint64_t* seed) const;
+
+  //! Sample a source site (without applying constraints)
+  //
+  //! Sample from the external source distribution
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Sampled site
   virtual SourceSite sample(uint64_t* seed) const = 0;
 
-  // Methods that can be overridden
-  virtual double strength() const { return 1.0; }
-
   static unique_ptr<Source> create(pugi::xml_node node);
+
+protected:
+  // Domain types
+  enum class DomainType { UNIVERSE, MATERIAL, CELL };
+
+  // Strategy used for rejecting sites when constraints are applied. KILL means
+  // that sites are always accepted but if they don't satisfy constraints, they
+  // are given weight 0. RESAMPLE means that a new source site will be sampled
+  // until constraints are met.
+  enum class RejectionStrategy { KILL, RESAMPLE };
+
+  // Indicates whether derived class already handles constraints
+  virtual bool constraints_applied() const { return false; }
+
+  // Methods for constraints
+  void read_constraints(pugi::xml_node node);
+  bool satisfies_spatial_constraints(Position r) const;
+  bool satisfies_energy_constraints(double E) const;
+  bool satisfies_time_constraints(double time) const;
+
+  // Data members
+  double strength_ {1.0};                  //!< Source strength
+  std::unordered_set<int32_t> domain_ids_; //!< Domains to reject from
+  DomainType domain_type_;                 //!< Domain type for rejection
+  std::pair<double, double> time_bounds_ {-std::numeric_limits<double>::max(),
+    std::numeric_limits<double>::max()}; //!< time limits
+  std::pair<double, double> energy_bounds_ {
+    0, std::numeric_limits<double>::max()}; //!< energy limits
+  bool only_fissionable_ {
+    false}; //!< Whether site must be in fissionable material
+  RejectionStrategy rejection_strategy_ {
+    RejectionStrategy::RESAMPLE}; //!< Procedure for rejecting
 };
 
 //==============================================================================
@@ -73,7 +127,6 @@ public:
 
   // Properties
   ParticleType particle_type() const { return particle_; }
-  double strength() const override { return strength_; }
 
   // Make observing pointers available
   SpatialDistribution* space() const { return space_.get(); }
@@ -81,19 +134,17 @@ public:
   Distribution* energy() const { return energy_.get(); }
   Distribution* time() const { return time_.get(); }
 
-private:
-  // Domain types
-  enum class DomainType { UNIVERSE, MATERIAL, CELL };
+protected:
+  // Indicates whether derived class already handles constraints
+  bool constraints_applied() const override { return true; }
 
+private:
   // Data members
   ParticleType particle_ {ParticleType::neutron}; //!< Type of particle emitted
-  double strength_ {1.0};                         //!< Source strength
   UPtrSpace space_;                               //!< Spatial distribution
   UPtrAngle angle_;                               //!< Angular distribution
   UPtrDist energy_;                               //!< Energy distribution
   UPtrDist time_;                                 //!< Time distribution
-  DomainType domain_type_;                        //!< Domain type for rejection
-  std::unordered_set<int32_t> domain_ids_;        //!< Domains to reject from
 };
 
 //==============================================================================
@@ -107,9 +158,12 @@ public:
   explicit FileSource(const std::string& path);
 
   // Methods
-  SourceSite sample(uint64_t* seed) const override;
   void load_sites_from_file(
     const std::string& path); //!< Load source sites from file
+
+protected:
+  SourceSite sample(uint64_t* seed) const override;
+
 private:
   vector<SourceSite> sites_; //!< Source sites from a file
 };
@@ -124,16 +178,17 @@ public:
   CompiledSourceWrapper(pugi::xml_node node);
   ~CompiledSourceWrapper();
 
+  double strength() const override { return compiled_source_->strength(); }
+
+  void setup(const std::string& path, const std::string& parameters);
+
+protected:
   // Defer implementation to custom source library
   SourceSite sample(uint64_t* seed) const override
   {
     return compiled_source_->sample(seed);
   }
 
-  double strength() const override { return compiled_source_->strength(); }
-
-  void setup(const std::string& path, const std::string& parameters);
-
 private:
   void* shared_library_; //!< library from dlopen
   unique_ptr<Source> compiled_source_;
@@ -164,6 +219,9 @@ public:
     return sources_.size() == 1 ? sources_[0] : sources_[i];
   }
 
+protected:
+  bool constraints_applied() const override { return true; }
+
 private:
   // Data members
   unique_ptr<MeshSpatial> space_;           //!< Mesh spatial

include/openmc/tallies/filter.h

@@ -33,6 +33,7 @@ enum class FilterType {
   MATERIAL,
   MATERIALFROM,
   MESH,
+  MESHBORN,
   MESH_SURFACE,
   MU,
   PARTICLE,

include/openmc/tallies/filter_meshborn.h

@@ -0,0 +1,26 @@
+#ifndef OPENMC_TALLIES_FILTER_MESHBORN_H
+#define OPENMC_TALLIES_FILTER_MESHBORN_H
+
+#include <cstdint>
+
+#include "openmc/position.h"
+#include "openmc/tallies/filter_mesh.h"
+
+namespace openmc {
+
+class MeshBornFilter : public MeshFilter {
+public:
+  //----------------------------------------------------------------------------
+  // Methods
+
+  std::string type_str() const override { return "meshborn"; }
+  FilterType type() const override { return FilterType::MESHBORN; }
+
+  void get_all_bins(const Particle& p, TallyEstimator estimator,
+    FilterMatch& match) const override;
+
+  std::string text_label(int bin) const override;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_MESHBORN_H

include/openmc/tallies/trigger.h

@@ -23,6 +23,7 @@ enum class TriggerMetric {
 struct Trigger {
   TriggerMetric metric; //!< The type of uncertainty (e.g. std dev) measured
   double threshold;     //!< Uncertainty value below which trigger is satisfied
+  bool ignore_zeros;    //!< Whether to allow zero tally bins to be ignored
   int score_index;      //!< Index of the relevant score in the tally's arrays
 };
 

include/openmc/timer.h

@@ -31,6 +31,7 @@ extern Timer time_event_advance_particle;
 extern Timer time_event_surface_crossing;
 extern Timer time_event_collision;
 extern Timer time_event_death;
+extern Timer time_update_src;
 
 } // namespace simulation
 

openmc/arithmetic.py

@@ -54,8 +54,7 @@ class CrossScore:
         return str(other) == str(self)
 
     def __repr__(self):
-        return '({} {} {})'.format(self.left_score, self.binary_op,
-                                   self.right_score)
+        return f'({self.left_score} {self.binary_op} {self.right_score})'
 
     @property
     def left_score(self):
@@ -188,9 +187,9 @@ class CrossFilter:
 
     Parameters
     ----------
-    left_filter : Filter or CrossFilter
+    left_filter : openmc.Filter or CrossFilter
         The left filter in the outer product
-    right_filter : Filter or CrossFilter
+    right_filter : openmc.Filter or CrossFilter
         The right filter in the outer product
     binary_op : str
         The tally arithmetic binary operator (e.g., '+', '-', etc.) used to
@@ -200,9 +199,9 @@ class CrossFilter:
     ----------
     type : str
         The type of the crossfilter (e.g., 'energy / energy')
-    left_filter : Filter or CrossFilter
+    left_filter : openmc.Filter or CrossFilter
         The left filter in the outer product
-    right_filter : Filter or CrossFilter
+    right_filter : openmc.Filter or CrossFilter
         The right filter in the outer product
     binary_op : str
         The tally arithmetic binary operator (e.g., '+', '-', etc.) used to
@@ -271,7 +270,7 @@ class CrossFilter:
     def type(self):
         left_type = self.left_filter.type
         right_type = self.right_filter.type
-        return '({} {} {})'.format(left_type, self.binary_op, right_type)
+        return f'({left_type} {self.binary_op} {right_type})'
 
     @property
     def bins(self):
@@ -517,7 +516,7 @@ class AggregateFilter:
 
     Parameters
     ----------
-    aggregate_filter : Filter or CrossFilter
+    aggregate_filter : openmc.Filter or CrossFilter
         The filter included in the aggregation
     bins : Iterable of tuple
         The filter bins included in the aggregation
@@ -529,7 +528,7 @@ class AggregateFilter:
     ----------
     type : str
         The type of the aggregatefilter (e.g., 'sum(energy)', 'sum(cell)')
-    aggregate_filter : filter
+    aggregate_filter : openmc.Filter
         The filter included in the aggregation
     aggregate_op : str
         The tally aggregation operator (e.g., 'sum', 'avg', etc.) used

openmc/bounding_box.py

@@ -95,9 +95,23 @@ class BoundingBox:
         new |= other
         return new
 
-    def __contains__(self, point):
-        """Check whether or not a point is in the bounding box"""
-        return all(point > self.lower_left) and all(point < self.upper_right)
+    def __contains__(self, other):
+        """Check whether or not a point or another bounding box is in the bounding box.
+
+        For another bounding box to be in the parent it must lie fully inside of it.
+        """
+        # test for a single point
+        if isinstance(other, (tuple, list, np.ndarray)):
+            point = other
+            check_length("Point", point, 3, 3)
+            return all(point > self.lower_left) and all(point < self.upper_right)
+        elif isinstance(other, BoundingBox):
+            return all([p in self for p in [other.lower_left, other.upper_right]])
+        else:
+            raise TypeError(
+                f"Unable to determine if {other} is in the bounding box."
+                f" Expected a tuple or a bounding box, but {type(other)} given"
+            )
 
     @property
     def center(self) -> np.ndarray:

openmc/cell.py

@@ -343,8 +343,7 @@ class Cell(IDManagerMixin):
         if self.region is not None:
             return self.region.bounding_box
         else:
-            return BoundingBox(np.array([-np.inf, -np.inf, -np.inf]),
-                               np.array([np.inf, np.inf, np.inf]))
+            return BoundingBox.infinite()
 
     @property
     def num_instances(self):
@@ -427,15 +426,13 @@ class Cell(IDManagerMixin):
             instances
 
         """
+        if memo is None:
+            memo = set()
+        elif self in memo:
+            return {}
+        memo.add(self)
 
         cells = {}
-
-        if memo and self in memo:
-            return cells
-
-        if memo is not None:
-            memo.add(self)
-
         if self.fill_type in ('universe', 'lattice'):
             cells.update(self.fill.get_all_cells(memo))
 
@@ -466,7 +463,7 @@ class Cell(IDManagerMixin):
 
         return materials
 
-    def get_all_universes(self):
+    def get_all_universes(self, memo=None):
         """Return all universes that are contained within this one if any of
         its cells are filled with a universe or lattice.
 
@@ -477,18 +474,22 @@ class Cell(IDManagerMixin):
             :class:`Universe` instances
 
         """
+        if memo is None:
+            memo = set()
+        if self in memo:
+            return {}
+        memo.add(self)
 
         universes = {}
-
         if self.fill_type == 'universe':
             universes[self.fill.id] = self.fill
-            universes.update(self.fill.get_all_universes())
+            universes.update(self.fill.get_all_universes(memo))
         elif self.fill_type == 'lattice':
-            universes.update(self.fill.get_all_universes())
+            universes.update(self.fill.get_all_universes(memo))
 
         return universes
 
-    def clone(self, clone_materials=True, clone_regions=True,  memo=None):
+    def clone(self, clone_materials=True, clone_regions=True, memo=None):
         """Create a copy of this cell with a new unique ID, and clones
         the cell's region and fill.
 
@@ -594,7 +595,7 @@ class Cell(IDManagerMixin):
 
                # Make water blue
                water = openmc.Cell(fill=h2o)
-               universe.plot(..., colors={water: (0., 0., 1.))
+               water.plot(colors={water: (0., 0., 1.)})
         seed : int
             Seed for the random number generator
         openmc_exec : str
@@ -619,8 +620,10 @@ class Cell(IDManagerMixin):
 
         """
         # Create dummy universe but preserve used_ids
-        u = openmc.Universe(cells=[self], universe_id=openmc.Universe.next_id + 1)
+        next_id = openmc.Universe.next_id
+        u = openmc.Universe(cells=[self])
         openmc.Universe.used_ids.remove(u.id)
+        openmc.Universe.next_id = next_id
         return u.plot(*args, **kwargs)
 
     def create_xml_subelement(self, xml_element, memo=None):
@@ -677,10 +680,11 @@ class Cell(IDManagerMixin):
             # thus far.
             def create_surface_elements(node, element, memo=None):
                 if isinstance(node, Halfspace):
-                    if memo and node.surface in memo:
+                    if memo is None:
+                        memo = set()
+                    elif node.surface in memo:
                         return
-                    if memo is not None:
-                        memo.add(node.surface)
+                    memo.add(node.surface)
                     xml_element.append(node.surface.to_xml_element())
 
                 elif isinstance(node, Complement):

openmc/checkvalue.py

@@ -171,6 +171,29 @@ def check_length(name, value, length_min, length_max=None):
         raise ValueError(msg)
 
 
+def check_increasing(name: str, value, equality: bool = False):
+    """Ensure that a list's elements are strictly or loosely increasing.
+
+    Parameters
+    ----------
+    name : str
+        Description of value being checked
+    value : iterable
+        Object to check if increasing
+    equality : bool, optional
+        Whether equality is allowed. Defaults to False.
+
+    """
+    if equality:
+        if not np.all(np.diff(value) >= 0.0):
+            raise ValueError(f'Unable to set "{name}" to "{value}" since its '
+                             'elements must be increasing.')
+    elif not equality:
+        if not np.all(np.diff(value) > 0.0):
+            raise ValueError(f'Unable to set "{name}" to "{value}" since its '
+                             'elements must be strictly increasing.')
+
+
 def check_value(name, value, accepted_values):
     """Ensure that an object's value is contained in a set of acceptable values.
 

openmc/cmfd.py

@@ -135,7 +135,7 @@ class CMFDMesh:
         return outstr
 
     def _get_repr(self, list_var, label):
-        outstr = "\t{:<11} = ".format(label)
+        outstr = f"\t{label:<11} = "
         if list(list_var):
             outstr += ", ".join(str(i) for i in list_var)
         return outstr
@@ -242,9 +242,9 @@ class CMFDMesh:
 
         check_length('CMFD mesh grid', grid, grid_length)
         for i in range(grid_length):
-            check_type('CMFD mesh {}-grid'.format(dims[i]), grid[i], Iterable,
+            check_type(f'CMFD mesh {dims[i]}-grid', grid[i], Iterable,
                        Real)
-            check_greater_than('CMFD mesh {}-grid length'.format(dims[i]),
+            check_greater_than(f'CMFD mesh {dims[i]}-grid length',
                                len(grid[i]), 1)
         self._grid = [np.array(g) for g in grid]
         self._display_mesh_warning('rectilinear', 'CMFD mesh grid')
@@ -612,7 +612,7 @@ class CMFDRun:
         for key, value in display.items():
             check_value('display key', key,
                         ('balance', 'entropy', 'dominance', 'source'))
-            check_type("display['{}']".format(key), value, bool)
+            check_type(f"display['{key}']", value, bool)
             self._display[key] = value
 
     @downscatter.setter
@@ -928,7 +928,7 @@ class CMFDRun:
             with h5py.File(filename, 'a') as f:
                 if 'cmfd' not in f:
                     if openmc.lib.settings.verbosity >= 5:
-                        print(' Writing CMFD data to {}...'.format(filename))
+                        print(f' Writing CMFD data to {filename}...')
                         sys.stdout.flush()
                     cmfd_group = f.create_group("cmfd")
                     cmfd_group.attrs['cmfd_on'] = self._cmfd_on
@@ -982,14 +982,16 @@ class CMFDRun:
         temp_data = np.ones(len(loss_row))
         temp_loss = sparse.csr_matrix((temp_data, (loss_row, loss_col)),
                                       shape=(n, n))
+        temp_loss.sort_indices()
 
         # Pass coremap as 1-d array of 32-bit integers
         coremap = np.swapaxes(self._coremap, 0, 2).flatten().astype(np.int32)
 
-        args = temp_loss.indptr, len(temp_loss.indptr), \
-            temp_loss.indices, len(temp_loss.indices), n, \
+        return openmc.lib._dll.openmc_initialize_linsolver(
+            temp_loss.indptr.astype(np.int32), len(temp_loss.indptr),
+            temp_loss.indices.astype(np.int32), len(temp_loss.indices), n,
             self._spectral, coremap, self._use_all_threads
-        return openmc.lib._dll.openmc_initialize_linsolver(*args)
+        )
 
     def _write_cmfd_output(self):
         """Write CMFD output to buffer at the end of each batch"""
@@ -1132,12 +1134,12 @@ class CMFDRun:
         with h5py.File(filename, 'r') as f:
             if 'cmfd' not in f:
                 raise OpenMCError('Could not find CMFD parameters in ',
-                                  'file {}'.format(filename))
+                                  f'file {filename}')
             else:
                 # Overwrite CMFD values from statepoint
                 if (openmc.lib.master() and
                         openmc.lib.settings.verbosity >= 5):
-                    print(' Loading CMFD data from {}...'.format(filename))
+                    print(f' Loading CMFD data from {filename}...')
                     sys.stdout.flush()
                 cmfd_group = f['cmfd']
 
@@ -1407,8 +1409,7 @@ class CMFDRun:
                 # Get all data entries for particular row in matrix
                 data = matrix.data[matrix.indptr[row]:matrix.indptr[row+1]]
                 for i in range(len(cols)):
-                    fh.write('{:3d}, {:3d}, {:0.8f}\n'.format(
-                        row, cols[i], data[i]))
+                    fh.write(f'{row:3d}, {cols[i]:3d}, {data[i]:0.8f}\n')
 
         # Save matrix in scipy format
         sparse.save_npz(base_filename, matrix)
@@ -1585,6 +1586,7 @@ class CMFDRun:
         loss_row = self._loss_row
         loss_col = self._loss_col
         loss = sparse.csr_matrix((data, (loss_row, loss_col)), shape=(n, n))
+        loss.sort_indices()
         return loss
 
     def _build_prod_matrix(self, adjoint):
@@ -1611,6 +1613,7 @@ class CMFDRun:
         prod_row = self._prod_row
         prod_col = self._prod_col
         prod = sparse.csr_matrix((data, (prod_row, prod_col)), shape=(n, n))
+        prod.sort_indices()
         return prod
 
     def _execute_power_iter(self, loss, prod):
@@ -2307,8 +2310,8 @@ class CMFDRun:
                                    constant_values=_CMFD_NOACCEL)[:,:,1:]
 
         # Create empty row and column vectors to store for loss matrix
-        row = np.array([])
-        col = np.array([])
+        row = np.array([], dtype=int)
+        col = np.array([], dtype=int)
 
         # Store all indices used to populate production and loss matrix
         is_accel = self._coremap != _CMFD_NOACCEL

openmc/data/ace.py

@@ -143,7 +143,7 @@ def ascii_to_binary(ascii_file, binary_file):
             # that XSS will start at the second record
             nxs = [int(x) for x in ' '.join(lines[idx + 6:idx + 8]).split()]
             jxs = [int(x) for x in ' '.join(lines[idx + 8:idx + 12]).split()]
-            binary_file.write(struct.pack(str('=16i32i{}x'.format(record_length - 500)),
+            binary_file.write(struct.pack(str(f'=16i32i{record_length - 500}x'),
                                           *(nxs + jxs)))
 
             # Read/write XSS array. Null bytes are added to form a complete record
@@ -152,8 +152,7 @@ def ascii_to_binary(ascii_file, binary_file):
             start = idx + _ACE_HEADER_SIZE
             xss = np.fromstring(' '.join(lines[start:start + n_lines]), sep=' ')
             extra_bytes = record_length - ((len(xss)*8 - 1) % record_length + 1)
-            binary_file.write(struct.pack(str('={}d{}x'.format(
-                nxs[0], extra_bytes)), *xss))
+            binary_file.write(struct.pack(str(f'={nxs[0]}d{extra_bytes}x'), *xss))
 
             # Advance to next table in file
             idx += _ACE_HEADER_SIZE + n_lines
@@ -184,8 +183,7 @@ def get_table(filename, name=None):
         if lib.tables:
             return lib.tables[0]
         else:
-            raise ValueError('Could not find ACE table with name: {}'
-                             .format(name))
+            raise ValueError(f'Could not find ACE table with name: {name}')
 
 
 # The beginning of an ASCII ACE file consists of 12 lines that include the name,
@@ -295,14 +293,14 @@ class Library(EqualityMixin):
 
             if verbose:
                 kelvin = round(temperature * EV_PER_MEV / K_BOLTZMANN)
-                print("Loading nuclide {} at {} K".format(name, kelvin))
+                print(f"Loading nuclide {name} at {kelvin} K")
 
             # Read JXS
             jxs = list(struct.unpack(str('=32i'), ace_file.read(128)))
 
             # Read XSS
             ace_file.seek(start_position + recl_length)
-            xss = list(struct.unpack(str('={}d'.format(length)),
+            xss = list(struct.unpack(str(f'={length}d'),
                                      ace_file.read(length*8)))
 
             # Insert zeros at beginning of NXS, JXS, and XSS arrays so that the
@@ -393,7 +391,7 @@ class Library(EqualityMixin):
 
             if verbose:
                 kelvin = round(temperature * EV_PER_MEV / K_BOLTZMANN)
-                print("Loading nuclide {} at {} K".format(name, kelvin))
+                print(f"Loading nuclide {name} at {kelvin} K")
 
             # Insert zeros at beginning of NXS, JXS, and XSS arrays so that the
             # indexing will be the same as Fortran. This makes it easier to
@@ -455,8 +453,7 @@ class TableType(enum.Enum):
         for member in cls:
             if suffix.endswith(member.value):
                 return member
-        raise ValueError("Suffix '{}' has no corresponding ACE table type."
-                         .format(suffix))
+        raise ValueError(f"Suffix '{suffix}' has no corresponding ACE table type.")
 
 
 class Table(EqualityMixin):
@@ -507,7 +504,7 @@ class Table(EqualityMixin):
         return TableType.from_suffix(xs[-1])
 
     def __repr__(self):
-        return "<ACE Table: {}>".format(self.name)
+        return f"<ACE Table: {self.name}>"
 
 
 def get_libraries_from_xsdir(path):

openmc/data/angle_energy.py

@@ -112,7 +112,6 @@ class AngleEnergy(EqualityMixin, ABC):
             distribution = openmc.data.NBodyPhaseSpace.from_ace(
                 ace, idx, rx.q_value)
         else:
-            raise ValueError("Unsupported ACE secondary energy "
-                             "distribution law {}".format(law))
+            raise ValueError(f"Unsupported ACE secondary energy distribution law {law}")
 
         return distribution

openmc/data/correlated.py

@@ -113,7 +113,7 @@ class CorrelatedAngleEnergy(AngleEnergy):
             HDF5 group to write to
 
         """
-        group.attrs['type'] = np.string_(self._name)
+        group.attrs['type'] = np.bytes_(self._name)
 
         dset = group.create_dataset('energy', data=self.energy)
         dset.attrs['interpolation'] = np.vstack((self.breakpoints,

openmc/data/decay.py

@@ -128,7 +128,7 @@ class FissionProductYields(EqualityMixin):
                     isomeric_state = int(values[4*j + 1])
                     name = ATOMIC_SYMBOL[Z] + str(A)
                     if isomeric_state > 0:
-                        name += '_m{}'.format(isomeric_state)
+                        name += f'_m{isomeric_state}'
                     yield_j = ufloat(values[4*j + 2], values[4*j + 3])
                     yields[name] = yield_j
 
@@ -257,9 +257,9 @@ class DecayMode(EqualityMixin):
                         Z += delta_Z
 
         if self._daughter_state > 0:
-            return '{}{}_m{}'.format(ATOMIC_SYMBOL[Z], A, self._daughter_state)
+            return f'{ATOMIC_SYMBOL[Z]}{A}_m{self._daughter_state}'
         else:
-            return '{}{}'.format(ATOMIC_SYMBOL[Z], A)
+            return f'{ATOMIC_SYMBOL[Z]}{A}'
 
     @property
     def parent(self):
@@ -350,10 +350,9 @@ class Decay(EqualityMixin):
         self.nuclide['mass_number'] = A
         self.nuclide['isomeric_state'] = metastable
         if metastable > 0:
-            self.nuclide['name'] = '{}{}_m{}'.format(ATOMIC_SYMBOL[Z], A,
-                                                     metastable)
+            self.nuclide['name'] = f'{ATOMIC_SYMBOL[Z]}{A}_m{metastable}'
         else:
-            self.nuclide['name'] = '{}{}'.format(ATOMIC_SYMBOL[Z], A)
+            self.nuclide['name'] = f'{ATOMIC_SYMBOL[Z]}{A}'
         self.nuclide['mass'] = items[1]  # AWR
         self.nuclide['excited_state'] = items[2]  # State of the original nuclide
         self.nuclide['stable'] = (items[4] == 1)  # Nucleus stability flag

openmc/data/effective_dose/dose.py

@@ -60,7 +60,7 @@ def dose_coefficients(particle, geometry='AP'):
     # Get all data for selected particle
     data = _DOSE_ICRP116.get(particle)
     if data is None:
-        raise ValueError("{} has no effective dose data".format(particle))
+        raise ValueError(f"{particle} has no effective dose data")
 
     # Determine index for selected geometry
     if particle in ('neutron', 'photon', 'proton', 'photon kerma'):

openmc/data/endf.c

@@ -14,7 +14,7 @@
 
 double cfloat_endf(const char* buffer, int n)
 {
-  char arr[12]; // 11 characters plus a null terminator
+  char arr[13]; // 11 characters plus e and a null terminator
   int j = 0; // current position in arr
   int found_significand = 0;
   int found_exponent = 0;

openmc/data/endf.py

@@ -449,8 +449,7 @@ class Evaluation:
 
     def __repr__(self):
         name = self.target['zsymam'].replace(' ', '')
-        return '<{} for {} {}>'.format(self.info['sublibrary'], name,
-                                       self.info['library'])
+        return f"<{self.info['sublibrary']} for {name} {self.info['library']}>"
 
     def _read_header(self):
         file_obj = io.StringIO(self.section[1, 451])

openmc/data/energy_distribution.py

@@ -53,8 +53,7 @@ class EnergyDistribution(EqualityMixin, ABC):
         elif energy_type == 'continuous':
             return ContinuousTabular.from_hdf5(group)
         else:
-            raise ValueError("Unknown energy distribution type: {}"
-                             .format(energy_type))
+            raise ValueError(f"Unknown energy distribution type: {energy_type}")
 
     @staticmethod
     def from_endf(file_obj, params):
@@ -277,7 +276,7 @@ class MaxwellEnergy(EnergyDistribution):
 
         """
 
-        group.attrs['type'] = np.string_('maxwell')
+        group.attrs['type'] = np.bytes_('maxwell')
         group.attrs['u'] = self.u
         self.theta.to_hdf5(group, 'theta')
 
@@ -410,7 +409,7 @@ class Evaporation(EnergyDistribution):
 
         """
 
-        group.attrs['type'] = np.string_('evaporation')
+        group.attrs['type'] = np.bytes_('evaporation')
         group.attrs['u'] = self.u
         self.theta.to_hdf5(group, 'theta')
 
@@ -556,7 +555,7 @@ class WattEnergy(EnergyDistribution):
 
         """
 
-        group.attrs['type'] = np.string_('watt')
+        group.attrs['type'] = np.bytes_('watt')
         group.attrs['u'] = self.u
         self.a.to_hdf5(group, 'a')
         self.b.to_hdf5(group, 'b')
@@ -728,7 +727,7 @@ class MadlandNix(EnergyDistribution):
 
         """
 
-        group.attrs['type'] = np.string_('madland-nix')
+        group.attrs['type'] = np.bytes_('madland-nix')
         group.attrs['efl'] = self.efl
         group.attrs['efh'] = self.efh
         self.tm.to_hdf5(group)
@@ -846,7 +845,7 @@ class DiscretePhoton(EnergyDistribution):
 
         """
 
-        group.attrs['type'] = np.string_('discrete_photon')
+        group.attrs['type'] = np.bytes_('discrete_photon')
         group.attrs['primary_flag'] = self.primary_flag
         group.attrs['energy'] = self.energy
         group.attrs['atomic_weight_ratio'] = self.atomic_weight_ratio
@@ -945,7 +944,7 @@ class LevelInelastic(EnergyDistribution):
 
         """
 
-        group.attrs['type'] = np.string_('level')
+        group.attrs['type'] = np.bytes_('level')
         group.attrs['threshold'] = self.threshold
         group.attrs['mass_ratio'] = self.mass_ratio
 
@@ -1074,7 +1073,7 @@ class ContinuousTabular(EnergyDistribution):
 
         """
 
-        group.attrs['type'] = np.string_('continuous')
+        group.attrs['type'] = np.bytes_('continuous')
 
         dset = group.create_dataset('energy', data=self.energy)
         dset.attrs['interpolation'] = np.vstack((self.breakpoints,

openmc/data/function.py

@@ -364,7 +364,7 @@ class Tabulated1D(Function1D):
         """
         dataset = group.create_dataset(name, data=np.vstack(
             [self.x, self.y]))
-        dataset.attrs['type'] = np.string_(type(self).__name__)
+        dataset.attrs['type'] = np.bytes_(type(self).__name__)
         dataset.attrs['breakpoints'] = self.breakpoints
         dataset.attrs['interpolation'] = self.interpolation
 
@@ -460,7 +460,7 @@ class Polynomial(np.polynomial.Polynomial, Function1D):
 
         """
         dataset = group.create_dataset(name, data=self.coef)
-        dataset.attrs['type'] = np.string_(type(self).__name__)
+        dataset.attrs['type'] = np.bytes_(type(self).__name__)
 
     @classmethod
     def from_hdf5(cls, dataset):
@@ -592,7 +592,7 @@ class Sum(Function1D):
 
         """
         sum_group = group.create_group(name)
-        sum_group.attrs['type'] = np.string_(type(self).__name__)
+        sum_group.attrs['type'] = np.bytes_(type(self).__name__)
         sum_group.attrs['n'] = len(self.functions)
         for i, f in enumerate(self.functions):
             f.to_hdf5(sum_group, f'func_{i+1}')

openmc/data/kalbach_mann.py

@@ -366,7 +366,7 @@ class KalbachMann(AngleEnergy):
             HDF5 group to write to
 
         """
-        group.attrs['type'] = np.string_('kalbach-mann')
+        group.attrs['type'] = np.bytes_('kalbach-mann')
 
         dset = group.create_dataset('energy', data=self.energy)
         dset.attrs['interpolation'] = np.vstack((self.breakpoints,

openmc/data/library.py

@@ -93,8 +93,7 @@ class DataLibrary(list):
                 materials = list(h5file)
         else:
             raise ValueError(
-                "File type {} not supported by {}"
-                .format(path.name, self.__class__.__name__))
+                f"File type {path.name} not supported by {self.__class__.__name__}")
 
         library = {'path': str(path), 'type': filetype, 'materials': materials}
         self.append(library)

openmc/data/multipole.py

@@ -194,9 +194,8 @@ def _vectfit_xs(energy, ce_xs, mts, rtol=1e-3, atol=1e-5, orders=None,
         test_xs_ref[i] = np.interp(test_energy, energy, ce_xs[i])
 
     if log:
-        print("  energy: {:.3e} to {:.3e} eV ({} points)".format(
-              energy[0], energy[-1], ne))
-        print("  error tolerance: rtol={}, atol={}".format(rtol, atol))
+        print(f"  energy: {energy[0]:.3e} to {energy[-1]:.3e} eV ({ne} points)")
+        print(f"  error tolerance: rtol={rtol}, atol={atol}")
 
     # transform xs (sigma) and energy (E) to f (sigma*E) and s (sqrt(E)) to be
     # compatible with the multipole representation
@@ -230,8 +229,8 @@ def _vectfit_xs(energy, ce_xs, mts, rtol=1e-3, atol=1e-5, orders=None,
         orders = list(range(lowest_order, highest_order + 1, 2))
 
     if log:
-        print("Found {} peaks".format(n_peaks))
-        print("Fitting orders from {} to {}".format(orders[0], orders[-1]))
+        print(f"Found {n_peaks} peaks")
+        print(f"Fitting orders from {orders[0]} to {orders[-1]}")
 
     # perform VF with increasing orders
     found_ideal = False
@@ -239,7 +238,7 @@ def _vectfit_xs(energy, ce_xs, mts, rtol=1e-3, atol=1e-5, orders=None,
     best_quality = best_ratio = -np.inf
     for i, order in enumerate(orders):
         if log:
-            print("Order={}({}/{})".format(order, i, len(orders)))
+            print(f"Order={order}({i}/{len(orders)})")
         # initial guessed poles
         poles_r = np.linspace(s[0], s[-1], order//2)
         poles = poles_r + poles_r*0.01j
@@ -249,7 +248,7 @@ def _vectfit_xs(energy, ce_xs, mts, rtol=1e-3, atol=1e-5, orders=None,
         # fitting iteration
         for i_vf in range(n_vf_iter):
             if log >= DETAILED_LOGGING:
-                print("VF iteration {}/{}".format(i_vf + 1, n_vf_iter))
+                print(f"VF iteration {i_vf + 1}/{n_vf_iter}")
 
             # call vf
             poles, residues, cf, f_fit, rms = vf.vectfit(f, s, poles, weight)
@@ -268,7 +267,7 @@ def _vectfit_xs(energy, ce_xs, mts, rtol=1e-3, atol=1e-5, orders=None,
             # re-calculate residues if poles changed
             if n_real_poles > 0:
                 if log >= DETAILED_LOGGING:
-                    print("  # real poles: {}".format(n_real_poles))
+                    print(f"  # real poles: {n_real_poles}")
                 new_poles, residues, cf, f_fit, rms = \
                       vf.vectfit(f, s, new_poles, weight, skip_pole=True)
 
@@ -296,10 +295,10 @@ def _vectfit_xs(energy, ce_xs, mts, rtol=1e-3, atol=1e-5, orders=None,
                 quality = -np.inf
 
             if log >= DETAILED_LOGGING:
-                print("  # poles: {}".format(new_poles.size))
-                print("  Max relative error: {:.3f}%".format(maxre*100))
-                print("  Satisfaction: {:.1f}%, {:.1f}%".format(ratio*100, ratio2*100))
-                print("  Quality: {:.2f}".format(quality))
+                print(f"  # poles: {new_poles.size}")
+                print(f"  Max relative error: {maxre * 100:.3f}%")
+                print(f"  Satisfaction: {ratio * 100:.1f}%, {ratio2 * 100:.1f}%")
+                print(f"  Quality: {quality:.2f}")
 
             if quality > best_quality:
                 if log >= DETAILED_LOGGING:
@@ -354,7 +353,7 @@ def _vectfit_xs(energy, ce_xs, mts, rtol=1e-3, atol=1e-5, orders=None,
     mp_residues = np.concatenate((best_residues[:, real_idx],
                                   best_residues[:, conj_idx]*2), axis=1)/1j
     if log:
-        print("Final number of poles: {}".format(mp_poles.size))
+        print(f"Final number of poles: {mp_poles.size}")
 
     if path_out:
         if not os.path.exists(path_out):
@@ -378,14 +377,14 @@ def _vectfit_xs(energy, ce_xs, mts, rtol=1e-3, atol=1e-5, orders=None,
             ax2.set_ylabel('relative error', color='r')
             ax2.tick_params('y', colors='r')
 
-            plt.title("MT {} vector fitted with {} poles".format(mt, mp_poles.size))
+            plt.title(f"MT {mt} vector fitted with {mp_poles.size} poles")
             fig.tight_layout()
             fig_file = os.path.join(path_out, "{:.0f}-{:.0f}_MT{}.png".format(
                                     energy[0], energy[-1], mt))
             plt.savefig(fig_file)
             plt.close()
             if log:
-                print("Saved figure: {}".format(fig_file))
+                print(f"Saved figure: {fig_file}")
 
     return (mp_poles, mp_residues)
 
@@ -423,7 +422,7 @@ def vectfit_nuclide(endf_file, njoy_error=5e-4, vf_pieces=None,
 
     # make 0K ACE data using njoy
     if log:
-        print("Running NJOY to get 0K point-wise data (error={})...".format(njoy_error))
+        print(f"Running NJOY to get 0K point-wise data (error={njoy_error})...")
 
     nuc_ce = IncidentNeutron.from_njoy(endf_file, temperatures=[0.0],
              error=njoy_error, broadr=False, heatr=False, purr=False)
@@ -477,9 +476,8 @@ def vectfit_nuclide(endf_file, njoy_error=5e-4, vf_pieces=None,
         mts = [2, 27]
 
     if log:
-        print("  MTs: {}".format(mts))
-        print("  Energy range: {:.3e} to {:.3e} eV ({} points)".format(
-              E_min, E_max, n_points))
+        print(f"  MTs: {mts}")
+        print(f"  Energy range: {E_min:.3e} to {E_max:.3e} eV ({n_points} points)")
 
     # ======================================================================
     # PERFORM VECTOR FITTING
@@ -500,7 +498,7 @@ def vectfit_nuclide(endf_file, njoy_error=5e-4, vf_pieces=None,
     # VF piece by piece
     for i_piece in range(vf_pieces):
         if log:
-            print("Vector fitting piece {}/{}...".format(i_piece + 1, vf_pieces))
+            print(f"Vector fitting piece {i_piece + 1}/{vf_pieces}...")
         # start E of this piece
         e_bound = (sqrt(E_min) + piece_width*(i_piece-0.5))**2
         if i_piece == 0 or sqrt(alpha*e_bound) < 4.0:
@@ -534,12 +532,12 @@ def vectfit_nuclide(endf_file, njoy_error=5e-4, vf_pieces=None,
         if not os.path.exists(path_out):
             os.makedirs(path_out)
         if not mp_filename:
-            mp_filename = "{}_mp.pickle".format(nuc_ce.name)
+            mp_filename = f"{nuc_ce.name}_mp.pickle"
         mp_filename = os.path.join(path_out, mp_filename)
         with open(mp_filename, 'wb') as f:
             pickle.dump(mp_data, f)
         if log:
-            print("Dumped multipole data to file: {}".format(mp_filename))
+            print(f"Dumped multipole data to file: {mp_filename}")
 
     return mp_data
 
@@ -605,9 +603,8 @@ def _windowing(mp_data, n_cf, rtol=1e-3, atol=1e-5, n_win=None, spacing=None,
 
     if log:
         print("Windowing:")
-        print("  config: # windows={}, spacing={}, CF order={}".format(
-               n_win, spacing, n_cf))
-        print("  error tolerance: rtol={}, atol={}".format(rtol, atol))
+        print(f"  config: # windows={n_win}, spacing={spacing}, CF order={n_cf}")
+        print(f"  error tolerance: rtol={rtol}, atol={atol}")
 
     # sort poles (and residues) by the real component of the pole
     for ip in range(n_pieces):
@@ -623,7 +620,7 @@ def _windowing(mp_data, n_cf, rtol=1e-3, atol=1e-5, n_win=None, spacing=None,
     win_data = []
     for iw in range(n_win):
         if log >= DETAILED_LOGGING:
-            print("Processing window {}/{}...".format(iw + 1, n_win))
+            print(f"Processing window {iw + 1}/{n_win}...")
 
         # inner window boundaries
         inbegin = sqrt(E_min) + spacing * iw
@@ -658,7 +655,7 @@ def _windowing(mp_data, n_cf, rtol=1e-3, atol=1e-5, n_win=None, spacing=None,
         lp = rp = center_pole_ind
         while True:
             if log >= DETAILED_LOGGING:
-                print("Trying poles {} to {}".format(lp, rp))
+                print(f"Trying poles {lp} to {rp}")
 
             # calculate the cross sections contributed by the windowed poles
             if rp > lp:
@@ -1108,7 +1105,7 @@ class WindowedMultipole(EqualityMixin):
         for n_w in np.unique(np.linspace(n_win_min, n_win_max, 20, dtype=int)):
             for n_cf in range(10, 1, -1):
                 if log:
-                    print("Testing N_win={} N_cf={}".format(n_w, n_cf))
+                    print(f"Testing N_win={n_w} N_cf={n_cf}")
 
                 # update arguments dictionary
                 kwargs.update(n_win=n_w, n_cf=n_cf)
@@ -1169,10 +1166,11 @@ class WindowedMultipole(EqualityMixin):
         sqrtE = sqrt(E)
         invE = 1.0 / E
 
-        # Locate us.  The i_window calc omits a + 1 present in F90 because of
-        # the 1-based vs. 0-based indexing.  Similarly startw needs to be
-        # decreased by 1.  endw does not need to be decreased because
-        # range(startw, endw) does not include endw.
+        # Locate us.  The i_window calc omits a + 1 present from the legacy
+        # Fortran version of OpenMC because of the 1-based vs. 0-based
+        # indexing.  Similarly startw needs to be decreased by 1.  endw does
+        # not need to be decreased because range(startw, endw) does not include
+        # endw.
         i_window = min(self.n_windows - 1,
                        int(np.floor((sqrtE - sqrt(self.E_min)) / self.spacing)))
         startw = self.windows[i_window, 0] - 1
@@ -1273,7 +1271,7 @@ class WindowedMultipole(EqualityMixin):
 
         # Open file and write version.
         with h5py.File(str(path), mode, libver=libver) as f:
-            f.attrs['filetype'] = np.string_('data_wmp')
+            f.attrs['filetype'] = np.bytes_('data_wmp')
             f.attrs['version'] = np.array(WMP_VERSION)
 
             g = f.create_group(self.name)

openmc/data/nbody.py

@@ -91,7 +91,7 @@ class NBodyPhaseSpace(AngleEnergy):
             HDF5 group to write to
 
         """
-        group.attrs['type'] = np.string_('nbody')
+        group.attrs['type'] = np.bytes_('nbody')
         group.attrs['total_mass'] = self.total_mass
         group.attrs['n_particles'] = self.n_particles
         group.attrs['atomic_weight_ratio'] = self.atomic_weight_ratio

openmc/data/neutron.py

@@ -122,10 +122,10 @@ class IncidentNeutron(EqualityMixin):
             if len(mts) > 0:
                 return self._get_redundant_reaction(mt, mts)
             else:
-                raise KeyError('No reaction with MT={}.'.format(mt))
+                raise KeyError(f'No reaction with MT={mt}.')
 
     def __repr__(self):
-        return "<IncidentNeutron: {}>".format(self.name)
+        return f"<IncidentNeutron: {self.name}>"
 
     def __iter__(self):
         return iter(self.reactions.values())
@@ -231,7 +231,7 @@ class IncidentNeutron(EqualityMixin):
 
     @property
     def temperatures(self):
-        return ["{}K".format(int(round(kT / K_BOLTZMANN))) for kT in self.kTs]
+        return [f"{int(round(kT / K_BOLTZMANN))}K" for kT in self.kTs]
 
     @property
     def atomic_symbol(self):
@@ -261,7 +261,7 @@ class IncidentNeutron(EqualityMixin):
         # Check if temprature already exists
         strT = data.temperatures[0]
         if strT in self.temperatures:
-            warn('Cross sections at T={} already exist.'.format(strT))
+            warn(f'Cross sections at T={strT} already exist.')
             return
 
         # Check that name matches
@@ -425,7 +425,7 @@ class IncidentNeutron(EqualityMixin):
 
         # Open file and write version
         with h5py.File(str(path), mode, libver=libver) as f:
-            f.attrs['filetype'] = np.string_('data_neutron')
+            f.attrs['filetype'] = np.bytes_('data_neutron')
             f.attrs['version'] = np.array(HDF5_VERSION)
 
             # Write basic data
@@ -461,7 +461,7 @@ class IncidentNeutron(EqualityMixin):
                     if not (photon_rx or rx.mt in keep_mts):
                         continue
 
-                rx_group = rxs_group.create_group('reaction_{:03}'.format(rx.mt))
+                rx_group = rxs_group.create_group(f'reaction_{rx.mt:03}')
                 rx.to_hdf5(rx_group)
 
                 # Write total nu data if available
@@ -593,7 +593,7 @@ class IncidentNeutron(EqualityMixin):
         zaid, xs = ace.name.split('.')
         if not xs.endswith('c'):
             raise TypeError(
-                "{} is not a continuous-energy neutron ACE table.".format(ace))
+                f"{ace} is not a continuous-energy neutron ACE table.")
         name, element, Z, mass_number, metastable = \
             get_metadata(int(zaid), metastable_scheme)
 
@@ -732,9 +732,9 @@ class IncidentNeutron(EqualityMixin):
         # Determine name
         element = ATOMIC_SYMBOL[atomic_number]
         if metastable > 0:
-            name = '{}{}_m{}'.format(element, mass_number, metastable)
+            name = f'{element}{mass_number}_m{metastable}'
         else:
-            name = '{}{}'.format(element, mass_number)
+            name = f'{element}{mass_number}'
 
         # Instantiate incident neutron data
         data = cls(name, atomic_number, mass_number, metastable,

openmc/data/njoy.py

@@ -188,7 +188,7 @@ def run(commands, tapein, tapeout, input_filename=None, stdout=False,
     with tempfile.TemporaryDirectory() as tmpdir:
         # Copy evaluations to appropriates 'tapes'
         for tape_num, filename in tapein.items():
-            tmpfilename = os.path.join(tmpdir, 'tape{}'.format(tape_num))
+            tmpfilename = os.path.join(tmpdir, f'tape{tape_num}')
             shutil.copy(str(filename), tmpfilename)
 
         # Start up NJOY process
@@ -216,7 +216,7 @@ def run(commands, tapein, tapeout, input_filename=None, stdout=False,
 
         # Copy output files back to original directory
         for tape_num, filename in tapeout.items():
-            tmpfilename = os.path.join(tmpdir, 'tape{}'.format(tape_num))
+            tmpfilename = os.path.join(tmpdir, f'tape{tape_num}')
             if os.path.isfile(tmpfilename):
                 shutil.move(tmpfilename, str(filename))
 
@@ -317,7 +317,7 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None,
     else:
         output_dir = Path(output_dir)
         if not output_dir.is_dir():
-            raise IOError("{} is not a directory".format(output_dir))
+            raise IOError(f"{output_dir} is not a directory")
 
     ev = evaluation if evaluation is not None else endf.Evaluation(filename)
     mat = ev.material
@@ -389,12 +389,12 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None,
             # Extend input with an ACER run for each temperature
             nace = nacer_in + 1 + 2*i
             ndir = nace + 1
-            ext = '{:02}'.format(i + 1)
+            ext = f'{i + 1:02}'
             commands += _TEMPLATE_ACER.format(**locals())
 
             # Indicate tapes to save for each ACER run
-            tapeout[nace] = output_dir / "ace_{:.1f}".format(temperature)
-            tapeout[ndir] = output_dir / "xsdir_{:.1f}".format(temperature)
+            tapeout[nace] = output_dir / f"ace_{temperature:.1f}"
+            tapeout[ndir] = output_dir / f"xsdir_{temperature:.1f}"
     commands += 'stop\n'
     run(commands, tapein, tapeout, **kwargs)
 
@@ -404,7 +404,7 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None,
         with ace.open('w') as ace_file, xsdir.open('w') as xsdir_file:
             for temperature in temperatures:
                 # Get contents of ACE file
-                text = (output_dir / "ace_{:.1f}".format(temperature)).read_text()
+                text = (output_dir / f"ace_{temperature:.1f}").read_text()
 
                 # If the target is metastable, make sure that ZAID in the ACE
                 # file reflects this by adding 400
@@ -417,13 +417,13 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None,
                 ace_file.write(text)
 
                 # Concatenate into destination xsdir file
-                xsdir_in = output_dir / "xsdir_{:.1f}".format(temperature)
+                xsdir_in = output_dir / f"xsdir_{temperature:.1f}"
                 xsdir_file.write(xsdir_in.read_text())
 
         # Remove ACE/xsdir files for each temperature
         for temperature in temperatures:
-            (output_dir / "ace_{:.1f}".format(temperature)).unlink()
-            (output_dir / "xsdir_{:.1f}".format(temperature)).unlink()
+            (output_dir / f"ace_{temperature:.1f}").unlink()
+            (output_dir / f"xsdir_{temperature:.1f}").unlink()
 
 
 def make_ace_thermal(filename, filename_thermal, temperatures=None,
@@ -480,7 +480,7 @@ def make_ace_thermal(filename, filename_thermal, temperatures=None,
     else:
         output_dir = Path(output_dir)
         if not output_dir.is_dir():
-            raise IOError("{} is not a directory".format(output_dir))
+            raise IOError(f"{output_dir} is not a directory")
 
     ev = evaluation if evaluation is not None else endf.Evaluation(filename)
     mat = ev.material
@@ -581,12 +581,12 @@ def make_ace_thermal(filename, filename_thermal, temperatures=None,
         # Extend input with an ACER run for each temperature
         nace = nthermal_acer_in + 1 + 2*i
         ndir = nace + 1
-        ext = '{:02}'.format(i + 1)
+        ext = f'{i + 1:02}'
         commands += _THERMAL_TEMPLATE_ACER.format(**locals())
 
         # Indicate tapes to save for each ACER run
-        tapeout[nace] = output_dir / "ace_{:.1f}".format(temperature)
-        tapeout[ndir] = output_dir / "xsdir_{:.1f}".format(temperature)
+        tapeout[nace] = output_dir / f"ace_{temperature:.1f}"
+        tapeout[ndir] = output_dir / f"xsdir_{temperature:.1f}"
     commands += 'stop\n'
     run(commands, tapein, tapeout, **kwargs)
 
@@ -595,13 +595,13 @@ def make_ace_thermal(filename, filename_thermal, temperatures=None,
     with ace.open('w') as ace_file, xsdir.open('w') as xsdir_file:
         # Concatenate ACE and xsdir files together
         for temperature in temperatures:
-            ace_in = output_dir / "ace_{:.1f}".format(temperature)
+            ace_in = output_dir / f"ace_{temperature:.1f}"
             ace_file.write(ace_in.read_text())
 
-            xsdir_in = output_dir / "xsdir_{:.1f}".format(temperature)
+            xsdir_in = output_dir / f"xsdir_{temperature:.1f}"
             xsdir_file.write(xsdir_in.read_text())
 
     # Remove ACE/xsdir files for each temperature
     for temperature in temperatures:
-        (output_dir / "ace_{:.1f}".format(temperature)).unlink()
-        (output_dir / "xsdir_{:.1f}".format(temperature)).unlink()
+        (output_dir / f"ace_{temperature:.1f}").unlink()
+        (output_dir / f"xsdir_{temperature:.1f}").unlink()

openmc/data/photon.py

@@ -451,10 +451,10 @@ class IncidentPhoton(EqualityMixin):
         if mt in self.reactions:
             return self.reactions[mt]
         else:
-            raise KeyError('No reaction with MT={}.'.format(mt))
+            raise KeyError(f'No reaction with MT={mt}.')
 
     def __repr__(self):
-        return "<IncidentPhoton: {}>".format(self.name)
+        return f"<IncidentPhoton: {self.name}>"
 
     def __iter__(self):
         return iter(self.reactions.values())
@@ -508,7 +508,7 @@ class IncidentPhoton(EqualityMixin):
         # Get atomic number based on name of ACE table
         zaid, xs = ace.name.split('.')
         if not xs.endswith('p'):
-            raise TypeError("{} is not a photoatomic transport ACE table.".format(ace))
+            raise TypeError(f"{ace} is not a photoatomic transport ACE table.")
         Z = get_metadata(int(zaid))[2]
 
         # Read each reaction
@@ -638,7 +638,7 @@ class IncidentPhoton(EqualityMixin):
             with h5py.File(filename, 'r') as f:
                 _COMPTON_PROFILES['pz'] = f['pz'][()]
                 for i in range(1, 101):
-                    group = f['{:03}'.format(i)]
+                    group = f[f'{i:03}']
                     num_electrons = group['num_electrons'][()]
                     binding_energy = group['binding_energy'][()]*EV_PER_MEV
                     J = group['J'][()]
@@ -713,7 +713,7 @@ class IncidentPhoton(EqualityMixin):
 
         # Check for necessary reactions
         for mt in (502, 504, 522):
-            assert mt in data, "Reaction {} not found".format(mt)
+            assert mt in data, f"Reaction {mt} not found"
 
         # Read atomic relaxation
         data.atomic_relaxation = AtomicRelaxation.from_hdf5(group['subshells'])
@@ -764,7 +764,7 @@ class IncidentPhoton(EqualityMixin):
         """
         with h5py.File(str(path), mode, libver=libver) as f:
             # Write filetype and version
-            f.attrs['filetype'] = np.string_('data_photon')
+            f.attrs['filetype'] = np.bytes_('data_photon')
             if 'version' not in f.attrs:
                 f.attrs['version'] = np.array(HDF5_VERSION)
 
@@ -836,7 +836,7 @@ class IncidentPhoton(EqualityMixin):
             filename = os.path.join(os.path.dirname(__file__), 'density_effect.h5')
             with h5py.File(filename, 'r') as f:
                 for i in range(1, 101):
-                    group = f['{:03}'.format(i)]
+                    group = f[f'{i:03}']
                     _BREMSSTRAHLUNG[i] = {
                         'I': group.attrs['I'],
                         'num_electrons': group['num_electrons'][()],
@@ -924,10 +924,9 @@ class PhotonReaction(EqualityMixin):
 
     def __repr__(self):
         if self.mt in _REACTION_NAME:
-            return "<Photon Reaction: MT={} {}>".format(
-                self.mt, _REACTION_NAME[self.mt][0])
+            return f"<Photon Reaction: MT={self.mt} {_REACTION_NAME[self.mt][0]}>"
         else:
-            return "<Photon Reaction: MT={}>".format(self.mt)
+            return f"<Photon Reaction: MT={self.mt}>"
 
     @property
     def anomalous_real(self):

openmc/data/product.py

@@ -124,8 +124,8 @@ class Product(EqualityMixin):
             HDF5 group to write to
 
         """
-        group.attrs['particle'] = np.string_(self.particle)
-        group.attrs['emission_mode'] = np.string_(self.emission_mode)
+        group.attrs['particle'] = np.bytes_(self.particle)
+        group.attrs['emission_mode'] = np.bytes_(self.emission_mode)
         if self.decay_rate > 0.0:
             group.attrs['decay_rate'] = self.decay_rate
 
@@ -135,7 +135,7 @@ class Product(EqualityMixin):
         # Write applicability/distribution
         group.attrs['n_distribution'] = len(self.distribution)
         for i, d in enumerate(self.distribution):
-            dgroup = group.create_group('distribution_{}'.format(i))
+            dgroup = group.create_group(f'distribution_{i}')
             if self.applicability:
                 self.applicability[i].to_hdf5(dgroup, 'applicability')
             d.to_hdf5(dgroup)
@@ -170,7 +170,7 @@ class Product(EqualityMixin):
         distribution = []
         applicability = []
         for i in range(n_distribution):
-            dgroup = group['distribution_{}'.format(i)]
+            dgroup = group[f'distribution_{i}']
             if 'applicability' in dgroup:
                 applicability.append(Tabulated1D.from_hdf5(
                     dgroup['applicability']))

openmc/data/reaction.py

@@ -56,13 +56,13 @@ REACTION_NAME = {1: '(n,total)', 2: '(n,elastic)', 4: '(n,level)',
                  301: 'heating', 444: 'damage-energy',
                  649: '(n,pc)', 699: '(n,dc)', 749: '(n,tc)', 799: '(n,3Hec)',
                  849: '(n,ac)', 891: '(n,2nc)', 901: 'heating-local'}
-REACTION_NAME.update({i: '(n,n{})'.format(i - 50) for i in range(51, 91)})
-REACTION_NAME.update({i: '(n,p{})'.format(i - 600) for i in range(600, 649)})
-REACTION_NAME.update({i: '(n,d{})'.format(i - 650) for i in range(650, 699)})
-REACTION_NAME.update({i: '(n,t{})'.format(i - 700) for i in range(700, 749)})
-REACTION_NAME.update({i: '(n,3He{})'.format(i - 750) for i in range(750, 799)})
-REACTION_NAME.update({i: '(n,a{})'.format(i - 800) for i in range(800, 849)})
-REACTION_NAME.update({i: '(n,2n{})'.format(i - 875) for i in range(875, 891)})
+REACTION_NAME.update({i: f'(n,n{i - 50})' for i in range(51, 91)})
+REACTION_NAME.update({i: f'(n,p{i - 600})' for i in range(600, 649)})
+REACTION_NAME.update({i: f'(n,d{i - 650})' for i in range(650, 699)})
+REACTION_NAME.update({i: f'(n,t{i - 700})' for i in range(700, 749)})
+REACTION_NAME.update({i: f'(n,3He{i - 750})' for i in range(750, 799)})
+REACTION_NAME.update({i: f'(n,a{i - 800})' for i in range(800, 849)})
+REACTION_NAME.update({i: f'(n,2n{i - 875})' for i in range(875, 891)})
 
 REACTION_MT = {name: mt for mt, name in REACTION_NAME.items()}
 REACTION_MT['fission'] = 18
@@ -119,7 +119,7 @@ def _get_products(ev, mt):
             p = Product('electron')
         else:
             Z, A = divmod(za, 1000)
-            p = Product('{}{}'.format(ATOMIC_SYMBOL[Z], A))
+            p = Product(f'{ATOMIC_SYMBOL[Z]}{A}')
 
         p.yield_ = yield_
 
@@ -557,9 +557,9 @@ def _get_activation_products(ev, rx):
             # Get GNDS name for product
             symbol = ATOMIC_SYMBOL[Z]
             if excited_state > 0:
-                name = '{}{}_e{}'.format(symbol, A, excited_state)
+                name = f'{symbol}{A}_e{excited_state}'
             else:
-                name = '{}{}'.format(symbol, A)
+                name = f'{symbol}{A}'
 
             p = Product(name)
             if mf == 9:
@@ -656,8 +656,7 @@ def _get_photon_products_ace(ace, rx):
             photon.yield_ = Tabulated1D(energy, yield_)
 
         else:
-            raise ValueError("MFTYPE must be 12, 13, 16. Got {0}".format(
-                mftype))
+            raise ValueError(f"MFTYPE must be 12, 13, 16. Got {mftype}")
 
         # ==================================================================
         # Photon energy distribution
@@ -846,9 +845,9 @@ class Reaction(EqualityMixin):
 
     def __repr__(self):
         if self.mt in REACTION_NAME:
-            return "<Reaction: MT={} {}>".format(self.mt, REACTION_NAME[self.mt])
+            return f"<Reaction: MT={self.mt} {REACTION_NAME[self.mt]}>"
         else:
-            return "<Reaction: MT={}>".format(self.mt)
+            return f"<Reaction: MT={self.mt}>"
 
     @property
     def center_of_mass(self):
@@ -920,9 +919,9 @@ class Reaction(EqualityMixin):
 
         group.attrs['mt'] = self.mt
         if self.mt in REACTION_NAME:
-            group.attrs['label'] = np.string_(REACTION_NAME[self.mt])
+            group.attrs['label'] = np.bytes_(REACTION_NAME[self.mt])
         else:
-            group.attrs['label'] = np.string_(self.mt)
+            group.attrs['label'] = np.bytes_(self.mt)
         group.attrs['Q_value'] = self.q_value
         group.attrs['center_of_mass'] = 1 if self.center_of_mass else 0
         group.attrs['redundant'] = 1 if self.redundant else 0
@@ -933,7 +932,7 @@ class Reaction(EqualityMixin):
                 threshold_idx = getattr(self.xs[T], '_threshold_idx', 0)
                 dset.attrs['threshold_idx'] = threshold_idx
         for i, p in enumerate(self.products):
-            pgroup = group.create_group('product_{}'.format(i))
+            pgroup = group.create_group(f'product_{i}')
             p.to_hdf5(pgroup)
 
     @classmethod
@@ -985,7 +984,7 @@ class Reaction(EqualityMixin):
 
         # Read reaction products
         for i in range(n_product):
-            pgroup = group['product_{}'.format(i)]
+            pgroup = group[f'product_{i}']
             rx.products.append(Product.from_hdf5(pgroup))
 
         return rx

openmc/data/resonance.py

@@ -830,7 +830,7 @@ class RMatrixLimited(ResonanceRange):
                 elif mt == 102:
                     columns.append('captureWidth')
                 else:
-                    columns.append('width (MT={})'.format(mt))
+                    columns.append(f'width (MT={mt})')
 
             # Create Pandas dataframe with resonance parameters
             parameters = pd.DataFrame.from_records(records, columns=columns)
@@ -896,7 +896,7 @@ class SpinGroup:
         self.parameters = parameters
 
     def __repr__(self):
-        return '<SpinGroup: Jpi={}{}>'.format(self.spin, self.parity)
+        return f'<SpinGroup: Jpi={self.spin}{self.parity}>'
 
 
 class Unresolved(ResonanceRange):

openmc/data/resonance_covariance.py

@@ -239,14 +239,14 @@ class ResonanceCovarianceRange:
         samples = []
 
         # Handling MLBW/SLBW sampling
+        rng = np.random.default_rng()
         if formalism == 'mlbw' or formalism == 'slbw':
             params = ['energy', 'neutronWidth', 'captureWidth', 'fissionWidth',
                       'competitiveWidth']
             param_list = params[:mpar]
             mean_array = parameters[param_list].values
             mean = mean_array.flatten()
-            par_samples = np.random.multivariate_normal(mean, cov,
-                                                        size=n_samples)
+            par_samples = rng.multivariate_normal(mean, cov, size=n_samples)
             spin = parameters['J'].values
             l_value = parameters['L'].values
             for sample in par_samples:
@@ -277,8 +277,7 @@ class ResonanceCovarianceRange:
             param_list = params[:mpar]
             mean_array = parameters[param_list].values
             mean = mean_array.flatten()
-            par_samples = np.random.multivariate_normal(mean, cov,
-                                                        size=n_samples)
+            par_samples = rng.multivariate_normal(mean, cov, size=n_samples)
             spin = parameters['J'].values
             l_value = parameters['L'].values
             for sample in par_samples:

openmc/data/thermal.py

@@ -90,7 +90,7 @@ _THERMAL_NAMES = {
 def _temperature_str(T):
     # round() normally returns an int when called with a single argument, but
     # numpy floats overload rounding to return another float
-    return "{}K".format(int(round(T)))
+    return f"{int(round(T))}K"
 
 
 def get_thermal_name(name):
@@ -221,7 +221,7 @@ class CoherentElastic(Function1D):
         """
         dataset = group.create_dataset(name, data=np.vstack(
             [self.bragg_edges, self.factors]))
-        dataset.attrs['type'] = np.string_(type(self).__name__)
+        dataset.attrs['type'] = np.bytes_(type(self).__name__)
 
     @classmethod
     def from_hdf5(cls, dataset):
@@ -294,7 +294,7 @@ class IncoherentElastic(Function1D):
         """
         data = np.array([self.bound_xs, self.debye_waller])
         dataset = group.create_dataset(name, data=data)
-        dataset.attrs['type'] = np.string_(type(self).__name__)
+        dataset.attrs['type'] = np.bytes_(type(self).__name__)
 
     @classmethod
     def from_hdf5(cls, dataset):
@@ -439,7 +439,7 @@ class ThermalScattering(EqualityMixin):
 
     def __repr__(self):
         if hasattr(self, 'name'):
-            return "<Thermal Scattering Data: {}>".format(self.name)
+            return f"<Thermal Scattering Data: {self.name}>"
         else:
             return "<Thermal Scattering Data>"
 
@@ -464,7 +464,7 @@ class ThermalScattering(EqualityMixin):
         """
         # Open file and write version
         with h5py.File(str(path), mode, libver=libver) as f:
-            f.attrs['filetype'] = np.string_('data_thermal')
+            f.attrs['filetype'] = np.bytes_('data_thermal')
             f.attrs['version'] = np.array(HDF5_VERSION)
 
             # Write basic data
@@ -506,7 +506,7 @@ class ThermalScattering(EqualityMixin):
         # Check if temprature already exists
         strT = data.temperatures[0]
         if strT in self.temperatures:
-            warn('S(a,b) data at T={} already exists.'.format(strT))
+            warn(f'S(a,b) data at T={strT} already exists.')
             return
 
         # Check that name matches
@@ -614,7 +614,7 @@ class ThermalScattering(EqualityMixin):
         # Get new name that is GND-consistent
         ace_name, xs = ace.name.split('.')
         if not xs.endswith('t'):
-            raise TypeError("{} is not a thermal scattering ACE table.".format(ace))
+            raise TypeError(f"{ace} is not a thermal scattering ACE table.")
         if name is None:
             name = get_thermal_name(ace_name)
 
@@ -698,8 +698,8 @@ class ThermalScattering(EqualityMixin):
                 # add an outgoing energy 0 eV that has a PDF of 0 (and of
                 # course, a CDF of 0 as well).
                 if eout_i.c[0] > 0.:
-                    eout_i.x = np.insert(eout_i.x, 0, 0.)
-                    eout_i.p = np.insert(eout_i.p, 0, 0.)
+                    eout_i._x = np.insert(eout_i.x, 0, 0.)
+                    eout_i._p = np.insert(eout_i.p, 0, 0.)
                     eout_i.c = np.insert(eout_i.c, 0, 0.)
 
                     # For this added outgoing energy (of 0 eV) we add a set of

openmc/data/thermal_angle_energy.py

@@ -43,7 +43,7 @@ class CoherentElasticAE(AngleEnergy):
             HDF5 group to write to
 
         """
-        group.attrs['type'] = np.string_('coherent_elastic')
+        group.attrs['type'] = np.bytes_('coherent_elastic')
         self.coherent_xs.to_hdf5(group, 'coherent_xs')
 
     @classmethod
@@ -104,7 +104,7 @@ class IncoherentElasticAE(AngleEnergy):
             HDF5 group to write to
 
         """
-        group.attrs['type'] = np.string_('incoherent_elastic')
+        group.attrs['type'] = np.bytes_('incoherent_elastic')
         group.create_dataset('debye_waller', data=self.debye_waller)
 
     @classmethod
@@ -146,7 +146,7 @@ class IncoherentElasticAEDiscrete(AngleEnergy):
             HDF5 group to write to
 
         """
-        group.attrs['type'] = np.string_('incoherent_elastic_discrete')
+        group.attrs['type'] = np.bytes_('incoherent_elastic_discrete')
         group.create_dataset('mu_out', data=self.mu_out)
 
     @classmethod
@@ -203,7 +203,7 @@ class IncoherentInelasticAEDiscrete(AngleEnergy):
             HDF5 group to write to
 
         """
-        group.attrs['type'] = np.string_('incoherent_inelastic_discrete')
+        group.attrs['type'] = np.bytes_('incoherent_inelastic_discrete')
         group.create_dataset('energy_out', data=self.energy_out)
         group.create_dataset('mu_out', data=self.mu_out)
         group.create_dataset('skewed', data=self.skewed)
@@ -266,7 +266,7 @@ class MixedElasticAE(AngleEnergy):
             HDF5 group to write to
 
         """
-        group.attrs['type'] = np.string_('mixed_elastic')
+        group.attrs['type'] = np.bytes_('mixed_elastic')
         coherent_group = group.create_group('coherent')
         self.coherent.to_hdf5(coherent_group)
         incoherent_group = group.create_group('incoherent')

openmc/data/uncorrelated.py

@@ -63,7 +63,7 @@ class UncorrelatedAngleEnergy(AngleEnergy):
             HDF5 group to write to
 
         """
-        group.attrs['type'] = np.string_('uncorrelated')
+        group.attrs['type'] = np.bytes_('uncorrelated')
         if self.angle is not None:
             angle_group = group.create_group('angle')
             self.angle.to_hdf5(angle_group)

openmc/deplete/abc.py

@@ -10,8 +10,6 @@ from collections.abc import Iterable, Callable
 from copy import deepcopy
 from inspect import signature
 from numbers import Real, Integral
-from contextlib import contextmanager
-import os
 from pathlib import Path
 import time
 from typing import Optional, Union, Sequence
@@ -22,6 +20,7 @@ from uncertainties import ufloat
 
 from openmc.checkvalue import check_value, check_type, check_greater_than, PathLike
 from openmc.mpi import comm
+from openmc.utility_funcs import change_directory
 from openmc import Material
 from .stepresult import StepResult
 from .chain import Chain
@@ -66,25 +65,6 @@ except AttributeError:
     # Can't set __doc__ on properties on Python 3.4
     pass
 
-@contextmanager
-def change_directory(output_dir):
-    """
-    Helper function for managing the current directory.
-
-    Parameters
-    ----------
-    output_dir : pathlib.Path
-        Directory to switch to.
-    """
-    orig_dir  = os.getcwd()
-    output_dir = Path(output_dir)
-    output_dir.mkdir(parents=True, exist_ok=True)
-    try:
-        os.chdir(output_dir)
-        yield
-    finally:
-        os.chdir(orig_dir)
-
 
 class TransportOperator(ABC):
     """Abstract class defining a transport operator
@@ -655,7 +635,7 @@ class Integrator(ABC):
                 days = watt_days_per_kg * kilograms / rate
                 seconds.append(days*_SECONDS_PER_DAY)
             else:
-                raise ValueError("Invalid timestep unit '{}'".format(unit))
+                raise ValueError(f"Invalid timestep unit '{unit}'")
 
         self.timesteps = np.asarray(seconds)
         self.source_rates = np.asarray(source_rates)
@@ -673,8 +653,7 @@ class Integrator(ABC):
                 self._solver = CRAM16
             else:
                 raise ValueError(
-                    "Solver {} not understood. Expected 'cram48' or "
-                    "'cram16'".format(solver))
+                    f"Solver {solver} not understood. Expected 'cram48' or 'cram16'")
         else:
             self.solver = solver
 
@@ -686,14 +665,13 @@ class Integrator(ABC):
     def solver(self, func):
         if not isinstance(func, Callable):
             raise TypeError(
-                "Solver must be callable, not {}".format(type(func)))
+                f"Solver must be callable, not {type(func)}")
         try:
             sig = signature(func)
         except ValueError:
             # Guard against callables that aren't introspectable, e.g.
             # fortran functions wrapped by F2PY
-            warn("Could not determine arguments to {}. Proceeding "
-                 "anyways".format(func))
+            warn(f"Could not determine arguments to {func}. Proceeding anyways")
             self._solver = func
             return
 
@@ -705,8 +683,7 @@ class Integrator(ABC):
         for ix, param in enumerate(sig.parameters.values()):
             if param.kind in {param.KEYWORD_ONLY, param.VAR_KEYWORD}:
                 raise ValueError(
-                    "Keyword arguments like {} at position {} are not "
-                    "allowed".format(ix, param))
+                    f"Keyword arguments like {ix} at position {param} are not allowed")
 
         self._solver = func
 

openmc/deplete/chain.py

@@ -141,7 +141,7 @@ def replace_missing(product, decay_data):
 
     # First check if ground state is available
     if state:
-        product = '{}{}'.format(symbol, A)
+        product = f'{symbol}{A}'
 
     # Find isotope with longest half-life
     half_life = 0.0
@@ -172,7 +172,7 @@ def replace_missing(product, decay_data):
                 Z += 1
             else:
                 Z -= 1
-        product = '{}{}'.format(openmc.data.ATOMIC_SYMBOL[Z], A)
+        product = f'{openmc.data.ATOMIC_SYMBOL[Z]}{A}'
 
     return product
 
@@ -417,7 +417,7 @@ class Chain:
                     if mts & reactions_available:
                         A = data.nuclide['mass_number'] + delta_A
                         Z = data.nuclide['atomic_number'] + delta_Z
-                        daughter = '{}{}'.format(openmc.data.ATOMIC_SYMBOL[Z], A)
+                        daughter = f'{openmc.data.ATOMIC_SYMBOL[Z]}{A}'
 
                         if daughter not in decay_data:
                             daughter = replace_missing(daughter, decay_data)
@@ -483,7 +483,7 @@ class Chain:
         if missing_daughter:
             print('The following decay modes have daughters with no decay data:')
             for mode in missing_daughter:
-                print('  {}'.format(mode))
+                print(f'  {mode}')
             print('')
 
         if missing_rx_product:
@@ -495,7 +495,7 @@ class Chain:
         if missing_fpy:
             print('The following fissionable nuclides have no fission product yields:')
             for parent, replacement in missing_fpy:
-                print('  {}, replaced with {}'.format(parent, replacement))
+                print(f'  {parent}, replaced with {replacement}')
             print('')
 
         if missing_fp:
@@ -597,9 +597,12 @@ class Chain:
         --------
         :meth:`get_default_fission_yields`
         """
-        matrix = defaultdict(float)
         reactions = set()
 
+        # Use DOK matrix as intermediate representation for matrix
+        n = len(self)
+        matrix = sp.dok_matrix((n, n))
+
         if fission_yields is None:
             fission_yields = self.get_default_fission_yields()
 
@@ -674,11 +677,8 @@ class Chain:
                 # Clear set of reactions
                 reactions.clear()
 
-        # Use DOK matrix as intermediate representation, then convert to CSC and return
-        n = len(self)
-        matrix_dok = sp.dok_matrix((n, n))
-        dict.update(matrix_dok, matrix)
-        return matrix_dok.tocsc()
+        # Return CSC representation instead of DOK
+        return matrix.tocsc()
 
     def form_rr_term(self, tr_rates, mats):
         """Function to form the transfer rate term matrices.
@@ -711,7 +711,9 @@ class Chain:
             Sparse matrix representing transfer term.
 
         """
-        matrix = defaultdict(float)
+        # Use DOK as intermediate representation
+        n = len(self)
+        matrix = sp.dok_matrix((n, n))
 
         for i, nuc in enumerate(self.nuclides):
             elm = re.split(r'\d+', nuc.name)[0]
@@ -737,10 +739,9 @@ class Chain:
                 else:
                     matrix[i, i] = 0.0
             #Nothing else is allowed
-        n = len(self)
-        matrix_dok = sp.dok_matrix((n, n))
-        dict.update(matrix_dok, matrix)
-        return matrix_dok.tocsc()
+
+        # Return CSC instead of DOK
+        return matrix.tocsc()
 
     def get_branch_ratios(self, reaction="(n,gamma)"):
         """Return a dictionary with reaction branching ratios
@@ -872,8 +873,7 @@ class Chain:
             if len(indexes) == 0:
                 if strict:
                     raise AttributeError(
-                        "Nuclide {} does not have {} reactions".format(
-                            parent, reaction))
+                        f"Nuclide {parent} does not have {reaction} reactions")
                 missing_reaction.add(parent)
                 continue
 
@@ -895,8 +895,7 @@ class Chain:
 
         if len(rxn_ix_map) == 0:
             raise IndexError(
-                "No {} reactions found in this {}".format(
-                    reaction, self.__class__.__name__))
+                f"No {reaction} reactions found in this {self.__class__.__name__}")
 
         if len(missing_parents) > 0:
             warn("The following nuclides were not found in {}: {}".format(
@@ -907,14 +906,14 @@ class Chain:
                  "{}".format(reaction, ", ".join(sorted(missing_reaction))))
 
         if len(missing_products) > 0:
-            tail = ("{} -> {}".format(k, v)
+            tail = (f"{k} -> {v}"
                     for k, v in sorted(missing_products.items()))
             warn("The following products were not found in the {} and "
                  "parents were unmodified: \n{}".format(
                      self.__class__.__name__, ", ".join(tail)))
 
         if len(bad_sums) > 0:
-            tail = ("{}: {:5.3f}".format(k, s)
+            tail = (f"{k}: {s:5.3f}"
                     for k, s in sorted(bad_sums.items()))
             warn("The following parent nuclides were given {} branch ratios "
                  "with a sum outside tolerance of 1 +/- {:5.3e}:\n{}".format(

openmc/deplete/coupled_operator.py

@@ -291,7 +291,7 @@ class CoupledOperator(OpenMCOperator):
         # on this process
         if comm.size != 1:
             prev_results = self.prev_res
-            self.prev_res = Results()
+            self.prev_res = Results(filename=None)
             mat_indexes = _distribute(range(len(self.burnable_mats)))
             for res_obj in prev_results:
                 new_res = res_obj.distribute(self.local_mats, mat_indexes)
@@ -518,7 +518,7 @@ class CoupledOperator(OpenMCOperator):
 
         """
         openmc.lib.statepoint_write(
-            "openmc_simulation_n{}.h5".format(step),
+            f"openmc_simulation_n{step}.h5",
             write_source=False)
 
     def finalize(self):

openmc/deplete/independent_operator.py

@@ -130,7 +130,6 @@ class IndependentOperator(OpenMCOperator):
         # Validate micro-xs parameters
         check_type('materials', materials, openmc.Materials)
         check_type('micros', micros, Iterable, MicroXS)
-        check_type('fluxes', fluxes, Iterable, float)
 
         if not (len(fluxes) == len(micros) == len(materials)):
             msg = (f'The length of fluxes ({len(fluxes)}) should be equal to '
@@ -241,7 +240,6 @@ class IndependentOperator(OpenMCOperator):
                    reduce_chain_level=reduce_chain_level,
                    fission_yield_opts=fission_yield_opts)
 
-
     @staticmethod
     def _consolidate_nuclides_to_material(nuclides, nuc_units, volume):
         """Puts nuclide list into an openmc.Materials object.
@@ -270,7 +268,6 @@ class IndependentOperator(OpenMCOperator):
         self.prev_res[-1].transfer_volumes(model)
         self.materials = model.materials
 
-
         # Store previous results in operator
         # Distribute reaction rates according to those tracked
         # on this process
@@ -282,7 +279,6 @@ class IndependentOperator(OpenMCOperator):
                 new_res = res_obj.distribute(self.local_mats, mat_indexes)
                 self.prev_res.append(new_res)
 
-
     def _get_nuclides_with_data(self, cross_sections: List[MicroXS]) -> Set[str]:
         """Finds nuclides with cross section data"""
         return set(cross_sections[0].nuclides)
@@ -412,6 +408,12 @@ class IndependentOperator(OpenMCOperator):
 
         self._update_materials_and_nuclides(vec)
 
+        # If the source rate is zero, return zero reaction rates
+        if source_rate == 0.0:
+            rates = self.reaction_rates.copy()
+            rates.fill(0.0)
+            return OperatorResult(ufloat(0.0, 0.0), rates)
+
         rates = self._calculate_reaction_rates(source_rate)
         keff = self._keff
 

openmc/deplete/microxs.py

@@ -13,24 +13,24 @@ import numpy as np
 
 from openmc.checkvalue import check_type, check_value, check_iterable_type, PathLike
 from openmc.exceptions import DataError
+from openmc.utility_funcs import change_directory
 from openmc import StatePoint
 from openmc.mgxs import GROUP_STRUCTURES
 from openmc.data import REACTION_MT
 import openmc
-from .abc import change_directory
 from .chain import Chain, REACTIONS
 from .coupled_operator import _find_cross_sections, _get_nuclides_with_data
 import openmc.lib
 
 _valid_rxns = list(REACTIONS)
 _valid_rxns.append('fission')
+_valid_rxns.append('damage-energy')
 
 
 def _resolve_chain_file_path(chain_file: str):
-    # Determine what reactions and nuclides are available in chain
     if chain_file is None:
         chain_file = openmc.config.get('chain_file')
-        if 'chain_file' in openmc.config:
+        if 'chain_file' not in openmc.config:
             raise DataError(
                 "No depletion chain specified and could not find depletion "
                 "chain in openmc.config['chain_file']"
@@ -284,38 +284,37 @@ class MicroXS:
         # Create 2D array for microscopic cross sections
         microxs_arr = np.zeros((len(nuclides), len(mts)))
 
-        with TemporaryDirectory() as tmpdir:
-            # Create a material with all nuclides
-            mat_all_nucs = openmc.Material()
-            for nuc in nuclides:
-                if nuc in nuclides_with_data:
-                    mat_all_nucs.add_nuclide(nuc, 1.0)
-            mat_all_nucs.set_density("atom/b-cm", 1.0)
+        # Create a material with all nuclides
+        mat_all_nucs = openmc.Material()
+        for nuc in nuclides:
+            if nuc in nuclides_with_data:
+                mat_all_nucs.add_nuclide(nuc, 1.0)
+        mat_all_nucs.set_density("atom/b-cm", 1.0)
 
-            # Create simple model containing the above material
-            surf1 = openmc.Sphere(boundary_type="vacuum")
-            surf1_cell = openmc.Cell(fill=mat_all_nucs, region=-surf1)
-            model = openmc.Model()
-            model.geometry = openmc.Geometry([surf1_cell])
-            model.settings = openmc.Settings(
-                particles=1, batches=1, output={'summary': False})
+        # Create simple model containing the above material
+        surf1 = openmc.Sphere(boundary_type="vacuum")
+        surf1_cell = openmc.Cell(fill=mat_all_nucs, region=-surf1)
+        model = openmc.Model()
+        model.geometry = openmc.Geometry([surf1_cell])
+        model.settings = openmc.Settings(
+            particles=1, batches=1, output={'summary': False})
 
-            with change_directory(tmpdir):
-                # Export model within temporary directory
-                model.export_to_model_xml()
+        with change_directory(tmpdir=True):
+            # Export model within temporary directory
+            model.export_to_model_xml()
 
-                with openmc.lib.run_in_memory(**init_kwargs):
-                    # For each nuclide and reaction, compute the flux-averaged
-                    # cross section
-                    for nuc_index, nuc in enumerate(nuclides):
-                        if nuc not in nuclides_with_data:
-                            continue
-                        lib_nuc = openmc.lib.nuclides[nuc]
-                        for mt_index, mt in enumerate(mts):
-                            xs = lib_nuc.collapse_rate(
-                                mt, temperature, energies, multigroup_flux
-                            )
-                            microxs_arr[nuc_index, mt_index] = xs
+            with openmc.lib.run_in_memory(**init_kwargs):
+                # For each nuclide and reaction, compute the flux-averaged
+                # cross section
+                for nuc_index, nuc in enumerate(nuclides):
+                    if nuc not in nuclides_with_data:
+                        continue
+                    lib_nuc = openmc.lib.nuclides[nuc]
+                    for mt_index, mt in enumerate(mts):
+                        xs = lib_nuc.collapse_rate(
+                            mt, temperature, energies, multigroup_flux
+                        )
+                        microxs_arr[nuc_index, mt_index] = xs
 
         return cls(microxs_arr, nuclides, reactions)
 

openmc/deplete/nuclide.py

@@ -275,7 +275,7 @@ class Nuclide:
             if parent is not None:
                 assert root is not None
                 fpy_elem = root.find(
-                    './/nuclide[@name="{}"]/neutron_fission_yields'.format(parent)
+                    f'.//nuclide[@name="{parent}"]/neutron_fission_yields'
                 )
                 if fpy_elem is None:
                     raise ValueError(
@@ -413,7 +413,7 @@ class Nuclide:
                     continue
                 msg = msg_func(
                     name=self.name, actual=sum_br, expected=1.0, tol=tolerance,
-                    prop="{} reaction branch ratios".format(rxn_type))
+                    prop=f"{rxn_type} reaction branch ratios")
                 if strict:
                     raise ValueError(msg)
                 elif quiet:
@@ -430,7 +430,7 @@ class Nuclide:
                 msg = msg_func(
                     name=self.name, actual=sum_yield,
                     expected=2.0, tol=tolerance,
-                    prop="fission yields (E = {:7.4e} eV)".format(energy))
+                    prop=f"fission yields (E = {energy:7.4e} eV)")
                 if strict:
                     raise ValueError(msg)
                 elif quiet:
@@ -695,8 +695,7 @@ class FissionYield(Mapping):
         return self * scalar
 
     def __repr__(self):
-        return "<{} containing {} products and yields>".format(
-            self.__class__.__name__, len(self))
+        return f"<{self.__class__.__name__} containing {len(self)} products and yields>"
 
     def __deepcopy__(self, memo):
         result = FissionYield(self.products, self.yields.copy())

openmc/deplete/stepresult.py

@@ -294,7 +294,7 @@ class StepResult:
         # Store concentration mat and nuclide dictionaries (along with volumes)
 
         handle.attrs['version'] = np.array(VERSION_RESULTS)
-        handle.attrs['filetype'] = np.string_('depletion results')
+        handle.attrs['filetype'] = np.bytes_('depletion results')
 
         mat_list = sorted(self.mat_to_hdf5_ind, key=int)
         nuc_list = sorted(self.index_nuc)

openmc/element.py

@@ -1,4 +1,5 @@
 import re
+import warnings
 
 import lxml.etree as ET
 
@@ -123,6 +124,10 @@ class Element(str):
         # Get the nuclides present in nature
         natural_nuclides = {name for name, abundance in natural_isotopes(self)}
 
+        # Issue warning if no existing nuclides
+        if len(natural_nuclides) == 0:
+            warnings.warn(f"No naturally occurring isotopes found for {self}.")
+
         # Create dict to store the expanded nuclides and abundances
         abundances = {}
 

openmc/examples.py

@@ -76,8 +76,10 @@ def pwr_pin_cell():
     model.settings.batches = 10
     model.settings.inactive = 5
     model.settings.particles = 100
-    model.settings.source = openmc.IndependentSource(space=openmc.stats.Box(
-        [-pitch/2, -pitch/2, -1], [pitch/2, pitch/2, 1], only_fissionable=True))
+    model.settings.source = openmc.IndependentSource(
+        space=openmc.stats.Box([-pitch/2, -pitch/2, -1], [pitch/2, pitch/2, 1]),
+        constraints={'fissionable': True}
+    )
 
     plot = openmc.Plot.from_geometry(model.geometry)
     plot.pixels = (300, 300)
@@ -527,8 +529,10 @@ def pwr_assembly():
     model.settings.batches = 10
     model.settings.inactive = 5
     model.settings.particles = 100
-    model.settings.source = openmc.IndependentSource(space=openmc.stats.Box(
-        [-pitch/2, -pitch/2, -1], [pitch/2, pitch/2, 1], only_fissionable=True))
+    model.settings.source = openmc.IndependentSource(
+        space=openmc.stats.Box([-pitch/2, -pitch/2, -1], [pitch/2, pitch/2, 1]),
+        constraints={'fissionable': True}
+    )
 
     plot = openmc.Plot()
     plot.origin = (0.0, 0.0, 0)

openmc/filter.py

@@ -1,3 +1,4 @@
+from __future__ import annotations
 from abc import ABCMeta
 from collections.abc import Iterable
 import hashlib
@@ -804,8 +805,8 @@ class MeshFilter(Filter):
     id : int
         Unique identifier for the filter
     translation : Iterable of float
-        This array specifies a vector that is used to translate (shift)
-        the mesh for this filter
+        This array specifies a vector that is used to translate (shift) the mesh
+        for this filter
     bins : list of tuple
         A list of mesh indices for each filter bin, e.g. [(1, 1, 1), (2, 1, 1),
         ...]
@@ -846,7 +847,6 @@ class MeshFilter(Filter):
         mesh_obj = kwargs['meshes'][mesh_id]
         filter_id = int(group.name.split('/')[-1].lstrip('filter '))
 
-
         out = cls(mesh_obj, filter_id=filter_id)
 
         translation = group.get('translation')
@@ -864,10 +864,10 @@ class MeshFilter(Filter):
         cv.check_type('filter mesh', mesh, openmc.MeshBase)
         self._mesh = mesh
         if isinstance(mesh, openmc.UnstructuredMesh):
-            if mesh.volumes is None:
-                self.bins = []
-            else:
+            if mesh.has_statepoint_data:
                 self.bins = list(range(len(mesh.volumes)))
+            else:
+                self.bins = []
         else:
             self.bins = list(mesh.indices)
 
@@ -972,7 +972,7 @@ class MeshFilter(Filter):
         return element
 
     @classmethod
-    def from_xml_element(cls, elem, **kwargs):
+    def from_xml_element(cls, elem: ET.Element, **kwargs) -> MeshFilter:
         mesh_id = int(get_text(elem, 'bins'))
         mesh_obj = kwargs['meshes'][mesh_id]
         filter_id = int(elem.get('id'))
@@ -984,6 +984,34 @@ class MeshFilter(Filter):
         return out
 
 
+class MeshBornFilter(MeshFilter):
+    """Filter events by the mesh cell a particle originated from.
+
+    Parameters
+    ----------
+    mesh : openmc.MeshBase
+        The mesh object that events will be tallied onto
+    filter_id : int
+        Unique identifier for the filter
+
+    Attributes
+    ----------
+    mesh : openmc.MeshBase
+        The mesh object that events will be tallied onto
+    id : int
+        Unique identifier for the filter
+    translation : Iterable of float
+        This array specifies a vector that is used to translate (shift)
+        the mesh for this filter
+    bins : list of tuple
+        A list of mesh indices for each filter bin, e.g. [(1, 1, 1), (2, 1, 1),
+        ...]
+    num_bins : Integral
+        The number of filter bins
+
+    """
+
+
 class MeshSurfaceFilter(MeshFilter):
     """Filter events by surface crossings on a mesh.
 
@@ -1349,6 +1377,10 @@ class EnergyFilter(RealFilter):
     """
     units = 'eV'
 
+    def __init__(self, values, filter_id=None):
+        cv.check_length('values', values, 2)
+        super().__init__(values, filter_id)
+
     def get_bin_index(self, filter_bin):
         # Use lower energy bound to find index for RealFilters
         deltas = np.abs(self.bins[:, 1] - filter_bin[1]) / filter_bin[1]

openmc/geometry.py

@@ -134,7 +134,7 @@ class Geometry:
 
         # Create XML representation
         element = ET.Element("geometry")
-        self.root_universe.create_xml_subelement(element, memo=set())
+        self.root_universe.create_xml_subelement(element)
 
         # Sort the elements in the file
         element[:] = sorted(element, key=lambda x: (
@@ -373,7 +373,7 @@ class Geometry:
 
         """
         if self.root_universe is not None:
-            return self.root_universe.get_all_cells(memo=set())
+            return self.root_universe.get_all_cells()
         else:
             return {}
 
@@ -417,7 +417,7 @@ class Geometry:
 
         """
         if self.root_universe is not None:
-            return self.root_universe.get_all_materials(memo=set())
+            return self.root_universe.get_all_materials()
         else:
             return {}
 
@@ -701,7 +701,7 @@ class Geometry:
             coeffs = tuple(round(surf._coefficients[k],
                                  self.surface_precision)
                            for k in surf._coeff_keys)
-            key = (surf._type,) + coeffs
+            key = (surf._type, surf._boundary_type) + coeffs
             redundancies[key].append(surf)
 
         redundant_surfaces = {replace.id: keep

openmc/lattice.py

@@ -170,11 +170,11 @@ class Lattice(IDManagerMixin, ABC):
         """
         cells = {}
 
-        if memo and self in memo:
+        if memo is None:
+            memo = set()
+        elif self in memo:
             return cells
-
-        if memo is not None:
-            memo.add(self)
+        memo.add(self)
 
         unique_universes = self.get_unique_universes()
 
@@ -194,6 +194,9 @@ class Lattice(IDManagerMixin, ABC):
 
         """
 
+        if memo is None:
+            memo = set()
+
         materials = {}
 
         # Append all Cells in each Cell in the Universe to the dictionary
@@ -203,7 +206,7 @@ class Lattice(IDManagerMixin, ABC):
 
         return materials
 
-    def get_all_universes(self):
+    def get_all_universes(self, memo=None):
         """Return all universes that are contained within the lattice
 
         Returns
@@ -213,11 +216,16 @@ class Lattice(IDManagerMixin, ABC):
             :class:`Universe` instances
 
         """
-
         # Initialize a dictionary of all Universes contained by the Lattice
         # in each nested Universe level
         all_universes = {}
 
+        if memo is None:
+            memo = set()
+        elif self in memo:
+            return all_universes
+        memo.add(self)
+
         # Get all unique Universes contained in each of the lattice cells
         unique_universes = self.get_unique_universes()
 
@@ -226,7 +234,7 @@ class Lattice(IDManagerMixin, ABC):
 
         # Append all Universes containing each cell to the dictionary
         for universe in unique_universes.values():
-            all_universes.update(universe.get_all_universes())
+            all_universes.update(universe.get_all_universes(memo))
 
         return all_universes
 
@@ -846,10 +854,11 @@ class RectLattice(Lattice):
 
         """
         # If the element already contains the Lattice subelement, then return
-        if memo and self in memo:
+        if memo is None:
+            memo = set()
+        elif self in memo:
             return
-        if memo is not None:
-            memo.add(self)
+        memo.add(self)
 
         # Make sure universes have been assigned
         if self.universes is None:
@@ -1417,10 +1426,11 @@ class HexLattice(Lattice):
 
     def create_xml_subelement(self, xml_element, memo=None):
         # If this subelement has already been written, return
-        if memo and self in memo:
+        if memo is None:
+            memo = set()
+        elif self in memo:
             return
-        if memo is not None:
-            memo.add(self)
+        memo.add(self)
 
         lattice_subelement = ET.Element("hex_lattice")
         lattice_subelement.set("id", str(self._id))
@@ -1564,6 +1574,7 @@ class HexLattice(Lattice):
                             alpha -= 1
                             if not lat.is_valid_index((x, alpha, z)):
                                 # Reached the bottom
+                                j += 1
                                 break
                     j += 1
             else:
@@ -1583,6 +1594,7 @@ class HexLattice(Lattice):
 
                         # Check if we've reached the bottom
                         if y == -n_rings:
+                            j += 1
                             break
 
                         while not lat.is_valid_index((alpha, y, z)):
@@ -1851,7 +1863,7 @@ class HexLattice(Lattice):
         largest_index = 6*(num_rings - 1)
         n_digits_index = len(str(largest_index))
         n_digits_ring = len(str(num_rings - 1))
-        str_form = '({{:{}}},{{:{}}})'.format(n_digits_ring, n_digits_index)
+        str_form = f'({{:{n_digits_ring}}},{{:{n_digits_index}}})'
         pad = ' '*(n_digits_index + n_digits_ring + 3)
 
         # Initialize the list for each row.
@@ -1956,7 +1968,7 @@ class HexLattice(Lattice):
         largest_index = 6*(num_rings - 1)
         n_digits_index = len(str(largest_index))
         n_digits_ring = len(str(num_rings - 1))
-        str_form = '({{:{}}},{{:{}}})'.format(n_digits_ring, n_digits_index)
+        str_form = f'({{:{n_digits_ring}}},{{:{n_digits_index}}})'
         pad = ' '*(n_digits_index + n_digits_ring + 3)
 
         # Initialize the list for each row.

openmc/lib/__init__.py

@@ -28,7 +28,7 @@ else:
 if os.environ.get('READTHEDOCS', None) != 'True':
     # Open shared library
     _filename = pkg_resources.resource_filename(
-        __name__, 'libopenmc.{}'.format(_suffix))
+        __name__, f'libopenmc.{_suffix}')
     _dll = CDLL(_filename)
 else:
     # For documentation builds, we don't actually have the shared library
@@ -54,6 +54,9 @@ def _libmesh_enabled():
 def _mcpl_enabled():
     return c_bool.in_dll(_dll, "MCPL_ENABLED").value
 
+def _uwuw_enabled():
+    return c_bool.in_dll(_dll, "UWUW_ENABLED").value
+
 
 from .error import *
 from .core import *

openmc/lib/cell.py

@@ -268,7 +268,7 @@ class Cell(_FortranObjectWithID):
             return rotation_data[9:]
         else:
             raise ValueError(
-                'Invalid size of rotation matrix: {}'.format(rot_size))
+                f'Invalid size of rotation matrix: {rot_size}')
 
     @rotation.setter
     def rotation(self, rotation_data):

openmc/lib/core.py

@@ -95,6 +95,11 @@ _dll.openmc_simulation_finalize.errcheck = _error_handler
 _dll.openmc_statepoint_write.argtypes = [c_char_p, POINTER(c_bool)]
 _dll.openmc_statepoint_write.restype = c_int
 _dll.openmc_statepoint_write.errcheck = _error_handler
+_dll.openmc_statepoint_load.argtypes = [c_char_p]
+_dll.openmc_statepoint_load.restype = c_int
+_dll.openmc_statepoint_load.errcheck = _error_handler
+_dll.openmc_statepoint_write.restype = c_int
+_dll.openmc_statepoint_write.errcheck = _error_handler
 _dll.openmc_global_bounding_box.argtypes = [POINTER(c_double),
                                             POINTER(c_double)]
 _dll.openmc_global_bounding_box.restype = c_int
@@ -568,6 +573,19 @@ def statepoint_write(filename=None, write_source=True):
     _dll.openmc_statepoint_write(filename, c_bool(write_source))
 
 
+def statepoint_load(filename: PathLike):
+    """Load a statepoint file.
+
+    Parameters
+    ----------
+    filename : path-like
+        Path to the statepoint to load.
+
+    """
+    filename = c_char_p(str(filename).encode())
+    _dll.openmc_statepoint_load(filename)
+
+
 @contextmanager
 def run_in_memory(**kwargs):
     """Provides context manager for calling OpenMC shared library functions.
@@ -611,7 +629,7 @@ class _DLLGlobal:
 
 class _FortranObject:
     def __repr__(self):
-        return "<{}(index={})>".format(type(self).__name__, self._index)
+        return f"<{type(self).__name__}(index={self._index})>"
 
 
 class _FortranObjectWithID(_FortranObject):
@@ -623,7 +641,7 @@ class _FortranObjectWithID(_FortranObject):
         self.id
 
     def __repr__(self):
-        return "<{}(id={})>".format(type(self).__name__, self.id)
+        return f"<{type(self).__name__}(id={self.id})>"
 
 
 @contextmanager

openmc/lib/error.py

@@ -37,5 +37,5 @@ def _error_handler(err, func, args):
         warn(msg)
     elif err < 0:
         if not msg:
-            msg = "Unknown error encountered (code {}).".format(err)
+            msg = f"Unknown error encountered (code {err})."
         raise exc.OpenMCError(msg)

openmc/lib/filter.py

@@ -20,7 +20,8 @@ __all__ = [
     'Filter', 'AzimuthalFilter', 'CellFilter', 'CellbornFilter', 'CellfromFilter',
     'CellInstanceFilter', 'CollisionFilter', 'DistribcellFilter', 'DelayedGroupFilter',
     'EnergyFilter', 'EnergyoutFilter', 'EnergyFunctionFilter', 'LegendreFilter',
-    'MaterialFilter', 'MaterialFromFilter', 'MeshFilter', 'MeshSurfaceFilter', 'MuFilter', 'ParticleFilter',
+    'MaterialFilter', 'MaterialFromFilter', 'MeshFilter', 'MeshBornFilter',
+    'MeshSurfaceFilter', 'MuFilter', 'ParticleFilter',
     'PolarFilter', 'SphericalHarmonicsFilter', 'SpatialLegendreFilter', 'SurfaceFilter',
     'UniverseFilter', 'ZernikeFilter', 'ZernikeRadialFilter', 'filters'
 ]
@@ -89,18 +90,36 @@ _dll.openmc_mesh_filter_get_mesh.errcheck = _error_handler
 _dll.openmc_mesh_filter_set_mesh.argtypes = [c_int32, c_int32]
 _dll.openmc_mesh_filter_set_mesh.restype = c_int
 _dll.openmc_mesh_filter_set_mesh.errcheck = _error_handler
-_dll.openmc_meshsurface_filter_get_mesh.argtypes = [c_int32, POINTER(c_int32)]
-_dll.openmc_meshsurface_filter_get_mesh.restype = c_int
-_dll.openmc_meshsurface_filter_get_mesh.errcheck = _error_handler
-_dll.openmc_meshsurface_filter_set_mesh.argtypes = [c_int32, c_int32]
-_dll.openmc_meshsurface_filter_set_mesh.restype = c_int
-_dll.openmc_meshsurface_filter_set_mesh.errcheck = _error_handler
 _dll.openmc_mesh_filter_get_translation.argtypes = [c_int32, POINTER(c_double*3)]
 _dll.openmc_mesh_filter_get_translation.restype = c_int
 _dll.openmc_mesh_filter_get_translation.errcheck = _error_handler
 _dll.openmc_mesh_filter_set_translation.argtypes = [c_int32, POINTER(c_double*3)]
 _dll.openmc_mesh_filter_set_translation.restype = c_int
 _dll.openmc_mesh_filter_set_translation.errcheck = _error_handler
+_dll.openmc_meshborn_filter_get_mesh.argtypes = [c_int32, POINTER(c_int32)]
+_dll.openmc_meshborn_filter_get_mesh.restype = c_int
+_dll.openmc_meshborn_filter_get_mesh.errcheck = _error_handler
+_dll.openmc_meshborn_filter_set_mesh.argtypes = [c_int32, c_int32]
+_dll.openmc_meshborn_filter_set_mesh.restype = c_int
+_dll.openmc_meshborn_filter_set_mesh.errcheck = _error_handler
+_dll.openmc_meshborn_filter_get_translation.argtypes = [c_int32, POINTER(c_double*3)]
+_dll.openmc_meshborn_filter_get_translation.restype = c_int
+_dll.openmc_meshborn_filter_get_translation.errcheck = _error_handler
+_dll.openmc_meshborn_filter_set_translation.argtypes = [c_int32, POINTER(c_double*3)]
+_dll.openmc_meshborn_filter_set_translation.restype = c_int
+_dll.openmc_meshborn_filter_set_translation.errcheck = _error_handler
+_dll.openmc_meshsurface_filter_get_mesh.argtypes = [c_int32, POINTER(c_int32)]
+_dll.openmc_meshsurface_filter_get_mesh.restype = c_int
+_dll.openmc_meshsurface_filter_get_mesh.errcheck = _error_handler
+_dll.openmc_meshsurface_filter_set_mesh.argtypes = [c_int32, c_int32]
+_dll.openmc_meshsurface_filter_set_mesh.restype = c_int
+_dll.openmc_meshsurface_filter_set_mesh.errcheck = _error_handler
+_dll.openmc_meshsurface_filter_get_translation.argtypes = [c_int32, POINTER(c_double*3)]
+_dll.openmc_meshsurface_filter_get_translation.restype = c_int
+_dll.openmc_meshsurface_filter_get_translation.errcheck = _error_handler
+_dll.openmc_meshsurface_filter_set_translation.argtypes = [c_int32, POINTER(c_double*3)]
+_dll.openmc_meshsurface_filter_set_translation.restype = c_int
+_dll.openmc_meshsurface_filter_set_translation.errcheck = _error_handler
 _dll.openmc_new_filter.argtypes = [c_char_p, POINTER(c_int32)]
 _dll.openmc_new_filter.restype = c_int
 _dll.openmc_new_filter.errcheck = _error_handler
@@ -347,6 +366,34 @@ class MaterialFromFilter(Filter):
 
 
 class MeshFilter(Filter):
+    """Mesh filter stored internally.
+
+    This class exposes a Mesh filter that is stored internally in the OpenMC
+    library. To obtain a view of a Mesh filter with a given ID, use the
+    :data:`openmc.lib.filters` mapping.
+
+    Parameters
+    ----------
+    mesh : openmc.lib.Mesh
+        Mesh to use for the filter
+    uid : int or None
+        Unique ID of the Mesh filter
+    new : bool
+        When `index` is None, this argument controls whether a new object is
+        created or a view of an existing object is returned.
+    index : int
+        Index in the `filters` array.
+
+    Attributes
+    ----------
+    filter_type : str
+        Type of filter
+    mesh : openmc.lib.Mesh
+        Mesh used for the filter
+    translation : Iterable of float
+        3-D coordinates of the translation vector
+
+    """
     filter_type = 'mesh'
 
     def __init__(self, mesh=None, uid=None, new=True, index=None):
@@ -375,7 +422,92 @@ class MeshFilter(Filter):
         _dll.openmc_mesh_filter_set_translation(self._index, (c_double*3)(*translation))
 
 
+class MeshBornFilter(Filter):
+    """MeshBorn filter stored internally.
+
+    This class exposes a MeshBorn filter that is stored internally in the OpenMC
+    library. To obtain a view of a MeshBorn filter with a given ID, use the
+    :data:`openmc.lib.filters` mapping.
+
+    Parameters
+    ----------
+    mesh : openmc.lib.Mesh
+        Mesh to use for the filter
+    uid : int or None
+        Unique ID of the MeshBorn filter
+    new : bool
+        When `index` is None, this argument controls whether a new object is
+        created or a view of an existing object is returned.
+    index : int
+        Index in the `filters` array.
+
+    Attributes
+    ----------
+    filter_type : str
+        Type of filter
+    mesh : openmc.lib.Mesh
+        Mesh used for the filter
+    translation : Iterable of float
+        3-D coordinates of the translation vector
+
+    """
+    filter_type = 'meshborn'
+
+    def __init__(self, mesh=None, uid=None, new=True, index=None):
+        super().__init__(uid, new, index)
+        if mesh is not None:
+            self.mesh = mesh
+
+    @property
+    def mesh(self):
+        index_mesh = c_int32()
+        _dll.openmc_meshborn_filter_get_mesh(self._index, index_mesh)
+        return _get_mesh(index_mesh.value)
+
+    @mesh.setter
+    def mesh(self, mesh):
+        _dll.openmc_meshborn_filter_set_mesh(self._index, mesh._index)
+
+    @property
+    def translation(self):
+        translation = (c_double*3)()
+        _dll.openmc_meshborn_filter_get_translation(self._index, translation)
+        return tuple(translation)
+
+    @translation.setter
+    def translation(self, translation):
+        _dll.openmc_meshborn_filter_set_translation(self._index, (c_double*3)(*translation))
+
+
 class MeshSurfaceFilter(Filter):
+    """MeshSurface filter stored internally.
+
+    This class exposes a MeshSurface filter that is stored internally in the
+    OpenMC library. To obtain a view of a MeshSurface filter with a given ID,
+    use the :data:`openmc.lib.filters` mapping.
+
+    Parameters
+    ----------
+    mesh : openmc.lib.Mesh
+        Mesh to use for the filter
+    uid : int or None
+        Unique ID of the MeshSurface filter
+    new : bool
+        When `index` is None, this argument controls whether a new object is
+        created or a view of an existing object is returned.
+    index : int
+        Index in the `filters` array.
+
+    Attributes
+    ----------
+    filter_type : str
+        Type of filter
+    mesh : openmc.lib.Mesh
+        Mesh used for the filter
+    translation : Iterable of float
+        3-D coordinates of the translation vector
+
+    """
     filter_type = 'meshsurface'
 
     def __init__(self, mesh=None, uid=None, new=True, index=None):
@@ -396,12 +528,12 @@ class MeshSurfaceFilter(Filter):
     @property
     def translation(self):
         translation = (c_double*3)()
-        _dll.openmc_mesh_filter_get_translation(self._index, translation)
+        _dll.openmc_meshsurface_filter_get_translation(self._index, translation)
         return tuple(translation)
 
     @translation.setter
     def translation(self, translation):
-        _dll.openmc_mesh_filter_set_translation(self._index, (c_double*3)(*translation))
+        _dll.openmc_meshsurface_filter_set_translation(self._index, (c_double*3)(*translation))
 
 
 class MuFilter(Filter):
@@ -507,6 +639,7 @@ _FILTER_TYPE_MAP = {
     'material': MaterialFilter,
     'materialfrom': MaterialFromFilter,
     'mesh': MeshFilter,
+    'meshborn': MeshBornFilter,
     'meshsurface': MeshSurfaceFilter,
     'mu': MuFilter,
     'particle': ParticleFilter,

openmc/lib/mesh.py

@@ -15,7 +15,10 @@ from .error import _error_handler
 from .material import Material
 from .plot import _Position
 
-__all__ = ['RegularMesh', 'RectilinearMesh', 'CylindricalMesh', 'SphericalMesh', 'UnstructuredMesh', 'meshes']
+__all__ = [
+    'Mesh', 'RegularMesh', 'RectilinearMesh', 'CylindricalMesh',
+    'SphericalMesh', 'UnstructuredMesh', 'meshes'
+]
 
 
 class _MaterialVolume(Structure):
@@ -553,6 +556,7 @@ class CylindricalMesh(Mesh):
         _dll.openmc_cylindrical_mesh_set_grid(self._index, r_grid, nr, phi_grid,
                                               nphi, z_grid, nz)
 
+
 class SphericalMesh(Mesh):
     """SphericalMesh stored internally.
 

openmc/lib/plot.py

@@ -31,7 +31,7 @@ class _Position(Structure):
         elif idx == 2:
             return self.z
         else:
-            raise IndexError("{} index is invalid for _Position".format(idx))
+            raise IndexError(f"{idx} index is invalid for _Position")
 
     def __setitem__(self, idx, val):
         if idx == 0:
@@ -41,10 +41,10 @@ class _Position(Structure):
         elif idx == 2:
             self.z = val
         else:
-            raise IndexError("{} index is invalid for _Position".format(idx))
+            raise IndexError(f"{idx} index is invalid for _Position")
 
     def __repr__(self):
-        return "({}, {}, {})".format(self.x, self.y, self.z)
+        return f"({self.x}, {self.y}, {self.z})"
 
 
 class _PlotBase(Structure):
@@ -127,7 +127,7 @@ class _PlotBase(Structure):
         elif self.basis_ == 3:
             return 'yz'
 
-        raise ValueError("Plot basis {} is invalid".format(self.basis_))
+        raise ValueError(f"Plot basis {self.basis_} is invalid")
 
     @basis.setter
     def basis(self, basis):
@@ -135,7 +135,7 @@ class _PlotBase(Structure):
             valid_bases = ('xy', 'xz', 'yz')
             basis = basis.lower()
             if basis not in valid_bases:
-                raise ValueError("{} is not a valid plot basis.".format(basis))
+                raise ValueError(f"{basis} is not a valid plot basis.")
 
             if basis == 'xy':
                 self.basis_ = 1
@@ -148,12 +148,11 @@ class _PlotBase(Structure):
         if isinstance(basis, int):
             valid_bases = (1, 2, 3)
             if basis not in valid_bases:
-                raise ValueError("{} is not a valid plot basis.".format(basis))
+                raise ValueError(f"{basis} is not a valid plot basis.")
             self.basis_ = basis
             return
 
-        raise ValueError("{} of type {} is an"
-                         " invalid plot basis".format(basis, type(basis)))
+        raise ValueError(f"{basis} of type {type(basis)} is an invalid plot basis")
 
     @property
     def h_res(self):
@@ -199,14 +198,14 @@ class _PlotBase(Structure):
         out_str = ["-----",
                    "Plot:",
                    "-----",
-                   "Origin: {}".format(self.origin),
-                   "Width: {}".format(self.width),
-                   "Height: {}".format(self.height),
-                   "Basis: {}".format(self.basis),
-                   "HRes: {}".format(self.h_res),
-                   "VRes: {}".format(self.v_res),
-                   "Color Overlaps: {}".format(self.color_overlaps),
-                   "Level: {}".format(self.level)]
+                   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)
 
 

openmc/lib/settings.py

@@ -53,7 +53,7 @@ class _Settings:
                 current_idx.value = idx
                 break
         else:
-            raise ValueError('Invalid run mode: {}'.format(mode))
+            raise ValueError(f'Invalid run mode: {mode}')
 
     @property
     def path_statepoint(self):

openmc/material.py

@@ -152,7 +152,7 @@ class Material(IDManagerMixin):
 
         for nuclide, percent, percent_type in self._nuclides:
             string += '{: <16}'.format('\t{}'.format(nuclide))
-            string += '=\t{: <12} [{}]\n'.format(percent, percent_type)
+            string += f'=\t{percent: <12} [{percent_type}]\n'
 
         if self._macroscopic is not None:
             string += '{: <16}\n'.format('\tMacroscopic Data')
@@ -469,8 +469,7 @@ class Material(IDManagerMixin):
                 raise ValueError('No volume information found for material ID={}.'
                                  .format(self.id))
         else:
-            raise ValueError('No volume information found for material ID={}.'
-                             .format(self.id))
+            raise ValueError(f'No volume information found for material ID={self.id}.')
 
     def set_density(self, units: str, density: Optional[float] = None):
         """Set the density of the material
@@ -500,7 +499,7 @@ class Material(IDManagerMixin):
                       '"sum" unit'.format(self.id)
                 raise ValueError(msg)
 
-            cv.check_type('the density for Material ID="{}"'.format(self.id),
+            cv.check_type(f'the density for Material ID="{self.id}"',
                           density, Real)
             self._density = density
 
@@ -743,20 +742,18 @@ class Material(IDManagerMixin):
             el = element.lower()
             element = openmc.data.ELEMENT_SYMBOL.get(el)
             if element is None:
-                msg = 'Element name "{}" not recognised'.format(el)
+                msg = f'Element name "{el}" not recognised'
                 raise ValueError(msg)
         else:
             if element[0].islower():
-                msg = 'Element name "{}" should start with an uppercase ' \
-                      'letter'.format(element)
+                msg = f'Element name "{element}" should start with an uppercase letter'
                 raise ValueError(msg)
             if len(element) == 2 and element[1].isupper():
-                msg = 'Element name "{}" should end with a lowercase ' \
-                      'letter'.format(element)
+                msg = f'Element name "{element}" should end with a lowercase letter'
                 raise ValueError(msg)
             # skips the first entry of ATOMIC_SYMBOL which is n for neutron
             if element not in list(openmc.data.ATOMIC_SYMBOL.values())[1:]:
-                msg = 'Element name "{}" not recognised'.format(element)
+                msg = f'Element name "{element}" not recognised'
                 raise ValueError(msg)
 
         if self._macroscopic is not None:
@@ -847,8 +844,7 @@ class Material(IDManagerMixin):
             for token in row:
                 if token.isalpha():
                     if token == "n" or token not in openmc.data.ATOMIC_NUMBER:
-                        msg = 'Formula entry {} not an element symbol.' \
-                              .format(token)
+                        msg = f'Formula entry {token} not an element symbol.'
                         raise ValueError(msg)
                 elif token not in ['(', ')', ''] and not token.isdigit():
                         msg = 'Formula must be made from a sequence of ' \
@@ -1373,8 +1369,7 @@ class Material(IDManagerMixin):
                 subelement.set("value", str(self._density))
             subelement.set("units", self._density_units)
         else:
-            raise ValueError('Density has not been set for material {}!'
-                             .format(self.id))
+            raise ValueError(f'Density has not been set for material {self.id}!')
 
         if self._macroscopic is None:
             # Create nuclide XML subelements
@@ -1480,7 +1475,7 @@ class Material(IDManagerMixin):
 
         # Create the new material with the desired name
         if name is None:
-            name = '-'.join(['{}({})'.format(m.name, f) for m, f in
+            name = '-'.join([f'{m.name}({f})' for m, f in
                              zip(materials, fracs)])
         new_mat = openmc.Material(name=name)
 

openmc/mesh.py

@@ -3,10 +3,12 @@ import typing
 import warnings
 from abc import ABC, abstractmethod
 from collections.abc import Iterable
+from functools import wraps
 from math import pi, sqrt, atan2
 from numbers import Integral, Real
 from pathlib import Path
-from typing import Optional, Sequence, Tuple
+import tempfile
+from typing import Optional, Sequence, Tuple, List
 
 import h5py
 import lxml.etree as ET
@@ -15,6 +17,7 @@ import numpy as np
 import openmc
 import openmc.checkvalue as cv
 from openmc.checkvalue import PathLike
+from openmc.utility_funcs import change_directory
 from ._xml import get_text
 from .mixin import IDManagerMixin
 from .surface import _BOUNDARY_TYPES
@@ -39,7 +42,8 @@ class MeshBase(IDManagerMixin, ABC):
     bounding_box : openmc.BoundingBox
         Axis-aligned bounding box of the mesh as defined by the upper-right and
         lower-left coordinates.
-
+    indices : Iterable of tuple
+        An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1), (2, 1, 1), ...]
     """
 
     next_id = 1
@@ -66,6 +70,11 @@ class MeshBase(IDManagerMixin, ABC):
     def bounding_box(self) -> openmc.BoundingBox:
         return openmc.BoundingBox(self.lower_left, self.upper_right)
 
+    @property
+    @abstractmethod
+    def indices(self):
+        pass
+
     def __repr__(self):
         string = type(self).__name__ + '\n'
         string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
@@ -138,6 +147,94 @@ class MeshBase(IDManagerMixin, ABC):
         else:
             raise ValueError(f'Unrecognized mesh type "{mesh_type}" found.')
 
+    def get_homogenized_materials(
+            self,
+            model: openmc.Model,
+            n_samples: int = 10_000,
+            prn_seed: Optional[int] = None,
+            **kwargs
+    ) -> List[openmc.Material]:
+        """Generate homogenized materials over each element in a mesh.
+
+        .. versionadded:: 0.14.1
+
+        Parameters
+        ----------
+        model : openmc.Model
+            Model containing materials to be homogenized and the associated
+            geometry.
+        n_samples : int
+            Number of samples in each mesh element.
+        prn_seed : int, optional
+            Pseudorandom number generator (PRNG) seed; if None, one will be
+            generated randomly.
+        **kwargs
+            Keyword-arguments passed to :func:`openmc.lib.init`.
+
+        Returns
+        -------
+        list of openmc.Material
+            Homogenized material in each mesh element
+
+        """
+        import openmc.lib
+
+        with change_directory(tmpdir=True):
+            # In order to get mesh into model, we temporarily replace the
+            # tallies with a single mesh tally using the current mesh
+            original_tallies = model.tallies
+            new_tally = openmc.Tally()
+            new_tally.filters = [openmc.MeshFilter(self)]
+            new_tally.scores = ['flux']
+            model.tallies = [new_tally]
+
+            # Export model to XML
+            model.export_to_model_xml()
+
+            # Get material volume fractions
+            openmc.lib.init(**kwargs)
+            mesh = openmc.lib.tallies[new_tally.id].filters[0].mesh
+            mat_volume_by_element = [
+                [
+                    (mat.id if mat is not None else None, volume)
+                    for mat, volume in mat_volume_list
+                ]
+                for mat_volume_list in mesh.material_volumes(n_samples, prn_seed)
+            ]
+            openmc.lib.finalize()
+
+            # Restore original tallies
+            model.tallies = original_tallies
+
+        # Create homogenized material for each element
+        materials = model.geometry.get_all_materials()
+        homogenized_materials = []
+        for mat_volume_list in mat_volume_by_element:
+            material_ids, volumes = [list(x) for x in zip(*mat_volume_list)]
+            total_volume = sum(volumes)
+
+            # Check for void material and remove
+            try:
+                index_void = material_ids.index(None)
+            except ValueError:
+                pass
+            else:
+                material_ids.pop(index_void)
+                volumes.pop(index_void)
+
+            # Compute volume fractions
+            volume_fracs = np.array(volumes) / total_volume
+
+            # Get list of materials and mix 'em up!
+            mats = [materials[uid] for uid in material_ids]
+            homogenized_mat = openmc.Material.mix_materials(
+                mats, volume_fracs, 'vo'
+            )
+            homogenized_mat.volume = total_volume
+            homogenized_materials.append(homogenized_mat)
+
+        return homogenized_materials
+
 
 class StructuredMesh(MeshBase):
     """A base class for structured mesh functionality
@@ -1723,6 +1820,8 @@ class SphericalMesh(StructuredMesh):
     @r_grid.setter
     def r_grid(self, grid):
         cv.check_type('mesh r_grid', grid, Iterable, Real)
+        cv.check_length('mesh r_grid', grid, 2)
+        cv.check_increasing('mesh r_grid', grid)
         self._r_grid = np.asarray(grid, dtype=float)
 
     @property
@@ -1732,6 +1831,8 @@ class SphericalMesh(StructuredMesh):
     @theta_grid.setter
     def theta_grid(self, grid):
         cv.check_type('mesh theta_grid', grid, Iterable, Real)
+        cv.check_length('mesh theta_grid', grid, 2)
+        cv.check_increasing('mesh theta_grid', grid)
         self._theta_grid = np.asarray(grid, dtype=float)
 
     @property
@@ -1741,6 +1842,8 @@ class SphericalMesh(StructuredMesh):
     @phi_grid.setter
     def phi_grid(self, grid):
         cv.check_type('mesh phi_grid', grid, Iterable, Real)
+        cv.check_length('mesh phi_grid', grid, 2)
+        cv.check_increasing('mesh phi_grid', grid)
         self._phi_grid = np.asarray(grid, dtype=float)
 
     @property
@@ -1914,6 +2017,17 @@ class SphericalMesh(StructuredMesh):
         return arr
 
 
+def require_statepoint_data(func):
+    @wraps(func)
+    def wrapper(self: UnstructuredMesh, *args, **kwargs):
+        if not self._has_statepoint_data:
+            raise AttributeError(f'The "{func.__name__}" property requires '
+                                 'information about this mesh to be loaded '
+                                 'from a statepoint file.')
+        return func(self, *args, **kwargs)
+    return wrapper
+
+
 class UnstructuredMesh(MeshBase):
     """A 3D unstructured mesh
 
@@ -1934,6 +2048,11 @@ class UnstructuredMesh(MeshBase):
         Name of the mesh
     length_multiplier: float
         Constant multiplier to apply to mesh coordinates
+    options : str, optional
+        Special options that control spatial search data structures used. This
+        is currently only used to set `parameters
+        <https://tinyurl.com/kdtree-params>`_ for MOAB's AdaptiveKDTree. If
+        None, OpenMC internally uses a default of "MAX_DEPTH=20;PLANE_SET=2;".
 
     Attributes
     ----------
@@ -1947,6 +2066,11 @@ class UnstructuredMesh(MeshBase):
         Multiplicative factor to apply to mesh coordinates
     library : {'moab', 'libmesh'}
         Mesh library used for the unstructured mesh tally
+    options : str
+        Special options that control spatial search data structures used. This
+        is currently only used to set `parameters
+        <https://tinyurl.com/kdtree-params>`_ for MOAB's AdaptiveKDTree. If
+        None, OpenMC internally uses a default of "MAX_DEPTH=20;PLANE_SET=2;".
     output : bool
         Indicates whether or not automatic tally output should be generated for
         this mesh
@@ -1980,7 +2104,8 @@ class UnstructuredMesh(MeshBase):
     _LINEAR_HEX = 1
 
     def __init__(self, filename: PathLike, library: str, mesh_id: Optional[int] = None,
-                 name: str = '', length_multiplier: float = 1.0):
+                 name: str = '', length_multiplier: float = 1.0,
+                 options: Optional[str] = None):
         super().__init__(mesh_id, name)
         self.filename = filename
         self._volumes = None
@@ -1990,6 +2115,8 @@ class UnstructuredMesh(MeshBase):
         self.library = library
         self._output = False
         self.length_multiplier = length_multiplier
+        self.options = options
+        self._has_statepoint_data = False
 
     @property
     def filename(self):
@@ -2010,6 +2137,16 @@ class UnstructuredMesh(MeshBase):
         self._library = lib
 
     @property
+    def options(self) -> Optional[str]:
+        return self._options
+
+    @options.setter
+    def options(self, options: Optional[str]):
+        cv.check_type('options', options, (str, type(None)))
+        self._options = options
+
+    @property
+    @require_statepoint_data
     def size(self):
         return self._size
 
@@ -2028,6 +2165,7 @@ class UnstructuredMesh(MeshBase):
         self._output = val
 
     @property
+    @require_statepoint_data
     def volumes(self):
         """Return Volumes for every mesh cell if
         populated by a StatePoint file
@@ -2046,26 +2184,32 @@ class UnstructuredMesh(MeshBase):
         self._volumes = volumes
 
     @property
+    @require_statepoint_data
     def total_volume(self):
         return np.sum(self.volumes)
 
     @property
+    @require_statepoint_data
     def vertices(self):
         return self._vertices
 
     @property
+    @require_statepoint_data
     def connectivity(self):
         return self._connectivity
 
     @property
+    @require_statepoint_data
     def element_types(self):
         return self._element_types
 
     @property
+    @require_statepoint_data
     def centroids(self):
         return np.array([self.centroid(i) for i in range(self.n_elements)])
 
     @property
+    @require_statepoint_data
     def n_elements(self):
         if self._n_elements is None:
             raise RuntimeError("No information about this mesh has "
@@ -2096,6 +2240,15 @@ class UnstructuredMesh(MeshBase):
     def n_dimension(self):
         return 3
 
+    @property
+    @require_statepoint_data
+    def indices(self):
+        return [(i,) for i in range(self.n_elements)]
+
+    @property
+    def has_statepoint_data(self) -> bool:
+        return self._has_statepoint_data
+
     def __repr__(self):
         string = super().__repr__()
         string += '{: <16}=\t{}\n'.format('\tFilename', self.filename)
@@ -2103,16 +2256,21 @@ class UnstructuredMesh(MeshBase):
         if self.length_multiplier != 1.0:
             string += '{: <16}=\t{}\n'.format('\tLength multiplier',
                                               self.length_multiplier)
+        if self.options is not None:
+            string += '{: <16}=\t{}\n'.format('\tOptions', self.options)
         return string
 
     @property
+    @require_statepoint_data
     def lower_left(self):
         return self.vertices.min(axis=0)
 
     @property
+    @require_statepoint_data
     def upper_right(self):
         return self.vertices.max(axis=0)
 
+    @require_statepoint_data
     def centroid(self, bin: int):
         """Return the vertex averaged centroid of an element
 
@@ -2255,8 +2413,13 @@ class UnstructuredMesh(MeshBase):
         mesh_id = int(group.name.split('/')[-1].lstrip('mesh '))
         filename = group['filename'][()].decode()
         library = group['library'][()].decode()
+        if 'options' in group.attrs:
+            options = group.attrs['options'].decode()
+        else:
+            options = None
 
-        mesh = cls(filename=filename, library=library, mesh_id=mesh_id)
+        mesh = cls(filename=filename, library=library, mesh_id=mesh_id, options=options)
+        mesh._has_statepoint_data = True
         vol_data = group['volumes'][()]
         mesh.volumes = np.reshape(vol_data, (vol_data.shape[0],))
         mesh.n_elements = mesh.volumes.size
@@ -2286,6 +2449,8 @@ class UnstructuredMesh(MeshBase):
         element.set("id", str(self._id))
         element.set("type", "unstructured")
         element.set("library", self._library)
+        if self.options is not None:
+            element.set('options', self.options)
         subelement = ET.SubElement(element, "filename")
         subelement.text = str(self.filename)
 
@@ -2312,8 +2477,9 @@ class UnstructuredMesh(MeshBase):
         filename = get_text(elem, 'filename')
         library = get_text(elem, 'library')
         length_multiplier = float(get_text(elem, 'length_multiplier', 1.0))
+        options = elem.get('options')
 
-        return cls(filename, library, mesh_id, '', length_multiplier)
+        return cls(filename, library, mesh_id, '', length_multiplier, options)
 
 
 def _read_meshes(elem):

openmc/mgxs/library.py

@@ -685,7 +685,7 @@ class Library:
 
         # Check that requested domain is included in library
         if mgxs_type not in self.mgxs_types:
-            msg = 'Unable to find MGXS type "{0}"'.format(mgxs_type)
+            msg = f'Unable to find MGXS type "{mgxs_type}"'
             raise ValueError(msg)
 
         return self.all_mgxs[domain_id][mgxs_type]
@@ -901,7 +901,7 @@ class Library:
         if not os.path.exists(directory):
             os.makedirs(directory)
 
-        full_filename = os.path.join(directory, '{}.pkl'.format(filename))
+        full_filename = os.path.join(directory, f'{filename}.pkl')
         full_filename = full_filename.replace(' ', '-')
 
         # Load and return pickled Library object