diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index ab2121b9c9..99d2155fad 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -93,10 +93,10 @@ jobs:
       RDMAV_FORK_SAFE: 1
 
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
 
       - name: Set up Python ${{ matrix.python-version }}
-        uses: actions/setup-python@v4
+        uses: actions/setup-python@v5
         with:
           python-version: ${{ matrix.python-version }}
 
@@ -128,7 +128,7 @@ jobs:
           $GITHUB_WORKSPACE/tools/ci/gha-install.sh
 
       - name: cache-xs
-        uses: actions/cache@v3
+        uses: actions/cache@v4
         with:
           path: |
             ~/nndc_hdf5
@@ -156,7 +156,7 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - name: Coveralls Finished
-        uses: coverallsapp/github-action@master
+        uses: coverallsapp/github-action@v2
         with:
           github-token: ${{ secrets.github_token }}
           parallel-finished: true
diff --git a/.github/workflows/dockerhub-publish-dagmc-libmesh.yml b/.github/workflows/dockerhub-publish-dagmc-libmesh.yml
index 3f0b5ab9be..813596953b 100644
--- a/.github/workflows/dockerhub-publish-dagmc-libmesh.yml
+++ b/.github/workflows/dockerhub-publish-dagmc-libmesh.yml
@@ -10,20 +10,20 @@ jobs:
     steps:
       -
         name: Set up QEMU
-        uses: docker/setup-qemu-action@v2
+        uses: docker/setup-qemu-action@v3
       -
         name: Set up Docker Buildx
-        uses: docker/setup-buildx-action@v2
+        uses: docker/setup-buildx-action@v3
       -
         name: Login to DockerHub
-        uses: docker/login-action@v2
+        uses: docker/login-action@v3
         with:
           username: ${{ secrets.DOCKERHUB_USERNAME }}
           password: ${{ secrets.DOCKERHUB_TOKEN }}
       -
         name: Build and push
         id: docker_build
-        uses: docker/build-push-action@v4
+        uses: docker/build-push-action@v5
         with:
           push: true
           tags: openmc/openmc:latest-dagmc-libmesh
diff --git a/.github/workflows/dockerhub-publish-dagmc.yml b/.github/workflows/dockerhub-publish-dagmc.yml
index b6ba2f075c..6757f77271 100644
--- a/.github/workflows/dockerhub-publish-dagmc.yml
+++ b/.github/workflows/dockerhub-publish-dagmc.yml
@@ -10,20 +10,20 @@ jobs:
     steps:
       -
         name: Set up QEMU
-        uses: docker/setup-qemu-action@v2
+        uses: docker/setup-qemu-action@v3
       -
         name: Set up Docker Buildx
-        uses: docker/setup-buildx-action@v2
+        uses: docker/setup-buildx-action@v3
       -
         name: Login to DockerHub
-        uses: docker/login-action@v2 
+        uses: docker/login-action@v3 
         with:
           username: ${{ secrets.DOCKERHUB_USERNAME }}
           password: ${{ secrets.DOCKERHUB_TOKEN }}
       -
         name: Build and push
         id: docker_build
-        uses: docker/build-push-action@v4
+        uses: docker/build-push-action@v5
         with:
           push: true
           tags: openmc/openmc:latest-dagmc
diff --git a/.github/workflows/dockerhub-publish-dev.yml b/.github/workflows/dockerhub-publish-dev.yml
index d2603ead4a..7a81363a78 100644
--- a/.github/workflows/dockerhub-publish-dev.yml
+++ b/.github/workflows/dockerhub-publish-dev.yml
@@ -10,20 +10,20 @@ jobs:
     steps:
       -
         name: Set up QEMU
-        uses: docker/setup-qemu-action@v2
+        uses: docker/setup-qemu-action@v3
       -
         name: Set up Docker Buildx
-        uses: docker/setup-buildx-action@v2
+        uses: docker/setup-buildx-action@v3
       -
         name: Login to DockerHub
-        uses: docker/login-action@v2 
+        uses: docker/login-action@v3 
         with:
           username: ${{ secrets.DOCKERHUB_USERNAME }}
           password: ${{ secrets.DOCKERHUB_TOKEN }}
       -
         name: Build and push
         id: docker_build
-        uses: docker/build-push-action@v4
+        uses: docker/build-push-action@v5
         with:
           push: true
           tags: openmc/openmc:develop
diff --git a/.github/workflows/dockerhub-publish-develop-dagmc-libmesh.yml b/.github/workflows/dockerhub-publish-develop-dagmc-libmesh.yml
index 354f0a020e..a219f2a91d 100644
--- a/.github/workflows/dockerhub-publish-develop-dagmc-libmesh.yml
+++ b/.github/workflows/dockerhub-publish-develop-dagmc-libmesh.yml
@@ -10,20 +10,20 @@ jobs:
     steps:
       -
         name: Set up QEMU
-        uses: docker/setup-qemu-action@v2
+        uses: docker/setup-qemu-action@v3
       -
         name: Set up Docker Buildx
-        uses: docker/setup-buildx-action@v2
+        uses: docker/setup-buildx-action@v3
       -
         name: Login to DockerHub
-        uses: docker/login-action@v2 
+        uses: docker/login-action@v3 
         with:
           username: ${{ secrets.DOCKERHUB_USERNAME }}
           password: ${{ secrets.DOCKERHUB_TOKEN }}
       -
         name: Build and push
         id: docker_build
-        uses: docker/build-push-action@v4
+        uses: docker/build-push-action@v5
         with:
           push: true
           tags: openmc/openmc:develop-dagmc-libmesh
diff --git a/.github/workflows/dockerhub-publish-develop-dagmc.yml b/.github/workflows/dockerhub-publish-develop-dagmc.yml
index 36ec7a337b..a901b8d3f0 100644
--- a/.github/workflows/dockerhub-publish-develop-dagmc.yml
+++ b/.github/workflows/dockerhub-publish-develop-dagmc.yml
@@ -10,20 +10,20 @@ jobs:
     steps:
       -
         name: Set up QEMU
-        uses: docker/setup-qemu-action@v2
+        uses: docker/setup-qemu-action@v3
       -
         name: Set up Docker Buildx
-        uses: docker/setup-buildx-action@v2
+        uses: docker/setup-buildx-action@v3
       -
         name: Login to DockerHub
-        uses: docker/login-action@v2 
+        uses: docker/login-action@v3 
         with:
           username: ${{ secrets.DOCKERHUB_USERNAME }}
           password: ${{ secrets.DOCKERHUB_TOKEN }}
       -
         name: Build and push
         id: docker_build
-        uses: docker/build-push-action@v4
+        uses: docker/build-push-action@v5
         with:
           push: true
           tags: openmc/openmc:develop-dagmc
diff --git a/.github/workflows/dockerhub-publish-develop-libmesh.yml b/.github/workflows/dockerhub-publish-develop-libmesh.yml
index a89417316b..22e9aa68fb 100644
--- a/.github/workflows/dockerhub-publish-develop-libmesh.yml
+++ b/.github/workflows/dockerhub-publish-develop-libmesh.yml
@@ -10,20 +10,20 @@ jobs:
     steps:
       -
         name: Set up QEMU
-        uses: docker/setup-qemu-action@v2
+        uses: docker/setup-qemu-action@v3
       -
         name: Set up Docker Buildx
-        uses: docker/setup-buildx-action@v2
+        uses: docker/setup-buildx-action@v3
       -
         name: Login to DockerHub
-        uses: docker/login-action@v2 
+        uses: docker/login-action@v3 
         with:
           username: ${{ secrets.DOCKERHUB_USERNAME }}
           password: ${{ secrets.DOCKERHUB_TOKEN }}
       -
         name: Build and push
         id: docker_build
-        uses: docker/build-push-action@v4
+        uses: docker/build-push-action@v5
         with:
           push: true
           tags: openmc/openmc:develop-libmesh
diff --git a/.github/workflows/dockerhub-publish-libmesh.yml b/.github/workflows/dockerhub-publish-libmesh.yml
index e592ccb8ef..843ce0f6fd 100644
--- a/.github/workflows/dockerhub-publish-libmesh.yml
+++ b/.github/workflows/dockerhub-publish-libmesh.yml
@@ -10,20 +10,20 @@ jobs:
     steps:
       -
         name: Set up QEMU
-        uses: docker/setup-qemu-action@v2
+        uses: docker/setup-qemu-action@v3
       -
         name: Set up Docker Buildx
-        uses: docker/setup-buildx-action@v2
+        uses: docker/setup-buildx-action@v3
       -
         name: Login to DockerHub
-        uses: docker/login-action@v2 
+        uses: docker/login-action@v3 
         with:
           username: ${{ secrets.DOCKERHUB_USERNAME }}
           password: ${{ secrets.DOCKERHUB_TOKEN }}
       -
         name: Build and push
         id: docker_build
-        uses: docker/build-push-action@v4
+        uses: docker/build-push-action@v5
         with:
           push: true
           tags: openmc/openmc:latest-libmesh
diff --git a/.github/workflows/dockerhub-publish-release-dagmc-libmesh.yml b/.github/workflows/dockerhub-publish-release-dagmc-libmesh.yml
index c90302af44..db62bb53e6 100644
--- a/.github/workflows/dockerhub-publish-release-dagmc-libmesh.yml
+++ b/.github/workflows/dockerhub-publish-release-dagmc-libmesh.yml
@@ -8,25 +8,25 @@ jobs:
   main:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - name: Set env
         run: echo "RELEASE_VERSION=${GITHUB_REF#refs/*/}" >> $GITHUB_ENV
       -
         name: Set up QEMU
-        uses: docker/setup-qemu-action@v2
+        uses: docker/setup-qemu-action@v3
       -
         name: Set up Docker Buildx
