From f4f1757e7900ed3eae5889a362fe78c1dd38d0f6 Mon Sep 17 00:00:00 2001
From: Perry <yrrepy@gmail.com>
Date: Sun, 29 Mar 2026 10:48:43 -0700
Subject: [PATCH 1/8] Add CRAM substeps with LU decomposition reuse
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

  Splitting a CRAM step into identical substeps and reusing LU
  factorizations improves accuracy for nuclides with large λ·dt
  while adding negligible computational cost.

  This particularly improves accuracy for nuclides approaching zero,
  where a single-step CRAM can undershoot to negative atom densities,
  and then oscillate about zero.

  - IPFCramSolver accepts substeps parameter (default 1)
  - substeps=1 uses original spsolve path (bitwise identical)
  - substeps>1 pre-computes LU via splu, reuses across substeps
  - Integrator/SIIntegrator accept substeps parameter
  - 4 new unit tests: default, bitwise match, two-half-steps,
    CRAM16 self-convergence

  Reference: Isotalo & Pusa, "Improving the Accuracy of the
  Chebyshev Rational Approximation Method Using Substeps,"
  Nucl. Sci. Eng., 183:1, 65-77 (2016).
---
 openmc/deplete/abc.py                 | 38 +++++++++++++--
 openmc/deplete/cram.py                | 42 +++++++++++++---
 tests/unit_tests/test_deplete_cram.py | 70 ++++++++++++++++++++++++++-
 3 files changed, 137 insertions(+), 13 deletions(-)

diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py
index 056f7c2737a..2cf1d335d5e 100644
--- a/openmc/deplete/abc.py
+++ b/openmc/deplete/abc.py
@@ -572,6 +572,14 @@ class Integrator(ABC):
         :attr:`solver`.
 
         .. versionadded:: 0.12
+    substeps : int, optional
+        Number of substeps per depletion interval. When greater than 1,
+        each interval is subdivided into `substeps` identical sub-intervals
+        and LU factorizations are reused across them, improving accuracy
+        for nuclides with large decay-constant × timestep products.
+        Only used when `solver` is ``cram48"`` or ``cram16``.
+
+        .. versionadded:: 0.15.3
     continue_timesteps : bool, optional
         Whether or not to treat the current solve as a continuation of a
         previous simulation. Defaults to `False`. When `False`, the depletion
@@ -631,6 +639,7 @@ def __init__(
             source_rates: Optional[Union[float, Sequence[float]]] = None,
             timestep_units: str = 's',
             solver: str = "cram48",
+            substeps: int = 1,
             continue_timesteps: bool = False,
         ):
         if continue_timesteps and operator.prev_res is None:
@@ -690,15 +699,23 @@ def __init__(
         if isinstance(solver, str):
             # Delay importing of cram module, which requires this file
             if solver == "cram48":
-                from .cram import CRAM48
-                self._solver = CRAM48
+                from .cram import CRAM48 as _default, Cram48Solver as _base
             elif solver == "cram16":
-                from .cram import CRAM16
-                self._solver = CRAM16
+                from .cram import CRAM16 as _default, Cram16Solver as _base
             else:
                 raise ValueError(
                     f"Solver {solver} not understood. Expected 'cram48' or 'cram16'")
+
+            if substeps > 1:
+                from .cram import IPFCramSolver
+                self._solver = IPFCramSolver(
+                    _base.alpha, _base.theta, _base.alpha0,
+                    substeps=substeps)
+            else:
+                self._solver = _default
         else:
+            if substeps > 1:
+                warn("substeps is ignored when a custom solver is provided")
             self.solver = solver
 
     @property
@@ -1104,6 +1121,14 @@ class SIIntegrator(Integrator):
         :attr:`solver`.
 
         .. versionadded:: 0.12
+    substeps : int, optional
+        Number of substeps per depletion interval. When greater than 1,
+        each interval is subdivided into `substeps` identical sub-intervals
+        and LU factorizations are reused across them, improving accuracy
+        for nuclides with large decay-constant × timestep products.
+        Only used when `solver` is ``cram48`` or ``cram16``.
+
+        .. versionadded:: 0.15.3
     continue_timesteps : bool, optional
         Whether or not to treat the current solve as a continuation of a
         previous simulation. Defaults to `False`. If `False`, all time
@@ -1158,13 +1183,16 @@ def __init__(
             timestep_units: str = 's',
             n_steps: int = 10,
             solver: str = "cram48",
+            substeps: int = 1,
             continue_timesteps: bool = False,
         ):
         check_type("n_steps", n_steps, Integral)
         check_greater_than("n_steps", n_steps, 0)
         super().__init__(
             operator, timesteps, power, power_density, source_rates,
-            timestep_units=timestep_units, solver=solver, continue_timesteps=continue_timesteps)
+            timestep_units=timestep_units, solver=solver,
+            substeps=substeps,
+            continue_timesteps=continue_timesteps)
         self.n_steps = n_steps
 
     def _get_bos_data_from_operator(self, step_index, step_power, n_bos):
diff --git a/openmc/deplete/cram.py b/openmc/deplete/cram.py
index cecc388f4c3..5492a27ea20 100644
--- a/openmc/deplete/cram.py
+++ b/openmc/deplete/cram.py
@@ -8,7 +8,7 @@
 import numpy as np
 import scipy.sparse.linalg as sla
 
-from openmc.checkvalue import check_type, check_length
+from openmc.checkvalue import check_type, check_length, check_greater_than
 from .abc import DepSystemSolver
 from .._sparse_compat import csc_array, eye_array
 
@@ -24,6 +24,12 @@ class IPFCramSolver(DepSystemSolver):
     Chebyshev Rational Approximation Method and Application to Burnup Equations
     <https://doi.org/10.13182/NSE15-26>`_," Nucl. Sci. Eng., 182:3, 297-318.
 
+    When `substeps` > 1, the time interval is split into `substeps` identical
+    sub-intervals and LU factorizations are reused across them, as described
+    in: A. Isotalo and M. Pusa, "`Improving the Accuracy of the Chebyshev
+    Rational Approximation Method Using Substeps
+    <https://doi.org/10.13182/NSE15-67>`_," Nucl. Sci. Eng., 183:1, 65-77.
+
     Parameters
     ----------
     alpha : numpy.ndarray
@@ -33,6 +39,8 @@ class IPFCramSolver(DepSystemSolver):
         Complex poles. Must have an equal size as ``alpha``.
     alpha0 : float
         Limit of the approximation at infinity
+    substeps : int, optional
+        Number of substeps for LU-reuse CRAM.
 
     Attributes
     ----------
@@ -43,17 +51,22 @@ class IPFCramSolver(DepSystemSolver):
         Complex poles :math:`\theta` of the rational approximation
     alpha0 : float
         Limit of the approximation at infinity
+    substeps : int
+        Number of substeps per depletion interval
 
     """
 
-    def __init__(self, alpha, theta, alpha0):
+    def __init__(self, alpha, theta, alpha0, substeps=1):
         check_type("alpha", alpha, np.ndarray, numbers.Complex)
         check_type("theta", theta, np.ndarray, numbers.Complex)
         check_length("theta", theta, alpha.size)
         check_type("alpha0", alpha0, numbers.Real)
+        check_type("substeps", substeps, numbers.Integral)
+        check_greater_than("substeps", substeps, 0)
         self.alpha = alpha
         self.theta = theta
         self.alpha0 = alpha0
+        self.substeps = substeps
 
     def __call__(self, A, n0, dt):
         """Solve depletion equations using IPF CRAM
