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 `_," 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 + `_," 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