-        uses: docker/setup-buildx-action@v2
+        uses: docker/setup-buildx-action@v3
       -
         name: Login to DockerHub
-        uses: docker/login-action@v2
+        uses: docker/login-action@v3
         with:
           username: ${{ secrets.DOCKERHUB_USERNAME }}
           password: ${{ secrets.DOCKERHUB_TOKEN }}
       -
         name: Build and push
         id: docker_build
-        uses: docker/build-push-action@v4
+        uses: docker/build-push-action@v5
         with:
           push: true
           tags: openmc/openmc:${{ env.RELEASE_VERSION }}-dagmc-libmesh
diff --git a/.github/workflows/dockerhub-publish-release-dagmc.yml b/.github/workflows/dockerhub-publish-release-dagmc.yml
index ef66f6eadd..de95937828 100644
--- a/.github/workflows/dockerhub-publish-release-dagmc.yml
+++ b/.github/workflows/dockerhub-publish-release-dagmc.yml
@@ -8,25 +8,25 @@ jobs:
   main:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - name: Set env
         run: echo "RELEASE_VERSION=${GITHUB_REF#refs/*/}" >> $GITHUB_ENV
       -
         name: Set up QEMU
-        uses: docker/setup-qemu-action@v2
+        uses: docker/setup-qemu-action@v3
       -
         name: Set up Docker Buildx
-        uses: docker/setup-buildx-action@v2
+        uses: docker/setup-buildx-action@v3
       -
         name: Login to DockerHub
-        uses: docker/login-action@v2 
+        uses: docker/login-action@v3 
         with:
           username: ${{ secrets.DOCKERHUB_USERNAME }}
           password: ${{ secrets.DOCKERHUB_TOKEN }}
       -
         name: Build and push
         id: docker_build
-        uses: docker/build-push-action@v4
+        uses: docker/build-push-action@v5
         with:
           push: true
           tags: openmc/openmc:${{ env.RELEASE_VERSION }}-dagmc
diff --git a/.github/workflows/dockerhub-publish-release-libmesh.yml b/.github/workflows/dockerhub-publish-release-libmesh.yml
index 72edfbc68e..e8ea98aebd 100644
--- a/.github/workflows/dockerhub-publish-release-libmesh.yml
+++ b/.github/workflows/dockerhub-publish-release-libmesh.yml
@@ -8,25 +8,25 @@ jobs:
   main:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - name: Set env
         run: echo "RELEASE_VERSION=${GITHUB_REF#refs/*/}" >> $GITHUB_ENV
       -
         name: Set up QEMU
-        uses: docker/setup-qemu-action@v2
+        uses: docker/setup-qemu-action@v3
       -
         name: Set up Docker Buildx
-        uses: docker/setup-buildx-action@v2
+        uses: docker/setup-buildx-action@v3
       -
         name: Login to DockerHub
-        uses: docker/login-action@v2
+        uses: docker/login-action@v3
         with:
           username: ${{ secrets.DOCKERHUB_USERNAME }}
           password: ${{ secrets.DOCKERHUB_TOKEN }}
       -
         name: Build and push
         id: docker_build
-        uses: docker/build-push-action@v4
+        uses: docker/build-push-action@v5
         with:
           push: true
           tags: openmc/openmc:${{ env.RELEASE_VERSION }}-libmesh
diff --git a/.github/workflows/dockerhub-publish-release.yml b/.github/workflows/dockerhub-publish-release.yml
index 22d21b4c96..fab030192b 100644
--- a/.github/workflows/dockerhub-publish-release.yml
+++ b/.github/workflows/dockerhub-publish-release.yml
@@ -8,25 +8,25 @@ jobs:
   main:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - name: Set env
         run: echo "RELEASE_VERSION=${GITHUB_REF#refs/*/}" >> $GITHUB_ENV
       -
         name: Set up QEMU
-        uses: docker/setup-qemu-action@v2
+        uses: docker/setup-qemu-action@v3
       -
         name: Set up Docker Buildx
-        uses: docker/setup-buildx-action@v2
+        uses: docker/setup-buildx-action@v3
       -
         name: Login to DockerHub
-        uses: docker/login-action@v2 
+        uses: docker/login-action@v3 
         with:
           username: ${{ secrets.DOCKERHUB_USERNAME }}
           password: ${{ secrets.DOCKERHUB_TOKEN }}
       -
         name: Build and push
         id: docker_build
-        uses: docker/build-push-action@v4
+        uses: docker/build-push-action@v5
         with:
           push: true
           tags: openmc/openmc:${{ env.RELEASE_VERSION }}
diff --git a/.github/workflows/dockerhub-publish.yml b/.github/workflows/dockerhub-publish.yml
index 00a3ba3164..fd51a9fa73 100644
--- a/.github/workflows/dockerhub-publish.yml
+++ b/.github/workflows/dockerhub-publish.yml
@@ -10,20 +10,20 @@ jobs:
     steps:
       -
         name: Set up QEMU
-        uses: docker/setup-qemu-action@v2
+        uses: docker/setup-qemu-action@v3
       -
         name: Set up Docker Buildx
-        uses: docker/setup-buildx-action@v2
+        uses: docker/setup-buildx-action@v3
       -
         name: Login to DockerHub
-        uses: docker/login-action@v2 
+        uses: docker/login-action@v3 
         with:
           username: ${{ secrets.DOCKERHUB_USERNAME }}
           password: ${{ secrets.DOCKERHUB_TOKEN }}
       -
         name: Build and push
         id: docker_build
-        uses: docker/build-push-action@v4
+        uses: docker/build-push-action@v5
         with:
           push: true
           tags: openmc/openmc:latest
diff --git a/.github/workflows/format-check.yml b/.github/workflows/format-check.yml
index 36ccece3cc..cef14ca2c8 100644
--- a/.github/workflows/format-check.yml
+++ b/.github/workflows/format-check.yml
@@ -13,7 +13,7 @@ jobs:
   cpp-linter:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - uses: cpp-linter/cpp-linter-action@v2
         id: linter
         env:
diff --git a/CMakeLists.txt b/CMakeLists.txt
index ec6ea61792..1280daf2b9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -193,12 +193,10 @@ endif()
 if(NOT MSVC)
 
 if(OPENMC_USE_OPENMP)
-  find_package(OpenMP)
-  if(OPENMP_FOUND)
-    # 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()
+  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)
diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst
index 9b855b2502..da71e4e600 100644
--- a/docs/source/io_formats/settings.rst
+++ b/docs/source/io_formats/settings.rst
@@ -276,6 +276,13 @@ 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. 
+
+  *Default*: 1000
 
 --------------------------------------
 ``<max_write_lost_particles>`` Element
diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst
index 705aecee9f..7b924e3749 100644
--- a/docs/source/pythonapi/base.rst
+++ b/docs/source/pythonapi/base.rst
@@ -35,6 +35,7 @@ Simulation Settings
    :nosignatures:
    :template: myfunction.rst
 
+   openmc.read_source_file
    openmc.write_source_file
    openmc.wwinp_to_wws
 
diff --git a/include/openmc/capi.h b/include/openmc/capi.h
index 063ed8d9fa..5f98152cd6 100644
--- a/include/openmc/capi.h
+++ b/include/openmc/capi.h
@@ -106,6 +106,7 @@ int openmc_mesh_filter_set_translation(int32_t index, double translation[3]);
 int openmc_mesh_get_id(int32_t index, int32_t* id);
 int openmc_mesh_set_id(int32_t index, int32_t id);
 int openmc_mesh_get_n_elements(int32_t index, size_t* n);
+int openmc_mesh_get_volumes(int32_t index, double* volumes);
 int openmc_mesh_material_volumes(int32_t index, int n_sample, int bin,
   int result_size, void* result, int* hits, uint64_t* seed);
 int openmc_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh);
diff --git a/openmc/cell.py b/openmc/cell.py
index 0af6ea2b45..998640930a 100644
--- a/openmc/cell.py
+++ b/openmc/cell.py
@@ -722,7 +722,7 @@ class Cell(IDManagerMixin):
             Dictionary mapping surface IDs to :class:`openmc.Surface` instances
         materials : dict
             Dictionary mapping material ID strings to :class:`openmc.Material`
-            instances (defined in :math:`openmc.Geometry.from_xml`)
+            instances (defined in :meth:`openmc.Geometry.from_xml`)
         get_universe : function
             Function returning universe (defined in
             :meth:`openmc.Geometry.from_xml`)
diff --git a/openmc/data/data.py b/openmc/data/data.py
index 609a901d9b..7caf31e26e 100644
--- a/openmc/data/data.py
+++ b/openmc/data/data.py
@@ -179,6 +179,93 @@ ATOMIC_SYMBOL = {0: 'n', 1: 'H', 2: 'He', 3: 'Li', 4: 'Be', 5: 'B', 6: 'C',
                  118: 'Og'}
 ATOMIC_NUMBER = {value: key for key, value in ATOMIC_SYMBOL.items()}
 