@@ -75,12 +88,27 @@ def __call__(self, A, n0, dt):
             Final compositions after ``dt``
 
         """
-        A = dt * csc_array(A, dtype=np.float64)
+        if self.substeps == 1:
+            A = dt * csc_array(A, dtype=np.float64)
+            y = n0.copy()
+            ident = eye_array(A.shape[0], format='csc')
+            for alpha, theta in zip(self.alpha, self.theta):
+                y += 2*np.real(alpha*sla.spsolve(A - theta*ident, y))
+            return y * self.alpha0
+
+        # Substep path: pre-compute LU factorizations, reuse across substeps
+        sub_dt = dt / self.substeps
+        A_sub = sub_dt * csc_array(A, dtype=np.float64)
+        ident = eye_array(A_sub.shape[0], format='csc')
+        lu_solvers = [sla.splu(A_sub - theta * ident)
+                      for theta in self.theta]
+
         y = n0.copy()
-        ident = eye_array(A.shape[0], format='csc')
-        for alpha, theta in zip(self.alpha, self.theta):
-            y += 2*np.real(alpha*sla.spsolve(A - theta*ident, y))
-        return y * self.alpha0
+        for _ in range(self.substeps):
+            for alpha, lu in zip(self.alpha, lu_solvers):
+                y += 2 * np.real(alpha * lu.solve(y))
+            y *= self.alpha0
+        return y
 
 
 # Coefficients for IPF Cram 16
diff --git a/tests/unit_tests/test_deplete_cram.py b/tests/unit_tests/test_deplete_cram.py
index 8987fbd7aa6..ec7b8add584 100644
--- a/tests/unit_tests/test_deplete_cram.py
+++ b/tests/unit_tests/test_deplete_cram.py
@@ -1,12 +1,14 @@
 """ Tests for cram.py
 
 Compares a few Mathematica matrix exponentials to CRAM16/CRAM48.
+Tests substep accuracy against self-converged reference solutions.
 """
 
 from pytest import approx
 import numpy as np
 import scipy.sparse as sp
-from openmc.deplete.cram import CRAM16, CRAM48
+from openmc.deplete.cram import (CRAM16, CRAM48, Cram16Solver, Cram48Solver,
+                                  IPFCramSolver)
 
 
 def test_CRAM16():
@@ -35,3 +37,69 @@ def test_CRAM48():
     z0 = np.array((0.904837418035960, 0.576799023327476))
 
     assert z == approx(z0)
+
+
+# --- Substep tests ---
+
+def test_substeps_default():
+    """Default substeps=1 on module-level solvers."""
+    assert Cram48Solver.substeps == 1
+    assert Cram16Solver.substeps == 1
+
+
+def test_substeps1_matches_original():
+    """substeps=1 must be bitwise identical to original spsolve path."""
+    x = np.array([1.0, 1.0])
+    mat = sp.csr_matrix([[-1.0, 0.0], [-2.0, -3.0]])
+    dt = 0.1
+
+    z_orig = CRAM48(mat, x, dt)
+    solver_1 = IPFCramSolver(Cram48Solver.alpha, Cram48Solver.theta,
+                              Cram48Solver.alpha0, substeps=1)
+    z_sub1 = solver_1(mat, x, dt)
+
+    np.testing.assert_array_equal(z_sub1, z_orig)
+
+
+def test_substeps2_matches_two_half_steps():
+    """substeps=2 must match two independent CRAM calls with dt/2."""
+    x = np.array([1.0, 1.0])
+    mat = sp.csr_matrix([[-1.0, 0.0], [-2.0, -3.0]])
+    dt = 1.0
+
+    # Two manual half-steps using original spsolve path
+    z_half = CRAM48(mat, x, dt / 2)
+    z_two = CRAM48(mat, z_half, dt / 2)
+
+    # Single call with substeps=2
+    solver_2 = IPFCramSolver(Cram48Solver.alpha, Cram48Solver.theta,
+                              Cram48Solver.alpha0, substeps=2)
+    z_sub2 = solver_2(mat, x, dt)
+
+    assert z_sub2 == approx(z_two, rel=1e-12)
+
+
+def test_substeps_self_convergence():
+    """Increasing substeps converges toward reference solution.
+
+    Uses CRAM16 (alpha0 ~ 2e-16) where substep convergence is visible.
+    CRAM48 (alpha0 ~ 2e-47) is already near machine precision for small
+    systems; its correctness is verified by the other substep tests.
+    """
+    mat = sp.csr_matrix([[-1.0, 0.0], [-2.0, -3.0]])
+    x = np.array([1.0, 1.0])
+    dt = 50  # lambda*dt = 50 and 150, stresses CRAM16
+
+    ref_solver = IPFCramSolver(Cram16Solver.alpha, Cram16Solver.theta,
+                                Cram16Solver.alpha0, substeps=128)
+    n_ref = ref_solver(mat, x, dt)
+
+    prev_err = np.inf
+    for s in [1, 2, 4, 8, 16]:
+        solver = IPFCramSolver(Cram16Solver.alpha, Cram16Solver.theta,
+                                Cram16Solver.alpha0, substeps=s)
+        n_s = solver(mat, x, dt)
+        err = np.linalg.norm(n_s - n_ref) / np.linalg.norm(n_ref)
+        assert err < prev_err, \
+            f"substeps={s} error {err:.2e} not less than previous {prev_err:.2e}"
+        prev_err = err

From 7f0823d5fd5e5943005b455cfc30307421f893b9 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Tue, 21 Apr 2026 16:13:54 -0500
Subject: [PATCH 2/8] Add setter for IPFCramSolver.substeps

---
 openmc/deplete/cram.py | 12 ++++++++++--
 1 file changed, 10 insertions(+), 2 deletions(-)

diff --git a/openmc/deplete/cram.py b/openmc/deplete/cram.py
index 5492a27ea20..a40dd7d6a68 100644
--- a/openmc/deplete/cram.py
+++ b/openmc/deplete/cram.py
@@ -61,13 +61,21 @@ def __init__(self, alpha, theta, alpha0, substeps=1):
         check_type("theta", theta, np.ndarray, numbers.Complex)
         check_length("theta", theta, alpha.size)
         check_type("alpha0", alpha0, numbers.Real)
-        check_type("substeps", substeps, numbers.Integral)
-        check_greater_than("substeps", substeps, 0)
         self.alpha = alpha
         self.theta = theta
         self.alpha0 = alpha0
         self.substeps = substeps
 
+    @property
+    def substeps(self):
+        return self._substeps
+
+    @substeps.setter
+    def substeps(self, value):
+        check_type("substeps", value, numbers.Integral)
+        check_greater_than("substeps", value, 0)
+        self._substeps = value
+
     def __call__(self, A, n0, dt):
         """Solve depletion equations using IPF CRAM
 

From b01a9d713932e7327a572e506d9ceb29f566ae36 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Wed, 22 Apr 2026 12:17:20 -0500
Subject: [PATCH 3/8] Update docstrings

---
 openmc/deplete/abc.py  | 8 ++++----
 openmc/deplete/cram.py | 6 ++++--
 2 files changed, 8 insertions(+), 6 deletions(-)

diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py
index 2cf1d335d5e..b97136fde70 100644
--- a/openmc/deplete/abc.py
+++ b/openmc/deplete/abc.py
@@ -577,9 +577,9 @@ class Integrator(ABC):
         each interval is subdivided into `substeps` identical sub-intervals
         and LU factorizations are reused across them, improving accuracy
         for nuclides with large decay-constant × timestep products.
-        Only used when `solver` is ``cram48"`` or ``cram16``.
+        Only used when `solver` is "cram48" or "cram16".
 
-        .. versionadded:: 0.15.3
+        .. versionadded:: 0.15.4
     continue_timesteps : bool, optional
         Whether or not to treat the current solve as a continuation of a
         previous simulation. Defaults to `False`. When `False`, the depletion
@@ -1126,9 +1126,9 @@ class SIIntegrator(Integrator):
         each interval is subdivided into `substeps` identical sub-intervals
         and LU factorizations are reused across them, improving accuracy
         for nuclides with large decay-constant × timestep products.
