discrete1.cram module

Chebyshev Rational Approximation Method (CRAM) matrix exponential.

This module evaluates n(t + dt) = exp(A * dt) @ n0 for the depletion (Bateman) system dn/dt = A n using the Chebyshev Rational Approximation Method. CRAM is well suited to burnup matrices because their eigenvalues cluster near the negative real axis, where the Chebyshev rational approximation of exp(x) is near-optimal. This makes it accurate and stable for the stiff matrices that arise from mixing short- and long-lived nuclides, without pruning the chain.

The approximation is evaluated through its incomplete-partial-fraction (IPF) form, applying the half-set of conjugate poles theta_k sequentially:

y <- n0
for each (alpha_k, theta_k):
    y <- y + 2 * Re[ alpha_k * (A dt - theta_k I)^-1 y ]
y <- alpha0 * y

where alpha_k are the residues and alpha0 is a final multiplicative factor. Each pole requires a single complex sparse linear solve. This is the form used by OpenMC/Serpent; note the running vector y (not n0) is fed into each successive solve.

Two orders are provided, selectable via the order argument:

  • order=16: CRAM-16 (Pusa, 2010).

  • order=48: CRAM-48 incomplete partial fractions (Pusa, 2016).

The order-48 set is more accurate and robust for large time steps; order-16 is cheaper. Coefficient values are transcribed from the canonical Pusa references as used by OpenMC/Serpent.

discrete1.cram.cram_expm(matrix, n0, dt, order=48, substeps=1)[source]

Evaluate exp(matrix * dt) @ n0 with CRAM.

Parameters:

  • matrix (scipy.sparse matrix or numpy.ndarray, shape (M, M)) – Burnup (transmutation) matrix A with units of 1/s. Sparse input is recommended; dense input is converted internally.

  • n0 (numpy.ndarray, shape (M,)) – Initial nuclide number densities.

  • dt (float) – Time step in seconds.

  • order (int, optional) – CRAM order, either 16 or 48 (default 48).

  • substeps (int, optional) – Number of equal substeps over dt (default 1). Substeps reuse sparse LU factorizations and reduce error on stiff chains.

Returns:

Number densities after dt, i.e. exp(A dt) n0. The result is real; the tiny imaginary residue from the complex solves is discarded.

Return type:

numpy.ndarray, shape (M,)

Notes

Uses the incomplete-partial-fraction (IPF) form: the half-set of conjugate poles is applied sequentially to the running vector, each pole solving the complex sparse system (A dt - theta_k I) y = y via a sparse LU factorization, with a final multiply by alpha0.