+DADZ = {
+    '(n,2nd)': (-3, -1),
+    '(n,2n)': (-1, 0),
+    '(n,3n)': (-2, 0),
+    '(n,na)': (-4, -2),
+    '(n,n3a)': (-12, -6),
+    '(n,2na)': (-5, -2),
+    '(n,3na)': (-6, -2),
+    '(n,np)': (-1, -1),
+    '(n,n2a)': (-8, -4),
+    '(n,2n2a)': (-9, -4),
+    '(n,nd)': (-2, -1),
+    '(n,nt)': (-3, -1),
+    '(n,n3He)': (-3, -2),
+    '(n,nd2a)': (-10, -5),
+    '(n,nt2a)': (-11, -5),
+    '(n,4n)': (-3, 0),
+    '(n,2np)': (-2, -1),
+    '(n,3np)': (-3, -1),
+    '(n,n2p)': (-2, -2),
+    '(n,npa)': (-5, -3),
+    '(n,gamma)': (1, 0),
+    '(n,p)': (0, -1),
+    '(n,d)': (-1, -1),
+    '(n,t)': (-2, -1),
+    '(n,3He)': (-2, -2),
+    '(n,a)': (-3, -2),
+    '(n,2a)': (-7, -4),
+    '(n,3a)': (-11, -6),
+    '(n,2p)': (-1, -2),
+    '(n,pa)': (-4, -3),
+    '(n,t2a)': (-10, -5),
+    '(n,d2a)': (-9, -5),
+    '(n,pd)': (-2, -2),
+    '(n,pt)': (-3, -2),
+    '(n,da)': (-5, -3),
+    '(n,5n)': (-4, 0),
+    '(n,6n)': (-5, 0),
+    '(n,2nt)': (-4, -1),
+    '(n,ta)': (-6, -3),
+    '(n,4np)': (-4, -1),
+    '(n,3nd)': (-4, -1),
+    '(n,nda)': (-6, -3),
+    '(n,2npa)': (-6, -3),
+    '(n,7n)': (-6, 0),
+    '(n,8n)': (-7, 0),
+    '(n,5np)': (-5, -1),
+    '(n,6np)': (-6, -1),
+    '(n,7np)': (-7, -1),
+    '(n,4na)': (-7, -2),
+    '(n,5na)': (-8, -2),
+    '(n,6na)': (-9, -2),
+    '(n,7na)': (-10, -2),
+    '(n,4nd)': (-5, -1),
+    '(n,5nd)': (-6, -1),
+    '(n,6nd)': (-7, -1),
+    '(n,3nt)': (-5, -1),
+    '(n,4nt)': (-6, -1),
+    '(n,5nt)': (-7, -1),
+    '(n,6nt)': (-8, -1),
+    '(n,2n3He)': (-4, -2),
+    '(n,3n3He)': (-5, -2),
+    '(n,4n3He)': (-6, -2),
+    '(n,3n2p)': (-4, -2),
+    '(n,3n2a)': (-10, -4),
+    '(n,3npa)': (-7, -3),
+    '(n,dt)': (-4, -2),
+    '(n,npd)': (-3, -2),
+    '(n,npt)': (-4, -2),
+    '(n,ndt)': (-5, -2),
+    '(n,np3He)': (-4, -3),
+    '(n,nd3He)': (-5, -3),
+    '(n,nt3He)': (-6, -3),
+    '(n,nta)': (-7, -3),
+    '(n,2n2p)': (-3, -2),
+    '(n,p3He)': (-4, -3),
+    '(n,d3He)': (-5, -3),
+    '(n,3Hea)': (-6, -4),
+    '(n,4n2p)': (-5, -2),
+    '(n,4n2a)': (-11, -4),
+    '(n,4npa)': (-8, -3),
+    '(n,3p)': (-2, -3),
+    '(n,n3p)': (-3, -3),
+    '(n,3n2pa)': (-8, -4),
+    '(n,5n2p)': (-6, -2),
+}
+
 # Values here are from the Committee on Data for Science and Technology
 # (CODATA) 2018 recommendation (https://physics.nist.gov/cuu/Constants/).
 
diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py
index b782ef830d..86275a54ab 100644
--- a/openmc/deplete/abc.py
+++ b/openmc/deplete/abc.py
@@ -3,6 +3,7 @@
 This module contains Abstract Base Classes for implementing operator, integrator, depletion system solver, and operator helper classes
 """
 
+from __future__ import annotations
 from abc import ABC, abstractmethod
 from collections import namedtuple, defaultdict
 from collections.abc import Iterable, Callable
@@ -13,17 +14,20 @@ from contextlib import contextmanager
 import os
 from pathlib import Path
 import time
+from typing import Optional, Union, Sequence
 from warnings import warn
 
-from numpy import nonzero, empty, asarray
+import numpy as np
 from uncertainties import ufloat
 
 from openmc.checkvalue import check_value, check_type, check_greater_than, PathLike
 from openmc.mpi import comm
+from openmc import Material
 from .stepresult import StepResult
 from .chain import Chain
 from .results import Results
 from .pool import deplete
+from .reaction_rates import ReactionRates
 from .transfer_rates import TransferRates
 from openmc import Material, Cell
 from .reactivity_control import (
@@ -181,7 +185,7 @@ class TransportOperator(ABC):
         pass
 
     @abstractmethod
-    def write_bos_data(self, step):
+    def write_bos_data(self, step: int):
         """Document beginning of step data for a given step
 
         Called at the beginning of a depletion step and at
@@ -220,7 +224,7 @@ class ReactionRateHelper(ABC):
 
     def __init__(self, n_nucs, n_react):
         self._nuclides = None
-        self._results_cache = empty((n_nucs, n_react))
+        self._results_cache = np.empty((n_nucs, n_react))
 
     @abstractmethod
     def generate_tallies(self, materials, scores):
@@ -237,7 +241,12 @@ class ReactionRateHelper(ABC):
         self._nuclides = nuclides
 
     @abstractmethod
-    def get_material_rates(self, mat_id, nuc_index, react_index):
+    def get_material_rates(
+        self,
+        mat_id: int,
+        nuc_index: Sequence[str],
+        react_index: Sequence[str]
+    ):
         """Return 2D array of [nuclide, reaction] reaction rates
 
         Parameters
@@ -250,7 +259,7 @@ class ReactionRateHelper(ABC):
             Ordering of reactions
         """
 
-    def divide_by_atoms(self, number):
+    def divide_by_atoms(self, number: Sequence[float]):
         """Normalize reaction rates by number of atoms
 
         Acts on the current material examined by :meth:`get_material_rates`
@@ -267,7 +276,7 @@ class ReactionRateHelper(ABC):
             normalized by the number of nuclides
         """
 
-        mask = nonzero(number)
+        mask = np.nonzero(number)
         results = self._results_cache
         for col in range(results.shape[1]):
             results[mask, col] /= number[mask]
@@ -299,7 +308,7 @@ class NormalizationHelper(ABC):
         """Reset state for normalization"""
 
     @abstractmethod
-    def prepare(self, chain_nucs, rate_index):
+    def prepare(self, chain_nucs: Sequence[str], rate_index: dict):
         """Perform work needed to obtain energy produced
 
         This method is called prior to calculating the reaction rates
@@ -338,7 +347,7 @@ class NormalizationHelper(ABC):
         self._nuclides = nuclides
 
     @abstractmethod
-    def factor(self, source_rate):
+    def factor(self, source_rate: float):
         """Return normalization factor
 
         Parameters
@@ -441,7 +450,7 @@ class FissionYieldHelper(ABC):
             in parallel mode.
         """
 
-    def update_tally_nuclides(self, nuclides):
+    def update_tally_nuclides(self, nuclides: Sequence[str]) -> list:
         """Return nuclides with non-zero densities and yield data
 
         Parameters
@@ -567,8 +576,16 @@ class Integrator(ABC):
 
     """
 
-    def __init__(self, operator, timesteps, power=None, power_density=None,
-                 source_rates=None, timestep_units='s', solver="cram48"):
+    def __init__(
+            self,
+            operator: TransportOperator,
+            timesteps: Sequence[float],
+            power: Optional[Union[float, Sequence[float]]] = None,
+            power_density: Optional[Union[float, Sequence[float]]] = None,
+            source_rates: Optional[Sequence[float]] = None,
+            timestep_units: str = 's',
+            solver: str = "cram48"
+        ):
         # Check number of stages previously used
         if operator.prev_res is not None:
             res = operator.prev_res[-1]
@@ -640,8 +657,8 @@ class Integrator(ABC):
             else:
                 raise ValueError("Invalid timestep unit '{}'".format(unit))
 
-        self.timesteps = asarray(seconds)
-        self.source_rates = asarray(source_rates)
+        self.timesteps = np.asarray(seconds)
+        self.source_rates = np.asarray(source_rates)
 
         self.transfer_rates = None
         self._reactivity_control = None
@@ -705,7 +722,14 @@ class Integrator(ABC):
         return time.time() - start, results
 
     @abstractmethod
-    def __call__(self, n, rates, dt, source_rate, i):
+    def __call__(
+        self,
+        n: Sequence[np.ndarray],
+        rates: ReactionRates,
+        dt: float,
+        source_rate: float,
+        i: int
+    ):
         """Perform the integration across one time step
 
         Parameters
@@ -758,7 +782,7 @@ class Integrator(ABC):
         self.operator.write_bos_data(step_index + self._i_res)
         return x, res
 
-    def _get_bos_data_from_restart(self, step_index, source_rate, bos_conc):
+    def _get_bos_data_from_restart(self, source_rate, bos_conc):
         """Get beginning of step concentrations, reaction rates from restart"""
         res = self.operator.prev_res[-1]
         # Depletion methods expect list of arrays
@@ -860,8 +884,14 @@ class Integrator(ABC):
 
         self.operator.finalize()
 
-    def add_transfer_rate(self, material, components, transfer_rate,
-                         transfer_rate_units='1/s', destination_material=None):
+    def add_transfer_rate(
+            self,
+            material: Union[str, int, Material],
+            components: Sequence[str],
+            transfer_rate: float,
+            transfer_rate_units: str = '1/s',
+            destination_material: Optional[Union[str, int, Material]] = None
+        ):
         """Add transfer rates to depletable material.
 
         Parameters
@@ -1001,9 +1031,17 @@ class SIIntegrator(Integrator):
 
     """
 