-        Only used when `solver` is ``cram48`` or ``cram16``.
+        Only used when `solver` is "cram48" or "cram16".
 
-        .. versionadded:: 0.15.3
+        .. versionadded:: 0.15.4
     continue_timesteps : bool, optional
         Whether or not to treat the current solve as a continuation of a
         previous simulation. Defaults to `False`. If `False`, all time
diff --git a/openmc/deplete/cram.py b/openmc/deplete/cram.py
index a40dd7d6a68..a5a4a1f1e1b 100644
--- a/openmc/deplete/cram.py
+++ b/openmc/deplete/cram.py
@@ -40,7 +40,9 @@ class IPFCramSolver(DepSystemSolver):
     alpha0 : float
         Limit of the approximation at infinity
     substeps : int, optional
-        Number of substeps for LU-reuse CRAM.
+        Number of substeps per depletion interval.
+
+        .. versionadded:: 0.15.4
 
     Attributes
     ----------
@@ -56,7 +58,7 @@ class IPFCramSolver(DepSystemSolver):
 
     """
 
-    def __init__(self, alpha, theta, alpha0, substeps=1):
+    def __init__(self, alpha, theta, alpha0, substeps: int = 1):
         check_type("alpha", alpha, np.ndarray, numbers.Complex)
         check_type("theta", theta, np.ndarray, numbers.Complex)
         check_length("theta", theta, alpha.size)

From e138de3e6d954d223c57ccddf07b375aecd32ef6 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Wed, 22 Apr 2026 13:40:18 -0500
Subject: [PATCH 4/8] Pass substeps through __call__

---
 openmc/deplete/abc.py                       | 69 ++++++++++++--------
 openmc/deplete/cram.py                      | 32 +++-------
 openmc/deplete/pool.py                      | 19 ++++--
 tests/unit_tests/test_deplete_cram.py       | 39 +++++------
 tests/unit_tests/test_deplete_integrator.py | 71 ++++++++++++++++++++-
 5 files changed, 151 insertions(+), 79 deletions(-)

diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py
index b97136fde70..1bfd7b1f16b 100644
--- a/openmc/deplete/abc.py
+++ b/openmc/deplete/abc.py
@@ -575,9 +575,8 @@ class Integrator(ABC):
     substeps : int, optional
         Number of substeps per depletion interval. When greater than 1,
         each interval is subdivided into `substeps` identical sub-intervals
-        and LU factorizations are reused across them, improving accuracy
+        and LU factorizations may be reused across them, improving accuracy
         for nuclides with large decay-constant × timestep products.
-        Only used when `solver` is "cram48" or "cram16".
 
         .. versionadded:: 0.15.4
     continue_timesteps : bool, optional
@@ -607,14 +606,19 @@ class Integrator(ABC):
         Function that will solve the Bateman equations
         :math:`\frac{\partial}{\partial t}\vec{n} = A_i\vec{n}_i` with a step
         size :math:`t_i`. Can be configured using the ``solver`` argument.
-        User-supplied functions are expected to have the following signature:
-        ``solver(A, n0, t) -> n1`` where
+        User-supplied functions may have either of the following signatures:
+        ``solver(A, n0, t) -> n1`` for the default single-step behavior or
+        ``solver(A, n0, t, substeps=1) -> n1`` to support non-default
+        substeps. When :attr:`substeps` is greater than one, OpenMC will call
+        the solver with the fourth argument and unsupported custom solvers will
+        raise at runtime. In either case,
 
             * ``A`` is a :class:`scipy.sparse.csc_array` making up the
               depletion matrix
             * ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions
               for a given material in atoms/cm3
             * ``t`` is a float of the time step size in seconds, and
+            * ``substeps`` is an optional integer number of substeps, and
             * ``n1`` is a :class:`numpy.ndarray` of compositions at the
               next time step. Expected to be of the same shape as ``n0``
 
@@ -661,6 +665,8 @@ def __init__(
         # Normalize timesteps and source rates
         seconds, source_rates = _normalize_timesteps(
             timesteps, source_rates, timestep_units, operator)
+        check_type("substeps", substeps, Integral)
+        check_greater_than("substeps", substeps, 0)
 
         if continue_timesteps:
             # Get timesteps and source rates from previous results
@@ -692,6 +698,7 @@ def __init__(
 
         self.timesteps = np.asarray(seconds)
         self.source_rates = np.asarray(source_rates)
+        self.substeps = substeps
 
         self.transfer_rates = None
         self.external_source_rates = None
@@ -699,23 +706,14 @@ def __init__(
         if isinstance(solver, str):
             # Delay importing of cram module, which requires this file
             if solver == "cram48":
-                from .cram import CRAM48 as _default, Cram48Solver as _base
+                from .cram import CRAM48 as _default
             elif solver == "cram16":
-                from .cram import CRAM16 as _default, Cram16Solver as _base
+                from .cram import CRAM16 as _default
             else:
                 raise ValueError(
                     f"Solver {solver} not understood. Expected 'cram48' or 'cram16'")
-
-            if substeps > 1:
-                from .cram import IPFCramSolver
-                self._solver = IPFCramSolver(
-                    _base.alpha, _base.theta, _base.alpha0,
-                    substeps=substeps)
-            else:
-                self._solver = _default
+            self._solver = _default
         else:
-            if substeps > 1:
-                warn("substeps is ignored when a custom solver is provided")
             self.solver = solver
 
     @property
@@ -736,23 +734,32 @@ def solver(self, func):
             self._solver = func
             return
 
-        # Inspect arguments
-        if len(sig.parameters) != 3:
-            raise ValueError("Function {} does not support three arguments: "
-                             "{!s}".format(func, sig))
+        params = list(sig.parameters.values())
 
-        for ix, param in enumerate(sig.parameters.values()):
-            if param.kind in {param.KEYWORD_ONLY, param.VAR_KEYWORD}:
+        # Inspect arguments
+        if len(params) not in {3, 4}:
+            raise ValueError(
+                "Function {} must support either three arguments "
+                "(A, n0, t) or four arguments (A, n0, t, substeps=1): {!s}"
+                .format(func, sig))
+
+        for ix, param in enumerate(params):
+            if param.kind in {param.KEYWORD_ONLY, param.VAR_KEYWORD,
+                              param.VAR_POSITIONAL}:
                 raise ValueError(
                     f"Keyword arguments like {ix} at position {param} are not allowed")
 
+        if len(params) == 4 and params[3].default != 1:
+            raise ValueError(
+                f"Fourth solver argument must default to 1, not {params[3].default}")
+
         self._solver = func
 
     def _timed_deplete(self, n, rates, dt, i=None, matrix_func=None):
         start = time.time()
         results = deplete(
             self._solver, self.chain, n, rates, dt, i, matrix_func,
-            self.transfer_rates, self.external_source_rates)
+            self.transfer_rates, self.external_source_rates, self.substeps)
         return time.time() - start, results
 
     @abstractmethod
@@ -1124,9 +1131,8 @@ class SIIntegrator(Integrator):
     substeps : int, optional
         Number of substeps per depletion interval. When greater than 1,
         each interval is subdivided into `substeps` identical sub-intervals
-        and LU factorizations are reused across them, improving accuracy
+        and LU factorizations may be reused across them, improving accuracy
         for nuclides with large decay-constant × timestep products.
-        Only used when `solver` is "cram48" or "cram16".
 
         .. versionadded:: 0.15.4
     continue_timesteps : bool, optional
@@ -1158,14 +1164,19 @@ class SIIntegrator(Integrator):
         Function that will solve the Bateman equations
         :math:`\frac{\partial}{\partial t}\vec{n} = A_i\vec{n}_i` with a step
         size :math:`t_i`. Can be configured using the ``solver`` argument.
-        User-supplied functions are expected to have the following signature:
-        ``solver(A, n0, t) -> n1`` where
+        User-supplied functions may have either of the following signatures:
+        ``solver(A, n0, t) -> n1`` for the default single-step behavior or
+        ``solver(A, n0, t, substeps=1) -> n1`` to support non-default
+        substeps. When :attr:`substeps` is greater than one, OpenMC will call
+        the solver with the fourth argument and unsupported custom solvers will
+        raise at runtime. In either case,
 
             * ``A`` is a :class:`scipy.sparse.csc_array` making up the
               depletion matrix
             * ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions
               for a given material in atoms/cm3
             * ``t`` is a float of the time step size in seconds, and
+            * ``substeps`` is an optional integer number of substeps, and
             * ``n1`` is a :class:`numpy.ndarray` of compositions at the
               next time step. Expected to be of the same shape as ``n0``
 
@@ -1322,7 +1333,7 @@ class DepSystemSolver(ABC):
     """
 
     @abstractmethod
