CRAM substeps

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):

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

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