-    def __init__(self, operator, timesteps, power=None, power_density=None,
-                 source_rates=None, timestep_units='s', n_steps=10,
-                 solver="cram48"):
+    def __init__(
+            self,
+            operator: TransportOperator,
+            timesteps: Sequence[float],
+            power: Optional[Union[float, Sequence[float]]] = None,
+            power_density: Optional[Union[float, Sequence[float]]] = None,
+            source_rates: Optional[Sequence[float]] = None,
+            timestep_units: str = 's',
+            n_steps: int = 10,
+            solver: str = "cram48"
+        ):
         check_type("n_steps", n_steps, Integral)
         check_greater_than("n_steps", n_steps, 0)
         super().__init__(
@@ -1050,7 +1088,7 @@ class SIIntegrator(Integrator):
                     if self.operator.prev_res is None:
                         n, res = self._get_bos_data_from_operator(i, p, n)
                     else:
-                        n, res = self._get_bos_data_from_restart(i, p, n)
+                        n, res = self._get_bos_data_from_restart(p, n)
                 else:
                     # Pull rates, k from previous iteration w/o
                     # re-running transport
diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py
index e06cb57253..6071689554 100644
--- a/openmc/deplete/chain.py
+++ b/openmc/deplete/chain.py
@@ -7,115 +7,109 @@ loaded from an .xml file and all the nuclides are linked together.
 from io import StringIO
 from itertools import chain
 import math
-import os
 import re
 from collections import defaultdict, namedtuple
 from collections.abc import Mapping, Iterable
 from numbers import Real, Integral
 from warnings import warn
 
-from openmc.checkvalue import check_type, check_greater_than
-from openmc.data import gnds_name, zam, DataLibrary
-from openmc.exceptions import DataError
-from .nuclide import FissionYieldDistribution
-
 import lxml.etree as ET
 import scipy.sparse as sp
 
+from openmc.checkvalue import check_type, check_greater_than
+from openmc.data import gnds_name, zam
+from .nuclide import FissionYieldDistribution, Nuclide
 import openmc.data
-from openmc._xml import clean_indentation
-from .nuclide import Nuclide, DecayTuple, ReactionTuple
 
 
-# tuple of (possible MT values, (dA, dZ), secondaries) where dA is the change in
-# the mass number and dZ is the change in the atomic number
-ReactionInfo = namedtuple('ReactionInfo', ('mts', 'dadz', 'secondaries'))
+# tuple of (possible MT values, secondaries)
+ReactionInfo = namedtuple('ReactionInfo', ('mts', 'secondaries'))
 
 REACTIONS = {
-    '(n,2nd)': ReactionInfo({11}, (-3, -1), ('H2',)),
-    '(n,2n)': ReactionInfo(set(chain([16], range(875, 892))), (-1, 0), ()),
-    '(n,3n)': ReactionInfo({17}, (-2, 0), ()),
-    '(n,na)': ReactionInfo({22}, (-4, -2), ('He4',)),
-    '(n,n3a)': ReactionInfo({23}, (-12, -6), ('He4', 'He4', 'He4')),
-    '(n,2na)': ReactionInfo({24}, (-5, -2), ('He4',)),
-    '(n,3na)': ReactionInfo({25}, (-6, -2), ('He4',)),
-    '(n,np)': ReactionInfo({28}, (-1, -1), ('H1',)),
-    '(n,n2a)': ReactionInfo({29}, (-8, -4), ('He4', 'He4')),
-    '(n,2n2a)': ReactionInfo({30}, (-9, -4), ('He4', 'He4')),
-    '(n,nd)': ReactionInfo({32}, (-2, -1), ('H2',)),
-    '(n,nt)': ReactionInfo({33}, (-3, -1), ('H3',)),
-    '(n,n3He)': ReactionInfo({34}, (-3, -2), ('He3',)),
-    '(n,nd2a)': ReactionInfo({35}, (-10, -5), ('H2', 'He4', 'He4')),
-    '(n,nt2a)': ReactionInfo({36}, (-11, -5), ('H3', 'He4', 'He4')),
-    '(n,4n)': ReactionInfo({37}, (-3, 0), ()),
-    '(n,2np)': ReactionInfo({41}, (-2, -1), ('H1',)),
-    '(n,3np)': ReactionInfo({42}, (-3, -1), ('H1',)),
-    '(n,n2p)': ReactionInfo({44}, (-2, -2), ('H1', 'H1')),
-    '(n,npa)': ReactionInfo({45}, (-5, -3), ('H1', 'He4')),
-    '(n,gamma)': ReactionInfo({102}, (1, 0), ()),
-    '(n,p)': ReactionInfo(set(chain([103], range(600, 650))), (0, -1), ('H1',)),
-    '(n,d)': ReactionInfo(set(chain([104], range(650, 700))), (-1, -1), ('H2',)),
-    '(n,t)': ReactionInfo(set(chain([105], range(700, 750))), (-2, -1), ('H3',)),
-    '(n,3He)': ReactionInfo(set(chain([106], range(750, 800))), (-2, -2), ('He3',)),
-    '(n,a)': ReactionInfo(set(chain([107], range(800, 850))), (-3, -2), ('He4',)),
-    '(n,2a)': ReactionInfo({108}, (-7, -4), ('He4', 'He4')),
-    '(n,3a)': ReactionInfo({109}, (-11, -6), ('He4', 'He4', 'He4')),
-    '(n,2p)': ReactionInfo({111}, (-1, -2), ('H1', 'H1')),
-    '(n,pa)': ReactionInfo({112}, (-4, -3), ('H1', 'He4')),
-    '(n,t2a)': ReactionInfo({113}, (-10, -5), ('H3', 'He4', 'He4')),
-    '(n,d2a)': ReactionInfo({114}, (-9, -5), ('H2', 'He4', 'He4')),
-    '(n,pd)': ReactionInfo({115}, (-2, -2), ('H1', 'H2')),
-    '(n,pt)': ReactionInfo({116}, (-3, -2), ('H1', 'H3')),
-    '(n,da)': ReactionInfo({117}, (-5, -3), ('H2', 'He4')),
-    '(n,5n)': ReactionInfo({152}, (-4, 0), ()),
-    '(n,6n)': ReactionInfo({153}, (-5, 0), ()),
-    '(n,2nt)': ReactionInfo({154}, (-4, -1), ('H3',)),
-    '(n,ta)': ReactionInfo({155}, (-6, -3), ('H3', 'He4')),
-    '(n,4np)': ReactionInfo({156}, (-4, -1), ('H1',)),
-    '(n,3nd)': ReactionInfo({157}, (-4, -1), ('H2',)),
-    '(n,nda)': ReactionInfo({158}, (-6, -3), ('H2', 'He4')),
-    '(n,2npa)': ReactionInfo({159}, (-6, -3), ('H1', 'He4')),
-    '(n,7n)': ReactionInfo({160}, (-6, 0), ()),
-    '(n,8n)': ReactionInfo({161}, (-7, 0), ()),
-    '(n,5np)': ReactionInfo({162}, (-5, -1), ('H1',)),
-    '(n,6np)': ReactionInfo({163}, (-6, -1), ('H1',)),
-    '(n,7np)': ReactionInfo({164}, (-7, -1), ('H1',)),
-    '(n,4na)': ReactionInfo({165}, (-7, -2), ('He4',)),
-    '(n,5na)': ReactionInfo({166}, (-8, -2), ('He4',)),
-    '(n,6na)': ReactionInfo({167}, (-9, -2), ('He4',)),
-    '(n,7na)': ReactionInfo({168}, (-10, -2), ('He4',)),
-    '(n,4nd)': ReactionInfo({169}, (-5, -1), ('H2',)),
-    '(n,5nd)': ReactionInfo({170}, (-6, -1), ('H2',)),
-    '(n,6nd)': ReactionInfo({171}, (-7, -1), ('H2',)),
-    '(n,3nt)': ReactionInfo({172}, (-5, -1), ('H3',)),
-    '(n,4nt)': ReactionInfo({173}, (-6, -1), ('H3',)),
-    '(n,5nt)': ReactionInfo({174}, (-7, -1), ('H3',)),
-    '(n,6nt)': ReactionInfo({175}, (-8, -1), ('H3',)),
-    '(n,2n3He)': ReactionInfo({176}, (-4, -2), ('He3',)),
-    '(n,3n3He)': ReactionInfo({177}, (-5, -2), ('He3',)),
-    '(n,4n3He)': ReactionInfo({178}, (-6, -2), ('He3',)),
-    '(n,3n2p)': ReactionInfo({179}, (-4, -2), ('H1', 'H1')),
-    '(n,3n2a)': ReactionInfo({180}, (-10, -4), ('He4', 'He4')),
-    '(n,3npa)': ReactionInfo({181}, (-7, -3), ('H1', 'He4')),
-    '(n,dt)': ReactionInfo({182}, (-4, -2), ('H2', 'H3')),
-    '(n,npd)': ReactionInfo({183}, (-3, -2), ('H1', 'H2')),
-    '(n,npt)': ReactionInfo({184}, (-4, -2), ('H1', 'H3')),
-    '(n,ndt)': ReactionInfo({185}, (-5, -2), ('H2', 'H3')),
-    '(n,np3He)': ReactionInfo({186}, (-4, -3), ('H1', 'He3')),
-    '(n,nd3He)': ReactionInfo({187}, (-5, -3), ('H2', 'He3')),
-    '(n,nt3He)': ReactionInfo({188}, (-6, -3), ('H3', 'He3')),
-    '(n,nta)': ReactionInfo({189}, (-7, -3), ('H3', 'He4')),
-    '(n,2n2p)': ReactionInfo({190}, (-3, -2), ('H1', 'H1')),
-    '(n,p3He)': ReactionInfo({191}, (-4, -3), ('H1', 'He3')),
-    '(n,d3He)': ReactionInfo({192}, (-5, -3), ('H2', 'He3')),
-    '(n,3Hea)': ReactionInfo({193}, (-6, -4), ('He3', 'He4')),
-    '(n,4n2p)': ReactionInfo({194}, (-5, -2), ('H1', 'H1')),
-    '(n,4n2a)': ReactionInfo({195}, (-11, -4), ('He4', 'He4')),
-    '(n,4npa)': ReactionInfo({196}, (-8, -3), ('H1', 'He4')),
-    '(n,3p)': ReactionInfo({197}, (-2, -3), ('H1', 'H1', 'H1')),
-    '(n,n3p)': ReactionInfo({198}, (-3, -3), ('H1', 'H1', 'H1')),
-    '(n,3n2pa)': ReactionInfo({199}, (-8, -4), ('H1', 'H1', 'He4')),
-    '(n,5n2p)': ReactionInfo({200}, (-6, -2), ('H1', 'H1')),
+    '(n,2nd)': ReactionInfo({11}, ('H2',)),
+    '(n,2n)': ReactionInfo(set(chain([16], range(875, 892))), ()),
+    '(n,3n)': ReactionInfo({17}, ()),
+    '(n,na)': ReactionInfo({22}, ('He4',)),
+    '(n,n3a)': ReactionInfo({23}, ('He4', 'He4', 'He4')),
+    '(n,2na)': ReactionInfo({24}, ('He4',)),
+    '(n,3na)': ReactionInfo({25}, ('He4',)),
+    '(n,np)': ReactionInfo({28}, ('H1',)),
+    '(n,n2a)': ReactionInfo({29}, ('He4', 'He4')),
+    '(n,2n2a)': ReactionInfo({30}, ('He4', 'He4')),
+    '(n,nd)': ReactionInfo({32}, ('H2',)),
+    '(n,nt)': ReactionInfo({33}, ('H3',)),
+    '(n,n3He)': ReactionInfo({34}, ('He3',)),
+    '(n,nd2a)': ReactionInfo({35}, ('H2', 'He4', 'He4')),
+    '(n,nt2a)': ReactionInfo({36}, ('H3', 'He4', 'He4')),
+    '(n,4n)': ReactionInfo({37}, ()),
+    '(n,2np)': ReactionInfo({41}, ('H1',)),
+    '(n,3np)': ReactionInfo({42}, ('H1',)),
+    '(n,n2p)': ReactionInfo({44}, ('H1', 'H1')),
+    '(n,npa)': ReactionInfo({45}, ('H1', 'He4')),
+    '(n,gamma)': ReactionInfo({102}, ()),
+    '(n,p)': ReactionInfo(set(chain([103], range(600, 650))), ('H1',)),
+    '(n,d)': ReactionInfo(set(chain([104], range(650, 700))), ('H2',)),
+    '(n,t)': ReactionInfo(set(chain([105], range(700, 750))), ('H3',)),
+    '(n,3He)': ReactionInfo(set(chain([106], range(750, 800))), ('He3',)),
+    '(n,a)': ReactionInfo(set(chain([107], range(800, 850))), ('He4',)),
+    '(n,2a)': ReactionInfo({108}, ('He4', 'He4')),
+    '(n,3a)': ReactionInfo({109}, ('He4', 'He4', 'He4')),
+    '(n,2p)': ReactionInfo({111}, ('H1', 'H1')),
+    '(n,pa)': ReactionInfo({112}, ('H1', 'He4')),
+    '(n,t2a)': ReactionInfo({113}, ('H3', 'He4', 'He4')),
+    '(n,d2a)': ReactionInfo({114}, ('H2', 'He4', 'He4')),
+    '(n,pd)': ReactionInfo({115}, ('H1', 'H2')),
+    '(n,pt)': ReactionInfo({116}, ('H1', 'H3')),
+    '(n,da)': ReactionInfo({117}, ('H2', 'He4')),
+    '(n,5n)': ReactionInfo({152}, ()),
+    '(n,6n)': ReactionInfo({153}, ()),
+    '(n,2nt)': ReactionInfo({154}, ('H3',)),
+    '(n,ta)': ReactionInfo({155}, ('H3', 'He4')),
+    '(n,4np)': ReactionInfo({156}, ('H1',)),
+    '(n,3nd)': ReactionInfo({157}, ('H2',)),
+    '(n,nda)': ReactionInfo({158}, ('H2', 'He4')),
+    '(n,2npa)': ReactionInfo({159}, ('H1', 'He4')),
+    '(n,7n)': ReactionInfo({160}, ()),
+    '(n,8n)': ReactionInfo({161}, ()),
+    '(n,5np)': ReactionInfo({162}, ('H1',)),
+    '(n,6np)': ReactionInfo({163}, ('H1',)),
+    '(n,7np)': ReactionInfo({164}, ('H1',)),
+    '(n,4na)': ReactionInfo({165}, ('He4',)),
+    '(n,5na)': ReactionInfo({166}, ('He4',)),
+    '(n,6na)': ReactionInfo({167}, ('He4',)),
+    '(n,7na)': ReactionInfo({168}, ('He4',)),
+    '(n,4nd)': ReactionInfo({169}, ('H2',)),
+    '(n,5nd)': ReactionInfo({170}, ('H2',)),
+    '(n,6nd)': ReactionInfo({171}, ('H2',)),
+    '(n,3nt)': ReactionInfo({172}, ('H3',)),
+    '(n,4nt)': ReactionInfo({173}, ('H3',)),
+    '(n,5nt)': ReactionInfo({174}, ('H3',)),
+    '(n,6nt)': ReactionInfo({175}, ('H3',)),
+    '(n,2n3He)': ReactionInfo({176}, ('He3',)),
+    '(n,3n3He)': ReactionInfo({177}, ('He3',)),
+    '(n,4n3He)': ReactionInfo({178}, ('He3',)),
+    '(n,3n2p)': ReactionInfo({179}, ('H1', 'H1')),
+    '(n,3n2a)': ReactionInfo({180}, ('He4', 'He4')),
+    '(n,3npa)': ReactionInfo({181}, ('H1', 'He4')),
+    '(n,dt)': ReactionInfo({182}, ('H2', 'H3')),
+    '(n,npd)': ReactionInfo({183}, ('H1', 'H2')),
+    '(n,npt)': ReactionInfo({184}, ('H1', 'H3')),
+    '(n,ndt)': ReactionInfo({185}, ('H2', 'H3')),
+    '(n,np3He)': ReactionInfo({186}, ('H1', 'He3')),
+    '(n,nd3He)': ReactionInfo({187}, ('H2', 'He3')),
+    '(n,nt3He)': ReactionInfo({188}, ('H3', 'He3')),
+    '(n,nta)': ReactionInfo({189}, ('H3', 'He4')),
+    '(n,2n2p)': ReactionInfo({190}, ('H1', 'H1')),
+    '(n,p3He)': ReactionInfo({191}, ('H1', 'He3')),
+    '(n,d3He)': ReactionInfo({192}, ('H2', 'He3')),
+    '(n,3Hea)': ReactionInfo({193}, ('He3', 'He4')),
+    '(n,4n2p)': ReactionInfo({194}, ('H1', 'H1')),
+    '(n,4n2a)': ReactionInfo({195}, ('He4', 'He4')),
+    '(n,4npa)': ReactionInfo({196}, ('H1', 'He4')),
+    '(n,3p)': ReactionInfo({197}, ('H1', 'H1', 'H1')),
+    '(n,n3p)': ReactionInfo({198}, ('H1', 'H1', 'H1')),
+    '(n,3n2pa)': ReactionInfo({199}, ('H1', 'H1', 'He4')),
+    '(n,5n2p)': ReactionInfo({200}, ('H1', 'H1')),
 }
 
 __all__ = ["Chain", "REACTIONS"]
@@ -418,9 +412,9 @@ class Chain:
             if parent in reactions:
                 reactions_available = set(reactions[parent].keys())
                 for name in transmutation_reactions:
-                    mts, changes, _ = REACTIONS[name]
+                    mts = REACTIONS[name].mts
+                    delta_A, delta_Z = openmc.data.DADZ[name]
                     if mts & reactions_available:
-                        delta_A, delta_Z = changes
                         A = data.nuclide['mass_number'] + delta_A
                         Z = data.nuclide['atomic_number'] + delta_Z
                         daughter = '{}{}'.format(openmc.data.ATOMIC_SYMBOL[Z], A)
diff --git a/openmc/filter_expansion.py b/openmc/filter_expansion.py
index b9860846fd..f8e677578f 100644
--- a/openmc/filter_expansion.py
+++ b/openmc/filter_expansion.py
@@ -53,6 +53,34 @@ class ExpansionFilter(Filter):
         order = int(elem.find('order').text)
         return cls(order, filter_id=filter_id)
 
+    def merge(self, other):
+        """Merge this filter with another.
+
+        This overrides the behavior of the parent Filter class, since its
+        merging technique is to take the union of the set of bins of each
+        filter. That technique does not apply to expansion filters, since the
+        argument should be the maximum filter order rather than the list of all
+        bins.
+
+        Parameters
+        ----------
+        other : openmc.Filter
+            Filter to merge with
+
+        Returns
+        -------
+        merged_filter : openmc.Filter
+            Filter resulting from the merge
+
+        """
+
+        if not self.can_merge(other):
+            msg = f'Unable to merge "{type(self)}" with "{type(other)}"'
+            raise ValueError(msg)
+
+        # Create a new filter with these bins and a new auto-generated ID
+        return type(self)(max(self.order, other.order))
+
 
 class LegendreFilter(ExpansionFilter):
     r"""Score Legendre expansion moments up to specified order.
diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py
index 599b750fb0..7f17dd81ff 100644
--- a/openmc/lib/mesh.py
+++ b/openmc/lib/mesh.py
@@ -39,6 +39,9 @@ _dll.openmc_mesh_set_id.errcheck = _error_handler
 _dll.openmc_mesh_get_n_elements.argtypes = [c_int32, POINTER(c_size_t)]
 _dll.openmc_mesh_get_n_elements.restype = c_int
 _dll.openmc_mesh_get_n_elements.errcheck = _error_handler
+_dll.openmc_mesh_get_volumes.argtypes = [c_int32, POINTER(c_double)]
+_dll.openmc_mesh_get_volumes.restype = c_int
+_dll.openmc_mesh_get_volumes.errcheck = _error_handler
 _dll.openmc_mesh_material_volumes.argtypes = [
     c_int32, c_int, c_int, c_int, POINTER(_MaterialVolume),
     POINTER(c_int), POINTER(c_uint64)]
@@ -149,11 +152,18 @@ class Mesh(_FortranObjectWithID):
         _dll.openmc_mesh_set_id(self._index, mesh_id)
 
     @property
-    def n_elements(self):
+    def n_elements(self) -> int:
         n = c_size_t()
         _dll.openmc_mesh_get_n_elements(self._index, n)
         return n.value
 
+    @property
+    def volumes(self) -> np.ndarray:
+        volumes = np.empty((self.n_elements,))
+        _dll.openmc_mesh_get_volumes(
+            self._index, volumes.ctypes.data_as(POINTER(c_double)))
+        return volumes
+
     def material_volumes(
             self,
             n_samples: int = 10_000,
@@ -276,6 +286,10 @@ class RegularMesh(Mesh):
         are given, it is assumed that the mesh is an x-y mesh.
     width : numpy.ndarray
         The width of mesh cells in each direction.
+    n_elements : int
+        Total number of mesh elements.
+    volumes : numpy.ndarray
+        Volume of each mesh element in [cm^3]
 
     """
     mesh_type = 'regular'
@@ -358,6 +372,10 @@ class RectilinearMesh(Mesh):
         The upper-right corner of the structrued mesh.
     width : numpy.ndarray
         The width of mesh cells in each direction.
+    n_elements : int
+        Total number of mesh elements.
+    volumes : numpy.ndarray
+        Volume of each mesh element in [cm^3]
 
     """
     mesh_type = 'rectilinear'
@@ -457,6 +475,10 @@ class CylindricalMesh(Mesh):
         The upper-right corner of the structrued mesh.
     width : numpy.ndarray
         The width of mesh cells in each direction.
+    n_elements : int
+        Total number of mesh elements.
+    volumes : numpy.ndarray
+        Volume of each mesh element in [cm^3]
 
     """
     mesh_type = 'cylindrical'
@@ -555,6 +577,10 @@ class SphericalMesh(Mesh):
         The upper-right corner of the structrued mesh.
     width : numpy.ndarray
         The width of mesh cells in each direction.
+    n_elements : int
+        Total number of mesh elements.
+    volumes : numpy.ndarray
+        Volume of each mesh element in [cm^3]
 
     """
     mesh_type = 'spherical'
diff --git a/openmc/material.py b/openmc/material.py
index d42c489aee..c31862a18a 100644
--- a/openmc/material.py
+++ b/openmc/material.py
@@ -571,7 +571,7 @@ class Material(IDManagerMixin):
 
         for component, params in components.items():
             cv.check_type('component', component, str)
-            if isinstance(params, float):
+            if isinstance(params, Real):
                 params = {'percent': params}
 
             else:
diff --git a/openmc/model/model.py b/openmc/model/model.py
index 5194fec1da..2be6fcefa3 100644
--- a/openmc/model/model.py
+++ b/openmc/model/model.py
@@ -449,7 +449,7 @@ class Model:
         # Create directory if required
         d = Path(directory)
         if not d.is_dir():
-            d.mkdir(parents=True)
+            d.mkdir(parents=True, exist_ok=True)
 
         self.settings.export_to_xml(d)
         self.geometry.export_to_xml(d, remove_surfs=remove_surfs)
@@ -488,13 +488,16 @@ class Model:
         # if the provided path doesn't end with the XML extension, assume the
         # input path is meant to be a directory. If the directory does not
         # exist, create it and place a 'model.xml' file there.
-        if not str(xml_path).endswith('.xml') and not xml_path.exists():
-            os.mkdir(xml_path)
+        if not str(xml_path).endswith('.xml'):
+            if not xml_path.exists():
+                xml_path.mkdir(parents=True, exist_ok=True)
+            elif not xml_path.is_dir():
+                raise FileExistsError(f"File exists and is not a directory: '{xml_path}'")
             xml_path /= 'model.xml'
         # if this is an XML file location and the file's parent directory does
         # not exist, create it before continuing
         elif not xml_path.parent.exists():
-            os.mkdir(xml_path.parent)
+            xml_path.parent.mkdir(parents=True, exist_ok=True)
 
         if remove_surfs:
             warnings.warn("remove_surfs kwarg will be deprecated soon, please "
@@ -709,9 +712,10 @@ class Model:
                     self.export_to_model_xml(**export_kwargs)
                 else:
                     self.export_to_xml(**export_kwargs)
+                path_input = export_kwargs.get("path", None)
                 openmc.run(particles, threads, geometry_debug, restart_file,
                            tracks, output, Path('.'), openmc_exec, mpi_args,
-                           event_based)
+                           event_based, path_input)
 
             # Get output directory and return the last statepoint written
             if self.settings.output and 'path' in self.settings.output:
diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py
index e4b7020d89..8f48b2f711 100644
--- a/openmc/model/surface_composite.py
+++ b/openmc/model/surface_composite.py
@@ -617,7 +617,7 @@ class RectangularParallelepiped(CompositeSurface):
 class XConeOneSided(CompositeSurface):
     """One-sided cone parallel the x-axis
 
-    A one-sided cone is composed of a normal cone surface and an "ambiguity"
+    A one-sided cone is composed of a normal cone surface and a "disambiguation"
     surface that eliminates the ambiguity as to which region of space is
     included. This class acts as a proper surface, meaning that unary `+` and
     `-` operators applied to it will produce a half-space. The negative side is
@@ -634,7 +634,9 @@ class XConeOneSided(CompositeSurface):
     z0 : float, optional
         z-coordinate of the apex. Defaults to 0.
     r2 : float, optional
-        Parameter related to the aperature. Defaults to 1.
+        Parameter related to the aperture [:math:`\\rm cm^2`].
+        It can be interpreted as the increase in the radius squared per cm along 
+        the cone's axis of revolution.
     up : bool
         Whether to select the side of the cone that extends to infinity in the
         positive direction of the coordinate axis (the positive half-space of
@@ -647,7 +649,7 @@ class XConeOneSided(CompositeSurface):
     cone : openmc.XCone
         Regular two-sided cone
     plane : openmc.XPlane
-        Ambiguity surface
+        Disambiguation surface
     up : bool
         Whether to select the side of the cone that extends to infinity in the
         positive direction of the coordinate axis (the positive half-space of
@@ -675,7 +677,7 @@ class XConeOneSided(CompositeSurface):
 class YConeOneSided(CompositeSurface):
     """One-sided cone parallel the y-axis
 
-    A one-sided cone is composed of a normal cone surface and an "ambiguity"
+    A one-sided cone is composed of a normal cone surface and a "disambiguation"
     surface that eliminates the ambiguity as to which region of space is
     included. This class acts as a proper surface, meaning that unary `+` and
     `-` operators applied to it will produce a half-space. The negative side is
@@ -692,7 +694,9 @@ class YConeOneSided(CompositeSurface):
     z0 : float, optional
         z-coordinate of the apex. Defaults to 0.
     r2 : float, optional
-        Parameter related to the aperature. Defaults to 1.
+        Parameter related to the aperture [:math:`\\rm cm^2`].
+        It can be interpreted as the increase in the radius squared per cm along 
+        the cone's axis of revolution.
     up : bool
         Whether to select the side of the cone that extends to infinity in the
         positive direction of the coordinate axis (the positive half-space of
@@ -705,7 +709,7 @@ class YConeOneSided(CompositeSurface):
     cone : openmc.YCone
         Regular two-sided cone
     plane : openmc.YPlane
-        Ambiguity surface
+        Disambiguation surface
     up : bool
         Whether to select the side of the cone that extends to infinity in the
         positive direction of the coordinate axis (the positive half-space of
@@ -727,7 +731,7 @@ class YConeOneSided(CompositeSurface):
 class ZConeOneSided(CompositeSurface):
     """One-sided cone parallel the z-axis
 
-    A one-sided cone is composed of a normal cone surface and an "ambiguity"
+    A one-sided cone is composed of a normal cone surface and a "disambiguation"
     surface that eliminates the ambiguity as to which region of space is
     included. This class acts as a proper surface, meaning that unary `+` and
     `-` operators applied to it will produce a half-space. The negative side is
@@ -744,7 +748,9 @@ class ZConeOneSided(CompositeSurface):
     z0 : float, optional
         z-coordinate of the apex. Defaults to 0.
     r2 : float, optional
-        Parameter related to the aperature. Defaults to 1.
+        Parameter related to the aperture [:math:`\\rm cm^2`].
+        It can be interpreted as the increase in the radius squared per cm along 
+        the cone's axis of revolution.
     up : bool
         Whether to select the side of the cone that extends to infinity in the
         positive direction of the coordinate axis (the positive half-space of
@@ -757,7 +763,7 @@ class ZConeOneSided(CompositeSurface):
     cone : openmc.ZCone
         Regular two-sided cone
     plane : openmc.ZPlane
-        Ambiguity surface
+        Disambiguation surface
     up : bool
         Whether to select the side of the cone that extends to infinity in the
         positive direction of the coordinate axis (the positive half-space of
diff --git a/openmc/plotter.py b/openmc/plotter.py
index 193ba9cbe6..849672fceb 100644
--- a/openmc/plotter.py
+++ b/openmc/plotter.py
@@ -58,6 +58,11 @@ ELEMENT_NAMES = list(openmc.data.ELEMENT_SYMBOL.values())[1:]
 def _get_legend_label(this, type):
     """Gets a label for the element or nuclide or material and reaction plotted"""
     if isinstance(this, str):
+        if type in openmc.data.DADZ:
+            z, a, m = openmc.data.zam(this)
+            da, dz = openmc.data.DADZ[type]
+            gnds_name = openmc.data.gnds_name(z + dz, a + da, m)
+            return f'{this} {type} {gnds_name}'
         return f'{this} {type}'
     elif this.name == '':
         return f'Material {this.id} {type}'
@@ -100,7 +105,7 @@ def _get_title(reactions):
 def plot_xs(reactions, divisor_types=None, temperature=294., axis=None,
             sab_name=None, ce_cross_sections=None, mg_cross_sections=None,
             enrichment=None, plot_CE=True, orders=None, divisor_orders=None,
-            **kwargs):
+            energy_axis_units="eV", **kwargs):
     """Creates a figure of continuous-energy cross sections for this item.
 
     Parameters
@@ -143,6 +148,10 @@ def plot_xs(reactions, divisor_types=None, temperature=294., axis=None,
     **kwargs :
         All keyword arguments are passed to
         :func:`matplotlib.pyplot.figure`.
+    energy_axis_units : {'eV', 'keV', 'MeV'}
+        Units used on the plot energy axis
+
+        .. versionadded:: 0.14.1
 
     Returns
     -------
@@ -156,6 +165,9 @@ def plot_xs(reactions, divisor_types=None, temperature=294., axis=None,
     import matplotlib.pyplot as plt
 
     cv.check_type("plot_CE", plot_CE, bool)
+    cv.check_value("energy_axis_units", energy_axis_units, {"eV", "keV", "MeV"})
+
+    axis_scaling_factor = {"eV": 1.0, "keV": 1e-3, "MeV": 1e-6}
 
     # Generate the plot
     if axis is None:
@@ -212,6 +224,8 @@ def plot_xs(reactions, divisor_types=None, temperature=294., axis=None,
                     if divisor_types[line] != 'unity':
                         types[line] += ' / ' + divisor_types[line]
 
+        E *= axis_scaling_factor[energy_axis_units]
+
         # Plot the data
         for i in range(len(data)):
             data[i, :] = np.nan_to_num(data[i, :])
@@ -227,9 +241,12 @@ def plot_xs(reactions, divisor_types=None, temperature=294., axis=None,
         ax.set_xscale('log')
         ax.set_yscale('log')
 
-    ax.set_xlabel('Energy [eV]')
+    ax.set_xlabel(f"Energy [{energy_axis_units}]")
     if plot_CE:
-        ax.set_xlim(_MIN_E, _MAX_E)
+        ax.set_xlim(
+            _MIN_E * axis_scaling_factor[energy_axis_units],
+            _MAX_E * axis_scaling_factor[energy_axis_units],
+        )
     else:
         ax.set_xlim(E[-1], E[0])
 
diff --git a/openmc/source.py b/openmc/source.py
index 1ba90c9213..441bf4e76b 100644
--- a/openmc/source.py
+++ b/openmc/source.py
@@ -885,3 +885,37 @@ def write_source_file(
     with h5py.File(filename, **kwargs) as fh:
         fh.attrs['filetype'] = np.string_("source")
         fh.create_dataset('source_bank', data=arr, dtype=source_dtype)
+
+
+def read_source_file(filename: PathLike) -> typing.List[SourceParticle]:
+    """Read a source file and return a list of source particles.
+
+    .. versionadded:: 0.14.1
+
+    Parameters
+    ----------
+    filename : str or path-like
+        Path to source file to read
+
+    Returns
+    -------
+    list of SourceParticle
+        Source particles read from file
+
+    See Also
+    --------
+    openmc.SourceParticle
+
+    """
+    with h5py.File(filename, 'r') as fh:
+        filetype = fh.attrs['filetype']
+        arr = fh['source_bank'][...]
+
+    if filetype != b'source':
+        raise ValueError(f'File {filename} is not a source file')
+
+    source_particles = []
+    for *params, particle in arr:
+        source_particles.append(SourceParticle(*params, ParticleType(particle)))
+
+    return source_particles
diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py
index b33f88af86..e7ff64257a 100644
--- a/openmc/stats/univariate.py
+++ b/openmc/stats/univariate.py
@@ -280,6 +280,9 @@ class Discrete(Univariate):
         Discrete distribution with low-importance points removed
 
         """
+        cv.check_less_than("tolerance", tolerance, 1.0, equality=True)
+        cv.check_greater_than("tolerance", tolerance, 0.0, equality=True)
+
         # Determine (reversed) sorted order of probabilities
         intensity = self.p * self.x
         index_sort = np.argsort(intensity)[::-1]
diff --git a/openmc/surface.py b/openmc/surface.py
index a16b6712d6..bc10f7e8a4 100644
--- a/openmc/surface.py
+++ b/openmc/surface.py
@@ -1748,6 +1748,11 @@ class Sphere(QuadricMixin, Surface):
 class Cone(QuadricMixin, Surface):
     """A conical surface parallel to the x-, y-, or z-axis.
 
+    .. Note::
+        This creates a double cone, which is two one-sided cones that meet at their apex.
+        For a one-sided cone see :class:`~openmc.model.XConeOneSided`,
+        :class:`~openmc.model.YConeOneSided`, and :class:`~openmc.model.ZConeOneSided`.
+
     Parameters
     ----------
     x0 : float, optional
@@ -1757,7 +1762,9 @@ class Cone(QuadricMixin, Surface):
     z0 : float, optional
         z-coordinate of the apex in [cm]. Defaults to 0.
     r2 : float, optional
-        Parameter related to the aperature. Defaults to 1.
+        Parameter related to the aperture [:math:`\\rm cm^2`].
+        It can be interpreted as the increase in the radius squared per cm along 
+        the cone's axis of revolution.
     dx : float, optional
         x-component of the vector representing the axis of the cone.
         Defaults to 0.
@@ -1791,7 +1798,7 @@ class Cone(QuadricMixin, Surface):
     z0 : float
         z-coordinate of the apex in [cm]
     r2 : float
-        Parameter related to the aperature
+        Parameter related to the aperature [cm^2]
     dx : float
         x-component of the vector representing the axis of the cone.
     dy : float
@@ -1900,6 +1907,10 @@ class XCone(QuadricMixin, Surface):
     """A cone parallel to the x-axis of the form :math:`(y - y_0)^2 + (z - z_0)^2 =
     r^2 (x - x_0)^2`.
 
+    .. Note::
+        This creates a double cone, which is two one-sided cones that meet at their apex.
+        For a one-sided cone see :class:`~openmc.model.XConeOneSided`.
+
     Parameters
     ----------
     x0 : float, optional
@@ -1909,7 +1920,9 @@ class XCone(QuadricMixin, Surface):
     z0 : float, optional
         z-coordinate of the apex in [cm]. Defaults to 0.
     r2 : float, optional
-        Parameter related to the aperature. Defaults to 1.
+        Parameter related to the aperture [:math:`\\rm cm^2`].
+        It can be interpreted as the increase in the radius squared per cm along 
+        the cone's axis of revolution.
     boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional
         Boundary condition that defines the behavior for particles hitting the
         surface. Defaults to transmissive boundary condition where particles
@@ -1995,6 +2008,10 @@ class YCone(QuadricMixin, Surface):
     """A cone parallel to the y-axis of the form :math:`(x - x_0)^2 + (z - z_0)^2 =
     r^2 (y - y_0)^2`.
 
+    .. Note::
+        This creates a double cone, which is two one-sided cones that meet at their apex.
+        For a one-sided cone see :class:`~openmc.model.YConeOneSided`.
+
     Parameters
     ----------
     x0 : float, optional
@@ -2004,7 +2021,9 @@ class YCone(QuadricMixin, Surface):
     z0 : float, optional
         z-coordinate of the apex in [cm]. Defaults to 0.
     r2 : float, optional
-        Parameter related to the aperature. Defaults to 1.
+        Parameter related to the aperture [:math:`\\rm cm^2`].
+        It can be interpreted as the increase in the radius squared per cm along 
+        the cone's axis of revolution.
     boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional
         Boundary condition that defines the behavior for particles hitting the
         surface. Defaults to transmissive boundary condition where particles
@@ -2090,6 +2109,10 @@ class ZCone(QuadricMixin, Surface):
     """A cone parallel to the z-axis of the form :math:`(x - x_0)^2 + (y - y_0)^2 =
     r^2 (z - z_0)^2`.
 
+    .. Note::
+        This creates a double cone, which is two one-sided cones that meet at their apex.
+        For a one-sided cone see :class:`~openmc.model.ZConeOneSided`.
+
     Parameters
     ----------
     x0 : float, optional
@@ -2099,7 +2122,9 @@ class ZCone(QuadricMixin, Surface):
     z0 : float, optional
         z-coordinate of the apex in [cm]. Defaults to 0.
     r2 : float, optional
-        Parameter related to the aperature. Defaults to 1.
+        Parameter related to the aperature [cm^2].
+        This is the square of the radius of the cone 1 cm from.
+        This can also be treated as the square of the slope of the cone relative to its axis.
     boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional
         Boundary condition that defines the behavior for particles hitting the
         surface. Defaults to transmissive boundary condition where particles
diff --git a/src/cell.cpp b/src/cell.cpp
index 00083d967b..a46f05687e 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -1272,6 +1272,9 @@ struct ParentCellStack {
   //! compute an instance for the provided distribcell index
   int32_t compute_instance(int32_t distribcell_index) const
   {
+    if (distribcell_index == C_NONE)
+      return 0;
+
     int32_t instance = 0;
     for (const auto& parent_cell : this->parent_cells_) {
       auto& cell = model::cells[parent_cell.cell_index];
diff --git a/src/mesh.cpp b/src/mesh.cpp
index f343b27e33..bc6781ad31 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -1881,6 +1881,17 @@ extern "C" int openmc_mesh_get_n_elements(int32_t index, size_t* n)
   return 0;
 }
 
+//! Get the volume of each element in the mesh
+extern "C" int openmc_mesh_get_volumes(int32_t index, double* volumes)
+{
+  if (int err = check_mesh(index))
+    return err;
+  for (int i = 0; i < model::meshes[index]->n_bins(); ++i) {
+    volumes[i] = model::meshes[index]->volume(i);
+  }
+  return 0;
+}
+
 extern "C" int openmc_mesh_material_volumes(int32_t index, int n_sample,
   int bin, int result_size, void* result, int* hits, uint64_t* seed)
 {
diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py
index e4e0231355..9753a53f73 100644
--- a/tests/unit_tests/test_deplete_chain.py
+++ b/tests/unit_tests/test_deplete_chain.py
@@ -5,6 +5,7 @@ from itertools import product
 from math import log
 import os
 from pathlib import Path
+import warnings
 
 import numpy as np
 from openmc.mpi import comm
@@ -438,9 +439,9 @@ def test_validate(simple_chain):
     simple_chain["C"].yield_data = {0.0253: {"A": 1.4, "B": 0.6}}
 
     assert simple_chain.validate(strict=True, tolerance=0.0)
-    with pytest.warns(None) as record:
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
         assert simple_chain.validate(strict=False, quiet=False, tolerance=0.0)
-    assert len(record) == 0
 
     # Mess up "earlier" nuclide's reactions
     decay_mode = simple_chain["A"].decay_modes.pop()
diff --git a/tests/unit_tests/test_deplete_nuclide.py b/tests/unit_tests/test_deplete_nuclide.py
index 6482d76a59..f2bb7d1b65 100644
--- a/tests/unit_tests/test_deplete_nuclide.py
+++ b/tests/unit_tests/test_deplete_nuclide.py
@@ -1,6 +1,7 @@
 """Tests for the openmc.deplete.Nuclide class."""
 
 import copy
+import warnings
 
 import lxml.etree as ET
 import numpy as np
@@ -277,9 +278,9 @@ def test_validate():
     }
 
     # nuclide is good and should have no warnings raise
-    with pytest.warns(None) as record:
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
         assert nuc.validate(strict=True, quiet=False, tolerance=0.0)
-    assert len(record) == 0
 
     # invalidate decay modes
     decay = nuc.decay_modes.pop()
diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py
index 6ae13ef76e..1310a6a7fa 100644
--- a/tests/unit_tests/test_lib.py
+++ b/tests/unit_tests/test_lib.py
@@ -566,6 +566,8 @@ def test_regular_mesh(lib_init):
     assert mesh.upper_right == pytest.approx(ur)
     assert mesh.width == pytest.approx(width)
 
+    np.testing.assert_allclose(mesh.volumes, 1.0)
+
     meshes = openmc.lib.meshes
     assert isinstance(meshes, Mapping)
     assert len(meshes) == 1
@@ -644,6 +646,8 @@ def test_rectilinear_mesh(lib_init):
             for k, diff_z in enumerate(np.diff(z_grid)):
                 assert np.all(mesh.width[i, j, k, :] == (10, 10, 10))
 
+    np.testing.assert_allclose(mesh.volumes, 1000.0)
+
     with pytest.raises(exc.AllocationError):
         mesh2 = openmc.lib.RectilinearMesh(mesh.id)
 
@@ -688,6 +692,9 @@ def test_cylindrical_mesh(lib_init):
             for k, _ in enumerate(np.diff(z_grid)):
                 assert np.allclose(mesh.width[i, j, k, :], (5, deg2rad(10), 10))
 
+    np.testing.assert_allclose(mesh.volumes[::2], 10/360 * pi * 5**2 * 10)
+    np.testing.assert_allclose(mesh.volumes[1::2], 10/360 * pi * (10**2 - 5**2) * 10)
+
     with pytest.raises(exc.AllocationError):
         mesh2 = openmc.lib.CylindricalMesh(mesh.id)
 
@@ -734,6 +741,13 @@ def test_spherical_mesh(lib_init):
             for k, _ in enumerate(np.diff(phi_grid)):
                 assert np.allclose(mesh.width[i, j, k, :], (5, deg2rad(10), deg2rad(10)))
 
+    dtheta = lambda d1, d2: np.cos(deg2rad(d1)) - np.cos(deg2rad(d2))
+    f = 1/3 * deg2rad(10.)
+    np.testing.assert_allclose(mesh.volumes[::4],  f * 5**3 * dtheta(0., 10.))
+    np.testing.assert_allclose(mesh.volumes[1::4], f * (10**3 - 5**3) * dtheta(0., 10.))
+    np.testing.assert_allclose(mesh.volumes[2::4], f * 5**3 * dtheta(10., 20.))
+    np.testing.assert_allclose(mesh.volumes[3::4], f * (10**3 - 5**3) * dtheta(10., 20.))
+
     with pytest.raises(exc.AllocationError):
         mesh2 = openmc.lib.SphericalMesh(mesh.id)
 
diff --git a/tests/unit_tests/test_model.py b/tests/unit_tests/test_model.py
index 24f136c8c4..16fa18d453 100644
--- a/tests/unit_tests/test_model.py
+++ b/tests/unit_tests/test_model.py
@@ -585,3 +585,9 @@ def test_single_xml_exec(run_in_tmpdir):
 
     with pytest.raises(RuntimeError, match='input_dir'):
         openmc.run(path_input='input_dir/pincell.xml')
+
+    # Make sure path can be specified with run
+    pincell_model.run(path='my_model.xml')
+
+    os.mkdir('subdir')
+    pincell_model.run(path='subdir')
diff --git a/tests/unit_tests/test_plotter.py b/tests/unit_tests/test_plotter.py
index 8a1147b4af..fd635d89d5 100644
--- a/tests/unit_tests/test_plotter.py
+++ b/tests/unit_tests/test_plotter.py
@@ -77,10 +77,19 @@ def test_plot_xs(this):
     from matplotlib.figure import Figure
     assert isinstance(openmc.plot_xs({this: ['total', 'elastic']}), Figure)
 
+
 def test_plot_xs_mat(test_mat):
     from matplotlib.figure import Figure
     assert isinstance(openmc.plot_xs({test_mat: ['total']}), Figure)
 
+
+@pytest.mark.parametrize("units", ["eV", "keV", "MeV"])
+def test_plot_xs_energy_axis(units):
+    plot = openmc.plot_xs({'Be9': ['(n,2n)']}, energy_axis_units=units)
+    axis_text = plot.get_axes()[0].get_xaxis().get_label().get_text()
+    assert axis_text == f'Energy [{units}]'
+
+
 def test_plot_axes_labels():
     # just nuclides
     axis_label = openmc.plotter._get_yaxis_label(
@@ -156,4 +165,4 @@ def test_get_title():
     mat1.set_density('g/cm3', 1)
     mat1.name = 'my_mat'
     title = openmc.plotter._get_title(reactions={mat1: [205]})
-    assert title == 'Cross Section Plot For my_mat'
\ No newline at end of file
+    assert title == 'Cross Section Plot For my_mat'
diff --git a/tests/unit_tests/test_source_file.py b/tests/unit_tests/test_source_file.py
index ec10260baf..d9fb3c1f90 100644
--- a/tests/unit_tests/test_source_file.py
+++ b/tests/unit_tests/test_source_file.py
@@ -45,6 +45,30 @@ def test_source_file(run_in_tmpdir):
     assert np.all(arr['particle'] == 0)
 
 
+    # Ensure sites read in are consistent
+    sites = openmc.read_source_file('test_source.h5')
+
+    assert filetype == b'source'
+    xs = np.array([site.r[0] for site in sites])
+    ys = np.array([site.r[1] for site in sites])
+    zs = np.array([site.r[2] for site in sites])
+    assert np.all((xs > 0.0) & (xs < 1.0))
+    assert np.all(ys == np.arange(1000))
+    assert np.all(zs == 0.0)
+    u = np.array([s.u for s in sites])
+    assert np.all(u[..., 0] == 0.0)
+    assert np.all(u[..., 1] == 0.0)
+    assert np.all(u[..., 2] == 1.0)
+    E = np.array([s.E for s in sites])
+    assert np.all(E == n - np.arange(n))
+    wgt = np.array([s.wgt for s in sites])
+    assert np.all(wgt == 1.0)
+    dgs = np.array([s.delayed_group for s in sites])
+    assert np.all(dgs == 0)
+    p_types = np.array([s.particle for s in sites])
+    assert np.all(p_types == 0)
+
+
 def test_wrong_source_attributes(run_in_tmpdir):
     # Create a source file with animal attributes
     source_dtype = np.dtype([
diff --git a/tests/unit_tests/test_stats.py b/tests/unit_tests/test_stats.py
index 8438535ff8..f5368e912f 100644
--- a/tests/unit_tests/test_stats.py
+++ b/tests/unit_tests/test_stats.py
@@ -89,6 +89,12 @@ def test_clip_discrete():
     d_same = d.clip(1e-6, inplace=True)
     assert d_same is d
 
+    with pytest.raises(ValueError):
+        d.clip(-1.)
+    
+    with pytest.raises(ValueError):
+        d.clip(5)
+
 
 def test_uniform():
     a, b = 10.0, 20.0