-    def __call__(self, A, n0, dt):
+    def __call__(self, A, n0, dt, substeps=1):
         """Solve the linear system of equations for depletion
 
         Parameters
@@ -1335,6 +1346,8 @@ def __call__(self, A, n0, dt):
             material or an atom density
         dt : float
             Time [s] of the specific interval to be solved
+        substeps : int, optional
+            Number of substeps to use when the solver supports substepping.
 
         Returns
         -------
diff --git a/openmc/deplete/cram.py b/openmc/deplete/cram.py
index a5a4a1f1e1b..f9e3ab5d21d 100644
--- a/openmc/deplete/cram.py
+++ b/openmc/deplete/cram.py
@@ -39,10 +39,6 @@ class IPFCramSolver(DepSystemSolver):
         Complex poles. Must have an equal size as ``alpha``.
     alpha0 : float
         Limit of the approximation at infinity
-    substeps : int, optional
-        Number of substeps per depletion interval.
-
-        .. versionadded:: 0.15.4
 
     Attributes
     ----------
@@ -53,12 +49,10 @@ class IPFCramSolver(DepSystemSolver):
         Complex poles :math:`\theta` of the rational approximation
     alpha0 : float
         Limit of the approximation at infinity
-    substeps : int
-        Number of substeps per depletion interval
 
     """
 
-    def __init__(self, alpha, theta, alpha0, substeps: int = 1):
+    def __init__(self, alpha, theta, alpha0):
         check_type("alpha", alpha, np.ndarray, numbers.Complex)
         check_type("theta", theta, np.ndarray, numbers.Complex)
         check_length("theta", theta, alpha.size)
@@ -66,19 +60,8 @@ def __init__(self, alpha, theta, alpha0, substeps: int = 1):
         self.alpha = alpha
         self.theta = theta
         self.alpha0 = alpha0
-        self.substeps = substeps
-
-    @property
-    def substeps(self):
-        return self._substeps
 
-    @substeps.setter
-    def substeps(self, value):
-        check_type("substeps", value, numbers.Integral)
-        check_greater_than("substeps", value, 0)
-        self._substeps = value
-
-    def __call__(self, A, n0, dt):
+    def __call__(self, A, n0, dt, substeps=1):
         """Solve depletion equations using IPF CRAM
 
         Parameters
@@ -91,6 +74,8 @@ def __call__(self, A, n0, dt):
             material or an atom density
         dt : float
             Time [s] of the specific interval to be solved
+        substeps : int, optional
+            Number of substeps per depletion interval.
 
         Returns
         -------
@@ -98,7 +83,10 @@ def __call__(self, A, n0, dt):
             Final compositions after ``dt``
 
         """
-        if self.substeps == 1:
+        check_type("substeps", substeps, numbers.Integral)
+        check_greater_than("substeps", substeps, 0)
+
+        if substeps == 1:
             A = dt * csc_array(A, dtype=np.float64)
             y = n0.copy()
             ident = eye_array(A.shape[0], format='csc')
@@ -107,14 +95,14 @@ def __call__(self, A, n0, dt):
             return y * self.alpha0
 
         # Substep path: pre-compute LU factorizations, reuse across substeps
-        sub_dt = dt / self.substeps
+        sub_dt = dt / substeps
         A_sub = sub_dt * csc_array(A, dtype=np.float64)
         ident = eye_array(A_sub.shape[0], format='csc')
         lu_solvers = [sla.splu(A_sub - theta * ident)
                       for theta in self.theta]
 
         y = n0.copy()
-        for _ in range(self.substeps):
+        for _ in range(substeps):
             for alpha, lu in zip(self.alpha, lu_solvers):
                 y += 2 * np.real(alpha * lu.solve(y))
             y *= self.alpha0
diff --git a/openmc/deplete/pool.py b/openmc/deplete/pool.py
index 58f90894b69..166a2399e0a 100644
--- a/openmc/deplete/pool.py
+++ b/openmc/deplete/pool.py
@@ -42,14 +42,17 @@ def _distribute(items):
         j += chunk_size
 
 def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
-            transfer_rates=None, external_source_rates=None, *matrix_args):
+            transfer_rates=None, external_source_rates=None, substeps=1,
+            *matrix_args):
     """Deplete materials using given reaction rates for a specified time
 
     Parameters
     ----------
     func : callable
         Function to use to get new compositions. Expected to have the signature
-        ``func(A, n0, t) -> n1``
+        ``func(A, n0, t) -> n1`` for default single-step behavior or
+        ``func(A, n0, t, substeps=1) -> n1`` to support non-default
+        substeps.
     chain : openmc.deplete.Chain
         Depletion chain
     n : list of numpy.ndarray
@@ -74,6 +77,8 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
         External source rates for continuous removal/feed.
 
         .. versionadded:: 0.15.3
+    substeps : int, optional
+        Number of substeps to pass to solvers that support substepping.
     matrix_args: Any, optional
         Additional arguments passed to matrix_func
 
@@ -164,7 +169,10 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
 
                 # Concatenate vectors of nuclides in one
                 n_multi = np.concatenate(n)
-                n_result = func(matrix, n_multi, dt)
+                if substeps == 1:
+                    n_result = func(matrix, n_multi, dt)
+                else:
+                    n_result = func(matrix, n_multi, dt, substeps)
 
                 # Split back the nuclide vector result into the original form
                 n_result = np.split(n_result, np.cumsum([len(i) for i in n])[:-1])
@@ -198,7 +206,10 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
                 matrix.resize(matrix.shape[1], matrix.shape[1])
                 n[i] = np.append(n[i], 1.0)
 
-    inputs = zip(matrices, n, repeat(dt))
+    if substeps == 1:
+        inputs = zip(matrices, n, repeat(dt))
+    else:
+        inputs = zip(matrices, n, repeat(dt), repeat(substeps))
 
     if USE_MULTIPROCESSING:
         with Pool(NUM_PROCESSES) as pool:
diff --git a/tests/unit_tests/test_deplete_cram.py b/tests/unit_tests/test_deplete_cram.py
index ec7b8add584..64cff3a8b73 100644
--- a/tests/unit_tests/test_deplete_cram.py
+++ b/tests/unit_tests/test_deplete_cram.py
@@ -1,12 +1,13 @@
-""" Tests for cram.py
+"""Tests for cram.py.
 
 Compares a few Mathematica matrix exponentials to CRAM16/CRAM48.
 Tests substep accuracy against self-converged reference solutions.
 """
 
-from pytest import approx
 import numpy as np
+import pytest
 import scipy.sparse as sp
+from pytest import approx
 from openmc.deplete.cram import (CRAM16, CRAM48, Cram16Solver, Cram48Solver,
                                   IPFCramSolver)
 
@@ -39,14 +40,6 @@ def test_CRAM48():
     assert z == approx(z0)
 
 
-# --- Substep tests ---
-
-def test_substeps_default():
-    """Default substeps=1 on module-level solvers."""
-    assert Cram48Solver.substeps == 1
-    assert Cram16Solver.substeps == 1
-
-
 def test_substeps1_matches_original():
     """substeps=1 must be bitwise identical to original spsolve path."""
     x = np.array([1.0, 1.0])
@@ -54,9 +47,7 @@ def test_substeps1_matches_original():
     dt = 0.1
 
     z_orig = CRAM48(mat, x, dt)
-    solver_1 = IPFCramSolver(Cram48Solver.alpha, Cram48Solver.theta,
-                              Cram48Solver.alpha0, substeps=1)
-    z_sub1 = solver_1(mat, x, dt)
+    z_sub1 = CRAM48(mat, x, dt, substeps=1)
 
     np.testing.assert_array_equal(z_sub1, z_orig)
 
@@ -72,13 +63,21 @@ def test_substeps2_matches_two_half_steps():
     z_two = CRAM48(mat, z_half, dt / 2)
 
     # Single call with substeps=2
-    solver_2 = IPFCramSolver(Cram48Solver.alpha, Cram48Solver.theta,
-                              Cram48Solver.alpha0, substeps=2)
-    z_sub2 = solver_2(mat, x, dt)
+    z_sub2 = CRAM48(mat, x, dt, substeps=2)
 
     assert z_sub2 == approx(z_two, rel=1e-12)
 
 
+@pytest.mark.parametrize("substeps", [0, -1])
+def test_invalid_substeps(substeps):
+    """substeps must be a positive integer at call time."""
+    x = np.array([1.0, 1.0])
+    mat = sp.csr_matrix([[-1.0, 0.0], [-2.0, -3.0]])
+
+    with pytest.raises(ValueError, match="substeps"):
+        CRAM48(mat, x, 0.1, substeps=substeps)
+
+
 def test_substeps_self_convergence():
     """Increasing substeps converges toward reference solution.
 
@@ -90,15 +89,11 @@ def test_substeps_self_convergence():
     x = np.array([1.0, 1.0])
     dt = 50  # lambda*dt = 50 and 150, stresses CRAM16
 
-    ref_solver = IPFCramSolver(Cram16Solver.alpha, Cram16Solver.theta,
-                                Cram16Solver.alpha0, substeps=128)
-    n_ref = ref_solver(mat, x, dt)
+    n_ref = CRAM16(mat, x, dt, substeps=128)
 
     prev_err = np.inf
     for s in [1, 2, 4, 8, 16]:
-        solver = IPFCramSolver(Cram16Solver.alpha, Cram16Solver.theta,
-                                Cram16Solver.alpha0, substeps=s)
-        n_s = solver(mat, x, dt)
+        n_s = CRAM16(mat, x, dt, substeps=s)
         err = np.linalg.norm(n_s - n_ref) / np.linalg.norm(n_ref)
         assert err < prev_err, \
             f"substeps={s} error {err:.2e} not less than previous {prev_err:.2e}"
diff --git a/tests/unit_tests/test_deplete_integrator.py b/tests/unit_tests/test_deplete_integrator.py
index 558cb434bee..b083c205f84 100644
--- a/tests/unit_tests/test_deplete_integrator.py
+++ b/tests/unit_tests/test_deplete_integrator.py
@@ -19,7 +19,7 @@
 from openmc.deplete import (
     ReactionRates, StepResult, Results, OperatorResult, PredictorIntegrator,
     CECMIntegrator, CF4Integrator, CELIIntegrator, EPCRK4Integrator,
-    LEQIIntegrator, SICELIIntegrator, SILEQIIntegrator, cram)
+    LEQIIntegrator, SICELIIntegrator, SILEQIIntegrator, cram, pool)
 
 from tests import dummy_operator
 
@@ -183,18 +183,36 @@ def test_bad_integrator_inputs():
     with pytest.raises(TypeError, match=".*callable.*NoneType"):
         PredictorIntegrator(op, timesteps, power=1, solver=None)
 
-    with pytest.raises(ValueError, match=".*arguments"):
+    with pytest.raises(ValueError, match="three arguments|four arguments"):
         PredictorIntegrator(op, timesteps, power=1, solver=mock_bad_solver_nargs)
 
+    with pytest.raises(ValueError, match="default to 1"):
+        PredictorIntegrator(op, timesteps, power=1,
+                            solver=mock_bad_solver_fourth_required)
+
+    with pytest.raises(ValueError, match="substeps"):
+        PredictorIntegrator(op, timesteps, power=1, substeps=0)
+
+    with pytest.raises(ValueError, match="substeps"):
+        PredictorIntegrator(op, timesteps, power=1, substeps=-1)
+
 
 def mock_good_solver(A, n, t):
-    pass
+    return n.copy()
+
+
+def mock_good_solver_substeps(A, n, t, substeps=1):
+    return n + substeps
 
 
 def mock_bad_solver_nargs(A, n):
     pass
 
 
+def mock_bad_solver_fourth_required(A, n, t, substeps):
+    pass
+
+
 @pytest.mark.parametrize("scheme", dummy_operator.SCHEMES)
 def test_integrator(run_in_tmpdir, scheme):
     """Test the integrators against their expected values"""
@@ -226,6 +244,11 @@ def test_integrator(run_in_tmpdir, scheme):
     integrator = bundle.solver(operator, [0.75], 1, solver="cram16")
     assert integrator.solver is cram.CRAM16
 
+    integrator = bundle.solver(operator, [0.75], 1, solver=cram.Cram48Solver,
+                               substeps=2)
+    assert integrator.solver is cram.Cram48Solver
+    assert integrator.substeps == 2
+
     integrator.solver = mock_good_solver
     assert integrator.solver is mock_good_solver
 
@@ -233,6 +256,48 @@ def test_integrator(run_in_tmpdir, scheme):
     integrator.solver = lfunc
     assert integrator.solver is lfunc
 
+    integrator.solver = mock_good_solver_substeps
+    assert integrator.solver is mock_good_solver_substeps
+
+
+def test_legacy_custom_solver_with_default_substeps(monkeypatch):
+    operator = dummy_operator.DummyOperator()
+    n = operator.initial_condition()
+    rates = operator(n, 1.0).rates
+    integrator = PredictorIntegrator(
+        operator, [0.75], power=1.0, solver=mock_good_solver)
+    monkeypatch.setattr(pool, "USE_MULTIPROCESSING", False)
+
+    _, result = integrator._timed_deplete(n, rates, 0.75)
+
+    np.testing.assert_array_equal(result[0], n[0])
+
+
+def test_legacy_custom_solver_raises_with_substeps(monkeypatch):
+    operator = dummy_operator.DummyOperator()
+    n = operator.initial_condition()
+    rates = operator(n, 1.0).rates
+    integrator = PredictorIntegrator(
+        operator, [0.75], power=1.0, solver=mock_good_solver, substeps=2)
+    monkeypatch.setattr(pool, "USE_MULTIPROCESSING", False)
+
+    with pytest.raises(TypeError):
+        integrator._timed_deplete(n, rates, 0.75)
+
+
+def test_substep_aware_custom_solver_receives_substeps(monkeypatch):
+    operator = dummy_operator.DummyOperator()
+    n = operator.initial_condition()
+    rates = operator(n, 1.0).rates
+    integrator = PredictorIntegrator(
+        operator, [0.75], power=1.0, solver=mock_good_solver_substeps,
+        substeps=3)
+    monkeypatch.setattr(pool, "USE_MULTIPROCESSING", False)
+
+    _, result = integrator._timed_deplete(n, rates, 0.75)
+
+    np.testing.assert_array_equal(result[0], n[0] + 3)
+
 
 @pytest.mark.parametrize("integrator", INTEGRATORS)
 def test_timesteps(integrator):

From c96b984f7752ae2452ad0bd2df19f1d24803b88d Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Wed, 22 Apr 2026 13:57:27 -0500
Subject: [PATCH 5/8] Custom solvers must implement 4-argument form

---
 openmc/deplete/abc.py                       | 24 +++++-------
 openmc/deplete/pool.py                      | 14 ++-----
 tests/unit_tests/test_deplete_integrator.py | 43 ++++++++++++++-------
 3 files changed, 42 insertions(+), 39 deletions(-)

diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py
index 1bfd7b1f16b..2f858d2b7bf 100644
--- a/openmc/deplete/abc.py
+++ b/openmc/deplete/abc.py
@@ -606,12 +606,9 @@ class Integrator(ABC):
         Function that will solve the Bateman equations
         :math:`\frac{\partial}{\partial t}\vec{n} = A_i\vec{n}_i` with a step
         size :math:`t_i`. Can be configured using the ``solver`` argument.
-        User-supplied functions may have either of the following signatures:
-        ``solver(A, n0, t) -> n1`` for the default single-step behavior or
-        ``solver(A, n0, t, substeps=1) -> n1`` to support non-default
-        substeps. When :attr:`substeps` is greater than one, OpenMC will call
-        the solver with the fourth argument and unsupported custom solvers will
-        raise at runtime. In either case,
+        User-supplied functions are expected to have the following signature:
+        ``solver(A, n0, t, substeps=1) -> n1``. Solvers that do not support
+        multiple substeps may raise when ``substeps > 1``.
 
             * ``A`` is a :class:`scipy.sparse.csc_array` making up the
               depletion matrix
@@ -737,10 +734,10 @@ def solver(self, func):
         params = list(sig.parameters.values())
 
         # Inspect arguments
-        if len(params) not in {3, 4}:
+        if len(params) != 4:
             raise ValueError(
-                "Function {} must support either three arguments "
-                "(A, n0, t) or four arguments (A, n0, t, substeps=1): {!s}"
+                "Function {} must support four arguments "
+                "(A, n0, t, substeps=1): {!s}"
                 .format(func, sig))
 
         for ix, param in enumerate(params):
@@ -1164,12 +1161,9 @@ class SIIntegrator(Integrator):
         Function that will solve the Bateman equations
         :math:`\frac{\partial}{\partial t}\vec{n} = A_i\vec{n}_i` with a step
         size :math:`t_i`. Can be configured using the ``solver`` argument.
-        User-supplied functions may have either of the following signatures:
-        ``solver(A, n0, t) -> n1`` for the default single-step behavior or
-        ``solver(A, n0, t, substeps=1) -> n1`` to support non-default
-        substeps. When :attr:`substeps` is greater than one, OpenMC will call
-        the solver with the fourth argument and unsupported custom solvers will
-        raise at runtime. In either case,
+        User-supplied functions are expected to have the following signature:
+        ``solver(A, n0, t, substeps=1) -> n1``. Solvers that do not support
+        multiple substeps may raise when ``substeps > 1``.
 
             * ``A`` is a :class:`scipy.sparse.csc_array` making up the
               depletion matrix
diff --git a/openmc/deplete/pool.py b/openmc/deplete/pool.py
index 166a2399e0a..19ad0ada503 100644
--- a/openmc/deplete/pool.py
+++ b/openmc/deplete/pool.py
@@ -50,9 +50,7 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
     ----------
     func : callable
         Function to use to get new compositions. Expected to have the signature
-        ``func(A, n0, t) -> n1`` for default single-step behavior or
-        ``func(A, n0, t, substeps=1) -> n1`` to support non-default
-        substeps.
+        ``func(A, n0, t, substeps=1) -> n1``.
     chain : openmc.deplete.Chain
         Depletion chain
     n : list of numpy.ndarray
@@ -169,10 +167,7 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
 
                 # Concatenate vectors of nuclides in one
                 n_multi = np.concatenate(n)
-                if substeps == 1:
-                    n_result = func(matrix, n_multi, dt)
-                else:
-                    n_result = func(matrix, n_multi, dt, substeps)
+                n_result = func(matrix, n_multi, dt, substeps)
 
                 # Split back the nuclide vector result into the original form
                 n_result = np.split(n_result, np.cumsum([len(i) for i in n])[:-1])
@@ -206,10 +201,7 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
                 matrix.resize(matrix.shape[1], matrix.shape[1])
                 n[i] = np.append(n[i], 1.0)
 
-    if substeps == 1:
-        inputs = zip(matrices, n, repeat(dt))
-    else:
-        inputs = zip(matrices, n, repeat(dt), repeat(substeps))
+    inputs = zip(matrices, n, repeat(dt), repeat(substeps))
 
     if USE_MULTIPROCESSING:
         with Pool(NUM_PROCESSES) as pool:
diff --git a/tests/unit_tests/test_deplete_integrator.py b/tests/unit_tests/test_deplete_integrator.py
index b083c205f84..ec886e70795 100644
--- a/tests/unit_tests/test_deplete_integrator.py
+++ b/tests/unit_tests/test_deplete_integrator.py
@@ -183,7 +183,7 @@ def test_bad_integrator_inputs():
     with pytest.raises(TypeError, match=".*callable.*NoneType"):
         PredictorIntegrator(op, timesteps, power=1, solver=None)
 
-    with pytest.raises(ValueError, match="three arguments|four arguments"):
+    with pytest.raises(ValueError, match="four arguments"):
         PredictorIntegrator(op, timesteps, power=1, solver=mock_bad_solver_nargs)
 
     with pytest.raises(ValueError, match="default to 1"):
@@ -197,7 +197,7 @@ def test_bad_integrator_inputs():
         PredictorIntegrator(op, timesteps, power=1, substeps=-1)
 
 
-def mock_good_solver(A, n, t):
+def mock_good_solver(A, n, t, substeps=1):
     return n.copy()
 
 
@@ -205,6 +205,12 @@ def mock_good_solver_substeps(A, n, t, substeps=1):
     return n + substeps
 
 
+def mock_unsupported_substeps_solver(A, n, t, substeps=1):
+    if substeps > 1:
+        raise NotImplementedError("substeps > 1 not supported")
+    return n.copy()
+
+
 def mock_bad_solver_nargs(A, n):
     pass
 
@@ -252,7 +258,7 @@ def test_integrator(run_in_tmpdir, scheme):
     integrator.solver = mock_good_solver
     assert integrator.solver is mock_good_solver
 
-    lfunc = lambda A, n, t: mock_good_solver(A, n, t)
+    lfunc = lambda A, n, t, substeps=1: mock_good_solver(A, n, t, substeps)
     integrator.solver = lfunc
     assert integrator.solver is lfunc
 
@@ -260,7 +266,7 @@ def test_integrator(run_in_tmpdir, scheme):
     assert integrator.solver is mock_good_solver_substeps
 
 
-def test_legacy_custom_solver_with_default_substeps(monkeypatch):
+def test_custom_solver_with_default_substeps(monkeypatch):
     operator = dummy_operator.DummyOperator()
     n = operator.initial_condition()
     rates = operator(n, 1.0).rates
@@ -273,30 +279,41 @@ def test_legacy_custom_solver_with_default_substeps(monkeypatch):
     np.testing.assert_array_equal(result[0], n[0])
 
 
-def test_legacy_custom_solver_raises_with_substeps(monkeypatch):
+def test_substep_aware_custom_solver_receives_substeps(monkeypatch):
     operator = dummy_operator.DummyOperator()
     n = operator.initial_condition()
     rates = operator(n, 1.0).rates
     integrator = PredictorIntegrator(
-        operator, [0.75], power=1.0, solver=mock_good_solver, substeps=2)
+        operator, [0.75], power=1.0, solver=mock_good_solver_substeps,
+        substeps=3)
     monkeypatch.setattr(pool, "USE_MULTIPROCESSING", False)
 
-    with pytest.raises(TypeError):
-        integrator._timed_deplete(n, rates, 0.75)
+    _, result = integrator._timed_deplete(n, rates, 0.75)
 
+    np.testing.assert_array_equal(result[0], n[0] + 3)
 
-def test_substep_aware_custom_solver_receives_substeps(monkeypatch):
+
+def test_custom_solver_propagates_substeps_error(monkeypatch):
     operator = dummy_operator.DummyOperator()
     n = operator.initial_condition()
     rates = operator(n, 1.0).rates
     integrator = PredictorIntegrator(
-        operator, [0.75], power=1.0, solver=mock_good_solver_substeps,
-        substeps=3)
+        operator, [0.75], power=1.0,
+        solver=mock_unsupported_substeps_solver, substeps=2)
     monkeypatch.setattr(pool, "USE_MULTIPROCESSING", False)
 
-    _, result = integrator._timed_deplete(n, rates, 0.75)
+    with pytest.raises(NotImplementedError, match="not supported"):
+        integrator._timed_deplete(n, rates, 0.75)
 
-    np.testing.assert_array_equal(result[0], n[0] + 3)
+
+def test_custom_solver_requires_four_args():
+    op = MagicMock()
+    op.prev_res = None
+    op.chain = None
+    op.heavy_metal = 1.0
+
+    with pytest.raises(ValueError, match="four arguments"):
+        PredictorIntegrator(op, [1], power=1, solver=mock_bad_solver_nargs)
 
 
 @pytest.mark.parametrize("integrator", INTEGRATORS)

From 92ba174c228ab534fe1fc923e312304a9664a0ca Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Wed, 22 Apr 2026 14:00:53 -0500
Subject: [PATCH 6/8] Update docstrings

---
 openmc/deplete/abc.py | 130 +++++++++++++++++++++---------------------
 1 file changed, 65 insertions(+), 65 deletions(-)

diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py
index 2f858d2b7bf..6a7a01a51be 100644
--- a/openmc/deplete/abc.py
+++ b/openmc/deplete/abc.py
@@ -540,17 +540,15 @@ class Integrator(ABC):
         iterable of float. Alternatively, units can be specified for each step
         by passing an iterable of (value, unit) tuples.
     power : float or iterable of float, optional
-        Power of the reactor in [W]. A single value indicates that
-        the power is constant over all timesteps. An iterable
-        indicates potentially different power levels for each timestep.
-        For a 2D problem, the power can be given in [W/cm] as long
-        as the "volume" assigned to a depletion material is actually
-        an area in [cm^2]. Either ``power``, ``power_density``, or
+        Power of the reactor in [W]. A single value indicates that the power is
+        constant over all timesteps. An iterable indicates potentially different
+        power levels for each timestep. For a 2D problem, the power can be given
+        in [W/cm] as long as the "volume" assigned to a depletion material is
+        actually an area in [cm^2]. Either ``power``, ``power_density``, or
         ``source_rates`` must be specified.
     power_density : float or iterable of float, optional
-        Power density of the reactor in [W/gHM]. It is multiplied by
-        initial heavy metal inventory to get total power if ``power``
-        is not specified.
+        Power density of the reactor in [W/gHM]. It is multiplied by initial
+        heavy metal inventory to get total power if ``power`` is not specified.
     source_rates : float or iterable of float, optional
         Source rate in [neutron/sec] or neutron flux in [neutron/s-cm^2] for
         each interval in :attr:`timesteps`
@@ -562,8 +560,8 @@ class Integrator(ABC):
         and 'MWd/kg' indicates that the values are given in burnup (MW-d of
         energy deposited per kilogram of initial heavy metal).
     solver : str or callable, optional
-        If a string, must be the name of the solver responsible for
-        solving the Bateman equations.  Current options are:
+        If a string, must be the name of the solver responsible for solving the
+        Bateman equations.  Current options are:
 
             * ``cram16`` - 16th order IPF CRAM
             * ``cram48`` - 48th order IPF CRAM [default]
@@ -573,21 +571,21 @@ class Integrator(ABC):
 
         .. versionadded:: 0.12
     substeps : int, optional
-        Number of substeps per depletion interval. When greater than 1,
-        each interval is subdivided into `substeps` identical sub-intervals
-        and LU factorizations may be reused across them, improving accuracy
-        for nuclides with large decay-constant × timestep products.
+        Number of substeps per depletion interval. When greater than 1, each
+        interval is subdivided into `substeps` identical sub-intervals and LU
+        factorizations may be reused across them, improving accuracy for
+        nuclides with large decay-constant × timestep products.
 
         .. versionadded:: 0.15.4
     continue_timesteps : bool, optional
         Whether or not to treat the current solve as a continuation of a
         previous simulation. Defaults to `False`. When `False`, the depletion
         steps provided are appended to any previous steps. If `True`, the
-        timesteps provided to the `Integrator` must exacly match any that
-        exist in the `prev_results` passed to the `Operator`. The `power`,
-        `power_density`, or `source_rates` must match as well. The
-        method of specifying `power`, `power_density`, or
-        `source_rates` should be the same as the initial run.
+        timesteps provided to the `Integrator` must exacly match any that exist
+        in the `prev_results` passed to the `Operator`. The `power`,
+        `power_density`, or `source_rates` must match as well. The method of
+        specifying `power`, `power_density`, or `source_rates` should be the
+        same as the initial run.
 
         .. versionadded:: 0.15.1
 
@@ -607,17 +605,19 @@ class Integrator(ABC):
         :math:`\frac{\partial}{\partial t}\vec{n} = A_i\vec{n}_i` with a step
         size :math:`t_i`. Can be configured using the ``solver`` argument.
         User-supplied functions are expected to have the following signature:
-        ``solver(A, n0, t, substeps=1) -> n1``. Solvers that do not support
-        multiple substeps may raise when ``substeps > 1``.
-
-            * ``A`` is a :class:`scipy.sparse.csc_array` making up the
-              depletion matrix
-            * ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions
-              for a given material in atoms/cm3
-            * ``t`` is a float of the time step size in seconds, and
+        ``solver(A, n0, t, substeps=1) -> n1``, where
+
+            * ``A`` is a :class:`scipy.sparse.csc_array` making up the depletion
+              matrix
+            * ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions for
+              a given material in atoms/cm3
+            * ``t`` is a float of the time step size in seconds
             * ``substeps`` is an optional integer number of substeps, and
-            * ``n1`` is a :class:`numpy.ndarray` of compositions at the
-              next time step. Expected to be of the same shape as ``n0``
+            * ``n1`` is a :class:`numpy.ndarray` of compositions at the next
+              time step. Expected to be of the same shape as ``n0``
+
+        Solvers that do not support multiple substeps should raise an exception
+        when ``substeps > 1``.
 
     transfer_rates : openmc.deplete.TransferRates
         Transfer rates for the depletion system used to model continuous
@@ -1090,17 +1090,15 @@ class SIIntegrator(Integrator):
         iterable of float. Alternatively, units can be specified for each step
         by passing an iterable of (value, unit) tuples.
     power : float or iterable of float, optional
-        Power of the reactor in [W]. A single value indicates that
-        the power is constant over all timesteps. An iterable
-        indicates potentially different power levels for each timestep.
-        For a 2D problem, the power can be given in [W/cm] as long
-        as the "volume" assigned to a depletion material is actually
-        an area in [cm^2]. Either ``power``, ``power_density``, or
+        Power of the reactor in [W]. A single value indicates that the power is
+        constant over all timesteps. An iterable indicates potentially different
+        power levels for each timestep. For a 2D problem, the power can be given
+        in [W/cm] as long as the "volume" assigned to a depletion material is
+        actually an area in [cm^2]. Either ``power``, ``power_density``, or
         ``source_rates`` must be specified.
     power_density : float or iterable of float, optional
-        Power density of the reactor in [W/gHM]. It is multiplied by
-        initial heavy metal inventory to get total power if ``power``
-        is not specified.
+        Power density of the reactor in [W/gHM]. It is multiplied by initial
+        heavy metal inventory to get total power if ``power`` is not specified.
     source_rates : float or iterable of float, optional
         Source rate in [neutron/sec] or neutron flux in [neutron/s-cm^2] for
         each interval in :attr:`timesteps`
@@ -1112,11 +1110,11 @@ class SIIntegrator(Integrator):
         that the values are given in burnup (MW-d of energy deposited per
         kilogram of initial heavy metal).
     n_steps : int, optional
-        Number of stochastic iterations per depletion interval.
-        Must be greater than zero. Default : 10
+        Number of stochastic iterations per depletion interval. Must be greater
+        than zero. Default : 10
     solver : str or callable, optional
-        If a string, must be the name of the solver responsible for
-        solving the Bateman equations.  Current options are:
+        If a string, must be the name of the solver responsible for solving the
+        Bateman equations.  Current options are:
 
             * ``cram16`` - 16th order IPF CRAM
             * ``cram48`` - 48th order IPF CRAM [default]
@@ -1126,22 +1124,22 @@ class SIIntegrator(Integrator):
 
         .. versionadded:: 0.12
     substeps : int, optional
-        Number of substeps per depletion interval. When greater than 1,
-        each interval is subdivided into `substeps` identical sub-intervals
-        and LU factorizations may be reused across them, improving accuracy
-        for nuclides with large decay-constant × timestep products.
+        Number of substeps per depletion interval. When greater than 1, each
+        interval is subdivided into `substeps` identical sub-intervals and LU
+        factorizations may be reused across them, improving accuracy for
+        nuclides with large decay-constant × timestep products.
 
         .. versionadded:: 0.15.4
     continue_timesteps : bool, optional
         Whether or not to treat the current solve as a continuation of a
-        previous simulation. Defaults to `False`. If `False`, all time
-        steps and source rates will be run in an append fashion and will run
-        after whatever time steps exist, if any. If `True`, the timesteps
-        provided to the `Integrator` must match exactly those that exist
-        in the `prev_results` passed to the `Opereator`. The `power`,
-        `power_density`, or `source_rates` must match as well. The
-        method of specifying `power`, `power_density`, or
-        `source_rates` should be the same as the initial run.
+        previous simulation. Defaults to `False`. If `False`, all time steps and
+        source rates will be run in an append fashion and will run after
+        whatever time steps exist, if any. If `True`, the timesteps provided to
+        the `Integrator` must match exactly those that exist in the
+        `prev_results` passed to the `Opereator`. The `power`, `power_density`,
+        or `source_rates` must match as well. The method of specifying `power`,
+        `power_density`, or `source_rates` should be the same as the initial
+        run.
 
         .. versionadded:: 0.15.1
 
@@ -1162,17 +1160,19 @@ class SIIntegrator(Integrator):
         :math:`\frac{\partial}{\partial t}\vec{n} = A_i\vec{n}_i` with a step
         size :math:`t_i`. Can be configured using the ``solver`` argument.
         User-supplied functions are expected to have the following signature:
-        ``solver(A, n0, t, substeps=1) -> n1``. Solvers that do not support
-        multiple substeps may raise when ``substeps > 1``.
-
-            * ``A`` is a :class:`scipy.sparse.csc_array` making up the
-              depletion matrix
-            * ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions
-              for a given material in atoms/cm3
-            * ``t`` is a float of the time step size in seconds, and
+        ``solver(A, n0, t, substeps=1) -> n1``, where
+
+            * ``A`` is a :class:`scipy.sparse.csc_array` making up the depletion
+              matrix
+            * ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions for
+              a given material in atoms/cm3
+            * ``t`` is a float of the time step size in seconds
             * ``substeps`` is an optional integer number of substeps, and
-            * ``n1`` is a :class:`numpy.ndarray` of compositions at the
-              next time step. Expected to be of the same shape as ``n0``
+            * ``n1`` is a :class:`numpy.ndarray` of compositions at the next
+              time step. Expected to be of the same shape as ``n0``
+
+        Solvers that do not support multiple substeps should raise an exception
+        when ``substeps > 1``.
 
         .. versionadded:: 0.12
 

From 1d4f554d9ab5eaaaa19cfaa3eb425ee5d0520295 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Wed, 22 Apr 2026 14:06:31 -0500
Subject: [PATCH 7/8] Restore setting of _solver

---
 openmc/deplete/abc.py | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py
index 6a7a01a51be..8ceb5f72874 100644
--- a/openmc/deplete/abc.py
+++ b/openmc/deplete/abc.py
@@ -703,13 +703,14 @@ def __init__(
         if isinstance(solver, str):
             # Delay importing of cram module, which requires this file
             if solver == "cram48":
-                from .cram import CRAM48 as _default
+                from .cram import CRAM48
+                self._solver = CRAM48
             elif solver == "cram16":
-                from .cram import CRAM16 as _default
+                from .cram import CRAM16
+                self._solver = CRAM16
             else:
                 raise ValueError(
                     f"Solver {solver} not understood. Expected 'cram48' or 'cram16'")
-            self._solver = _default
         else:
             self.solver = solver
 

From a69b78b19e1cdd0836b6066e31eda2cd7d04ba44 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Wed, 22 Apr 2026 14:25:15 -0500
Subject: [PATCH 8/8] Reduce redundant code in IPFCramSolver.__call__

---
 openmc/deplete/cram.py | 28 ++++++++++++----------------
 1 file changed, 12 insertions(+), 16 deletions(-)

diff --git a/openmc/deplete/cram.py b/openmc/deplete/cram.py
index f9e3ab5d21d..3594ffbe841 100644
--- a/openmc/deplete/cram.py
+++ b/openmc/deplete/cram.py
@@ -3,10 +3,11 @@
 Implements two different forms of CRAM for use in openmc.deplete.
 """
 
+from functools import partial
 import numbers
 
 import numpy as np
-import scipy.sparse.linalg as sla
+from scipy.sparse.linalg import spsolve, splu
 
 from openmc.checkvalue import check_type, check_length, check_greater_than
 from .abc import DepSystemSolver
@@ -86,25 +87,20 @@ def __call__(self, A, n0, dt, substeps=1):
         check_type("substeps", substeps, numbers.Integral)
         check_greater_than("substeps", substeps, 0)
 
+        step_dt = dt if substeps == 1 else dt / substeps
+        A = step_dt * csc_array(A, dtype=np.float64)
+        ident = eye_array(A.shape[0], format='csc')
+
         if substeps == 1:
-            A = dt * csc_array(A, dtype=np.float64)
-            y = n0.copy()
-            ident = eye_array(A.shape[0], format='csc')
-            for alpha, theta in zip(self.alpha, self.theta):
-                y += 2*np.real(alpha*sla.spsolve(A - theta*ident, y))
-            return y * self.alpha0
-
-        # Substep path: pre-compute LU factorizations, reuse across substeps
-        sub_dt = dt / substeps
-        A_sub = sub_dt * csc_array(A, dtype=np.float64)
-        ident = eye_array(A_sub.shape[0], format='csc')
-        lu_solvers = [sla.splu(A_sub - theta * ident)
-                      for theta in self.theta]
+            solvers = [partial(spsolve, A - theta * ident) for theta in self.theta]
+        else:
+            # Pre-compute LU factorizations and reuse them across substeps.
+            solvers = [splu(A - theta * ident).solve for theta in self.theta]
 
         y = n0.copy()
         for _ in range(substeps):
-            for alpha, lu in zip(self.alpha, lu_solvers):
-                y += 2 * np.real(alpha * lu.solve(y))
+            for alpha, solve in zip(self.alpha, solvers):
+                y += 2 * np.real(alpha * solve(y))
             y *= self.alpha0
         return y