From 07d3207d88ca5c0c8752ae81f0fb8794df298188 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Thu, 19 Sep 2019 08:50:42 -0500
Subject: [PATCH 1/3] Provide class-based CRAM solvers

CRAM16 and CRAM48 are now aliases for __call__ methods
on two new classes: openmc.deplete.cram.Cram16Solver
and Cram48Solver. These are two concrete subclasses of
openmc.deplete.abc.IPFCramSolver and
openmc.deplete.abc.DepSystemSolver abstract classes.

The primary benefit of the class based approach is that the
16th and 48th order have nearly identical implementations. The only
differences are the alpha, theta, and alpha_0 coefficients used.
This allows both the Cram16Solver and Cram48Solver to have unique
sets of coefficient vectors, while relying on the IPFSolver base
class implementation.

The implementation of IPF CRAM has been cleaned up. The NxN identity
matrix is not re-created each iteration. Given the alpha and
theta attributes on the IPFCramSolver instances, one can
iterate through the orders using a zip command. By forcing concrete
classes to declared alpha and theta vectors, the need to re-declare
the vectors at each entrance into CRAM48 is also removed.

Using these classes, a 10% speedup is observed depleting up to 500
materials on a single MPI process. The function-based approach
took 62s, while the new classes required 52s.
---
 openmc/deplete/abc.py  | 117 +++++++++++++++++-
 openmc/deplete/cram.py | 271 +++++++++++++++++++++++------------------
 2 files changed, 267 insertions(+), 121 deletions(-)

diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py
index 565b9e2dc..b42aa3e74 100644
--- a/openmc/deplete/abc.py
+++ b/openmc/deplete/abc.py
@@ -14,7 +14,9 @@ from copy import deepcopy
 from warnings import warn
 from numbers import Real, Integral
 
-from numpy import nonzero, empty, asarray
+from numpy import nonzero, empty, asarray, float64, real
+import scipy.sparse as sp
+import scipy.sparse.linalg as sla
 from uncertainties import ufloat
 
 from openmc.data import DataLibrary, JOULE_PER_EV
@@ -855,3 +857,116 @@ class SIIntegrator(Integrator):
             Results.save(self.operator, [conc], [res_list[-1]], [t, t],
                          p, self._i_res + len(self), proc_time)
             self.operator.write_bos_data(self._i_res + len(self))
+
+
+class DepSystemSolver(ABC):
+    r"""Abstract class for solving depletion equations
+
+    Responsible for solving
+
+    .. math::
+
+        \frac{\partial \vec{N}}{\partial t} = \bar{A}\vec{N}(t),
+
+    for :math:`0< t\leq t +\Delta t`, given :math:`\vec{N}(0) = \vec{N}_0`
+
+    """
+
+    @abstractmethod
+    def __call__(self, A, n0, dt):
+        """Solve the linear system of equations for depletion
+
+        Parameters
+        ----------
+        A : scipy.sparse.csr_matrix
+            Sparse transmutation matrix ``A[j, i]`` desribing rates at
+            which isotope ``i`` transmutes to isotope ``j``
+        n0 : numpy.ndarray
+            Initial compositions, typically given in number of atoms in some
+            material or an atom density
+        dt : float
+            Time [s] of the specific interval to be solved
+
+        Returns
+        -------
+        numpy.ndarray
+            Final compositions after ``dt``. Should be of identical shape
+            to ``n0``.
+
+        """
+
+
+class IPFCramSolver(DepSystemSolver):
+    r"""Abstract class that implements the IPF form of CRAM
+
+    Provides a :meth:`__call__` that utilizes an incomplete
+    partial factorization (IPF) for the Chebyshev Rational Approximation
+    Method (CRAM) [Pusa16]_
+
+    Concrete subclasses must provide two complex vectors :attr:`alpha`
+    and :attr:`theta` that make up the coefficients of the decompostion.
+    Vectors are expected to be of equal length ``N``.
+    Subclases are also expected to provide a coefficient :attr:`alpha0`
+    used in the final scaling step.
+
+    Attributes
+    ----------
+    alpha : numpy.ndarray
+        Complex residues of poles :attr:`theta` in the incomplete partial
+        factorization. Denoted as :math:`\tilde{\alpha}`
+    theta : numpy.ndarray
+        Complex poles :math:`\theta` of the rational approximation
+    alpha0 : float
+        Limit of the approximation at infinity
+
+    References
+    ----------
+
+    .. [Pusa16] M. Pusa, "Higher-Order Chebyshev Rational Approximation
+       Method and Application to Burnup Equations," Nuclear Science And
+       Engineering, 182:3,297-318
+       `DOI: 10.13182/NSE15-26 <https://doi.org/10.13182/NSE15-26>`_
+
+    """
+
+    @property
+    @abstractmethod
+    def alpha(self):
+        pass
+
+    @property
+    @abstractmethod
+    def theta(self):
+        pass
+
+    @property
+    @abstractmethod
+    def alpha0(self):
+        pass
+
+    def __call__(self, A, n0, dt):
+        """Solve depletion equations using IPF CRAM
+
+        Parameters
+        ----------
+        A : scipy.sparse.csr_matrix
+            Sparse transmutation matrix ``A[j, i]`` desribing rates at
+            which isotope ``i`` transmutes to isotope ``j``
+        n0 : numpy.ndarray
+            Initial compositions, typically given in number of atoms in some
+            material or an atom density
+        dt : float
+            Time [s] of the specific interval to be solved
+
+        Returns
+        -------
+        numpy.ndarray
+            Final compositions after ``dt``
+
+        """
+        A = sp.csr_matrix(A * dt, dtype=float64)
+        y = asarray(n0, dtype=float64)
+        ident = sp.eye(A.shape[0])
+        for alpha, theta in zip(self.alpha, self.theta):
+            y += 2*real(alpha*sla.spsolve(A - theta*ident, y))
+        return y * self.alpha0
diff --git a/openmc/deplete/cram.py b/openmc/deplete/cram.py
index 31049bb72..4fdba87bf 100644
--- a/openmc/deplete/cram.py
+++ b/openmc/deplete/cram.py
@@ -8,10 +8,12 @@ from multiprocessing import Pool
 import time
 
 import numpy as np
-import scipy.sparse as sp
-import scipy.sparse.linalg as sla
 
-from . import comm
+from .abc import IPFCramSolver
+
+__all__ = [
+    "deplete", "timed_deplete", "CRAM16", "CRAM48",
+    "Cram16Solver", "Cram48Solver"]
 
 
 def deplete(chain, x, rates, dt, matrix_func=None):
@@ -81,147 +83,176 @@ def timed_deplete(*args, **kwargs):
     return time.time() - start, results
 
 
-def CRAM16(A, n0, dt):
-    """Chebyshev Rational Approximation Method, order 16
+class Cram16Solver(IPFCramSolver):
+    r"""Solver implementing the 16th order IPF CRAM
 
-    Algorithm is the 16th order Chebyshev Rational Approximation Method,
-    implemented in the more stable `incomplete partial fraction (IPF)
-    <https://doi.org/10.13182/NSE15-26>`_ form.
+    Coefficients to :attr:`alpha`, :attr:`theta`, and :attr:`alpha0`
+    are from Table A.IV in [Pusa16]_.
 
-    Parameters
+    Attributes
     ----------
-    A : scipy.linalg.csr_matrix
-        Matrix to take exponent of.
-    n0 : numpy.array
-        Vector to operate a matrix exponent on.
-    dt : float
-        Time to integrate to.
-
-    Returns
-    -------
-    numpy.array
-        Results of the matrix exponent.
+    alpha : numpy.ndarray
+        Complex residues of poles :attr:`theta` in the incomplete partial
+        factorization. Denoted as :math:`\tilde{\alpha}`
+        in Algorithm 1 of [Pusa16]_
+    theta : numpy.ndarray
+        Complex poles :math:`\theta` of the rational approximation
+    alpha0 : float
+        Limit of the approximation at infinity
 
     """
+    alpha = np.array([
+        +5.464930576870210e+3 - 3.797983575308356e+4j,
+        +9.045112476907548e+1 - 1.115537522430261e+3j,
+        +2.344818070467641e+2 - 4.228020157070496e+2j,
+        +9.453304067358312e+1 - 2.951294291446048e+2j,
+        +7.283792954673409e+2 - 1.205646080220011e+5j,
+        +3.648229059594851e+1 - 1.155509621409682e+2j,
+        +2.547321630156819e+1 - 2.639500283021502e+1j,
+        +2.394538338734709e+1 - 5.650522971778156e+0j],
+        dtype=np.complex128)
 
-    alpha = np.array([+2.124853710495224e-16,
-                      +5.464930576870210e+3 - 3.797983575308356e+4j,
-                      +9.045112476907548e+1 - 1.115537522430261e+3j,
-                      +2.344818070467641e+2 - 4.228020157070496e+2j,
-                      +9.453304067358312e+1 - 2.951294291446048e+2j,
-                      +7.283792954673409e+2 - 1.205646080220011e+5j,
-                      +3.648229059594851e+1 - 1.155509621409682e+2j,
-                      +2.547321630156819e+1 - 2.639500283021502e+1j,
-                      +2.394538338734709e+1 - 5.650522971778156e+0j],
-                     dtype=np.complex128)
-    theta = np.array([+0.0,
-                      +3.509103608414918 + 8.436198985884374j,
-                      +5.948152268951177 + 3.587457362018322j,
-                      -5.264971343442647 + 16.22022147316793j,
-                      +1.419375897185666 + 10.92536348449672j,
-                      +6.416177699099435 + 1.194122393370139j,
-                      +4.993174737717997 + 5.996881713603942j,
-                      -1.413928462488886 + 13.49772569889275j,
-                      -10.84391707869699 + 19.27744616718165j],
-                     dtype=np.complex128)
-
-    n = A.shape[0]
+    theta = np.array([
+        +3.509103608414918 + 8.436198985884374j,
+        +5.948152268951177 + 3.587457362018322j,
+        -5.264971343442647 + 16.22022147316793j,
+        +1.419375897185666 + 10.92536348449672j,
+        +6.416177699099435 + 1.194122393370139j,
+        +4.993174737717997 + 5.996881713603942j,
+        -1.413928462488886 + 13.49772569889275j,
+        -10.84391707869699 + 19.27744616718165j],
+        dtype=np.complex128)
 
     alpha0 = 2.124853710495224e-16
 
-    k = 8
+    def __call__(self, A, n0, dt):
+        """Solve using 16th order IPF CRAM [Pusa16]_
 
-    y = np.array(n0, dtype=np.float64)
-    for l in range(1, k+1):
-        y = 2.0*np.real(alpha[l]*sla.spsolve(A*dt - theta[l]*sp.eye(n), y)) + y
+        Parameters
+        ----------
+        A : scipy.sparse.csr_matrix
+            Sparse transmutation matrix ``A[j, i]`` desribing rates at
+            which isotope ``i`` transmutes to isotope ``j``
+        n0 : numpy.ndarray
+            Initial compositions, typically given in number of atoms in some
+            material or an atom density
+        dt : float
+            Time [s] of the specific interval to be solved
 
-    y *= alpha0
-    return y
+        Returns
+        -------
+        numpy.ndarray
+            Final compositions after ``dt``
+
+        """
+        return super().__call__(A, n0, dt)
 
 
-def CRAM48(A, n0, dt):
-    """Chebyshev Rational Approximation Method, order 48
+class Cram48Solver(IPFCramSolver):
+    r"""Solver implementing the 48th order IPF CRAM
 
-    Algorithm is the 48th order Chebyshev Rational Approximation Method,
-    implemented in the more stable `incomplete partial fraction (IPF)
-    <https://doi.org/10.13182/NSE15-26>`_ form.
+    Coefficients to :attr:`alpha`, :attr:`theta`, and :attr:`alpha0`
+    are from Table A.XII in [Pusa16]_.
 
-    Parameters
+    Attributes
     ----------
-    A : scipy.linalg.csr_matrix
-        Matrix to take exponent of.
-    n0 : numpy.array
-        Vector to operate a matrix exponent on.
-    dt : float
-        Time to integrate to.
-
-    Returns
-    -------
-    numpy.array
-        Results of the matrix exponent.
+    alpha : numpy.ndarray
+        Complex residues of poles :attr:`theta` in the incomplete partial
+        factorization. Denoted as :math:`\tilde{\alpha}`
+        in Algorithm 1 of [Pusa16]_
+    theta : numpy.ndarray
+        Complex poles :math:`\theta` of the rational approximation
+    alpha0 : float
+        Limit of the approximation at infinity
 
     """
 
-    theta_r = np.array([-4.465731934165702e+1, -5.284616241568964e+0,
-                        -8.867715667624458e+0, +3.493013124279215e+0,
-                        +1.564102508858634e+1, +1.742097597385893e+1,
-                        -2.834466755180654e+1, +1.661569367939544e+1,
-                        +8.011836167974721e+0, -2.056267541998229e+0,
-                        +1.449208170441839e+1, +1.853807176907916e+1,
-                        +9.932562704505182e+0, -2.244223871767187e+1,
-                        +8.590014121680897e-1, -1.286192925744479e+1,
-                        +1.164596909542055e+1, +1.806076684783089e+1,
-                        +5.870672154659249e+0, -3.542938819659747e+1,
-                        +1.901323489060250e+1, +1.885508331552577e+1,
-                        -1.734689708174982e+1, +1.316284237125190e+1])
-    theta_i = np.array([+6.233225190695437e+1, +4.057499381311059e+1,
-                        +4.325515754166724e+1, +3.281615453173585e+1,
-                        +1.558061616372237e+1, +1.076629305714420e+1,
-                        +5.492841024648724e+1, +1.316994930024688e+1,
-                        +2.780232111309410e+1, +3.794824788914354e+1,
-                        +1.799988210051809e+1, +5.974332563100539e+0,
-                        +2.532823409972962e+1, +5.179633600312162e+1,
-                        +3.536456194294350e+1, +4.600304902833652e+1,
-                        +2.287153304140217e+1, +8.368200580099821e+0,
-                        +3.029700159040121e+1, +5.834381701800013e+1,
-                        +1.194282058271408e+0, +3.583428564427879e+0,
-                        +4.883941101108207e+1, +2.042951874827759e+1])
+    theta_r = np.array([
+        -4.465731934165702e+1, -5.284616241568964e+0,
+        -8.867715667624458e+0, +3.493013124279215e+0,
+        +1.564102508858634e+1, +1.742097597385893e+1,
+        -2.834466755180654e+1, +1.661569367939544e+1,
+        +8.011836167974721e+0, -2.056267541998229e+0,
+        +1.449208170441839e+1, +1.853807176907916e+1,
+        +9.932562704505182e+0, -2.244223871767187e+1,
+        +8.590014121680897e-1, -1.286192925744479e+1,
+        +1.164596909542055e+1, +1.806076684783089e+1,
+        +5.870672154659249e+0, -3.542938819659747e+1,
+        +1.901323489060250e+1, +1.885508331552577e+1,
+        -1.734689708174982e+1, +1.316284237125190e+1])
+
+    theta_i = np.array([
+        +6.233225190695437e+1, +4.057499381311059e+1,
+        +4.325515754166724e+1, +3.281615453173585e+1,
+        +1.558061616372237e+1, +1.076629305714420e+1,
+        +5.492841024648724e+1, +1.316994930024688e+1,
+        +2.780232111309410e+1, +3.794824788914354e+1,
+        +1.799988210051809e+1, +5.974332563100539e+0,
+        +2.532823409972962e+1, +5.179633600312162e+1,
+        +3.536456194294350e+1, +4.600304902833652e+1,
+        +2.287153304140217e+1, +8.368200580099821e+0,
+        +3.029700159040121e+1, +5.834381701800013e+1,
+        +1.194282058271408e+0, +3.583428564427879e+0,
+        +4.883941101108207e+1, +2.042951874827759e+1])
+
     theta = np.array(theta_r + theta_i * 1j, dtype=np.complex128)
 
-    alpha_r = np.array([+6.387380733878774e+2, +1.909896179065730e+2,
-                        +4.236195226571914e+2, +4.645770595258726e+2,
-                        +7.765163276752433e+2, +1.907115136768522e+3,
-                        +2.909892685603256e+3, +1.944772206620450e+2,
-                        +1.382799786972332e+5, +5.628442079602433e+3,
-                        +2.151681283794220e+2, +1.324720240514420e+3,
-                        +1.617548476343347e+4, +1.112729040439685e+2,
-                        +1.074624783191125e+2, +8.835727765158191e+1,
-                        +9.354078136054179e+1, +9.418142823531573e+1,
-                        +1.040012390717851e+2, +6.861882624343235e+1,
-                        +8.766654491283722e+1, +1.056007619389650e+2,
-                        +7.738987569039419e+1, +1.041366366475571e+2])
-    alpha_i = np.array([-6.743912502859256e+2, -3.973203432721332e+2,
-                        -2.041233768918671e+3, -1.652917287299683e+3,
-                        -1.783617639907328e+4, -5.887068595142284e+4,
-                        -9.953255345514560e+3, -1.427131226068449e+3,
-                        -3.256885197214938e+6, -2.924284515884309e+4,
-                        -1.121774011188224e+3, -6.370088443140973e+4,
-                        -1.008798413156542e+6, -8.837109731680418e+1,
-                        -1.457246116408180e+2, -6.388286188419360e+1,
-                        -2.195424319460237e+2, -6.719055740098035e+2,
-                        -1.693747595553868e+2, -1.177598523430493e+1,
-                        -4.596464999363902e+3, -1.738294585524067e+3,
-                        -4.311715386228984e+1, -2.777743732451969e+2])
+    alpha_r = np.array([
+        +6.387380733878774e+2, +1.909896179065730e+2,
+        +4.236195226571914e+2, +4.645770595258726e+2,
+        +7.765163276752433e+2, +1.907115136768522e+3,
+        +2.909892685603256e+3, +1.944772206620450e+2,
+        +1.382799786972332e+5, +5.628442079602433e+3,
+        +2.151681283794220e+2, +1.324720240514420e+3,
+        +1.617548476343347e+4, +1.112729040439685e+2,
+        +1.074624783191125e+2, +8.835727765158191e+1,
+        +9.354078136054179e+1, +9.418142823531573e+1,
+        +1.040012390717851e+2, +6.861882624343235e+1,
+        +8.766654491283722e+1, +1.056007619389650e+2,
+        +7.738987569039419e+1, +1.041366366475571e+2])
+
+    alpha_i = np.array([
+        -6.743912502859256e+2, -3.973203432721332e+2,
+        -2.041233768918671e+3, -1.652917287299683e+3,
+        -1.783617639907328e+4, -5.887068595142284e+4,
+        -9.953255345514560e+3, -1.427131226068449e+3,
+        -3.256885197214938e+6, -2.924284515884309e+4,
+        -1.121774011188224e+3, -6.370088443140973e+4,
+        -1.008798413156542e+6, -8.837109731680418e+1,
+        -1.457246116408180e+2, -6.388286188419360e+1,
+        -2.195424319460237e+2, -6.719055740098035e+2,
+        -1.693747595553868e+2, -1.177598523430493e+1,
+        -4.596464999363902e+3, -1.738294585524067e+3,
+        -4.311715386228984e+1, -2.777743732451969e+2])
+
     alpha = np.array(alpha_r + alpha_i * 1j, dtype=np.complex128)
-    n = A.shape[0]
+
+    del theta_i, theta_r, alpha_r, alpha_i
 
     alpha0 = 2.258038182743983e-47
 
-    k = 24
+    def __call__(self, A, n0, dt):
+        """Solve using 48th order IPF CRAM [Pusa16]_
 
-    y = np.array(n0, dtype=np.float64)
-    for l in range(k):
-        y = 2.0*np.real(alpha[l]*sla.spsolve(A*dt - theta[l]*sp.eye(n), y)) + y
+        Parameters
+        ----------
+        A : scipy.sparse.csr_matrix
+            Sparse transmutation matrix ``A[j, i]`` desribing rates at
+            which isotope ``i`` transmutes to isotope ``j``
+        n0 : numpy.ndarray
+            Initial compositions, typically given in number of atoms in some
+            material or an atom density
+        dt : float
+            Time [s] of the specific interval to be solved
 
-    y *= alpha0
-    return y
+        Returns
+        -------
+        numpy.ndarray
+            Final compositions after ``dt``
+
+        """
+        return super().__call__(A, n0, dt)
+
+
+CRAM16 = Cram16Solver().__call__
+CRAM48 = Cram48Solver().__call__

From df387e87dc069836a80cf62ed4839369fe66ca4c Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Thu, 19 Sep 2019 16:51:38 -0500
Subject: [PATCH 2/3] Document concrete and abstract depletion solvers

---
 docs/source/pythonapi/deplete.rst | 18 ++++++++++++++++--
 1 file changed, 16 insertions(+), 2 deletions(-)

diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst
index e18afd920..03f5e6b61 100644
--- a/docs/source/pythonapi/deplete.rst
+++ b/docs/source/pythonapi/deplete.rst
@@ -139,6 +139,18 @@ The following functions are used to solve the depletion equations, with
    cram.deplete
    cram.timed_deplete
 
+
+:func:`cram.CRAM16` and :func:`cram.CRAM48` are aliases to the ``__call__``
+methods for the following classes
+
+.. autosummary::
+   :toctree: generated
+   :nosignatures:
+   :template: mycallable.rst
+
+   cram.Cram16Solver
+   cram.Cram48Solver
+
 The following classes are used to help the :class:`openmc.deplete.Operator`
 compute quantities like effective fission yields, reaction rates, and
 total system energy.
@@ -185,8 +197,8 @@ OpenMC simulations back on to the :class:`abc.TransportOperator`
    abc.ReactionRateHelper
    abc.TalliedFissionYieldHelper
 
-Custom integrators can be developed by subclassing from the following abstract
-base classes:
+Custom integrators or depletion solvers can be developed by subclassing from
+the following abstract base classes:
 
 .. autosummary::
    :toctree: generated
@@ -195,3 +207,5 @@ base classes:
 
    abc.Integrator
    abc.SIIntegrator
+   abc.DepSystemSolver
+   abc.IPFCramSolver

From d112eddbc39470b9df8d2d5b8797d69e7d1002c1 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Thu, 26 Sep 2019 20:26:13 -0400
Subject: [PATCH 3/3] openmc.deplete.cram.IPFCramSolver accepts coefficients at
 init

Moved from deplete.abc into openmc.deplete.cram. CRAM16 and CRAM48
are aliases to __call__ methods for two instances of
IPFCramSolver: Cram16Solver and Cram48Solver, created using
16th and 48th order coefficients
---
 docs/source/pythonapi/deplete.rst |  12 --
 openmc/deplete/abc.py             |  80 +---------
 openmc/deplete/cram.py            | 247 +++++++++++++++---------------
 3 files changed, 125 insertions(+), 214 deletions(-)

diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst
index 03f5e6b61..82143e14a 100644
--- a/docs/source/pythonapi/deplete.rst
+++ b/docs/source/pythonapi/deplete.rst
@@ -139,18 +139,6 @@ The following functions are used to solve the depletion equations, with
    cram.deplete
    cram.timed_deplete
 
-
-:func:`cram.CRAM16` and :func:`cram.CRAM48` are aliases to the ``__call__``
-methods for the following classes
-
-.. autosummary::
-   :toctree: generated
-   :nosignatures:
-   :template: mycallable.rst
-
-   cram.Cram16Solver
-   cram.Cram48Solver
-
 The following classes are used to help the :class:`openmc.deplete.Operator`
 compute quantities like effective fission yields, reaction rates, and
 total system energy.
diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py
index b42aa3e74..200552870 100644
--- a/openmc/deplete/abc.py
+++ b/openmc/deplete/abc.py
@@ -14,9 +14,7 @@ from copy import deepcopy
 from warnings import warn
 from numbers import Real, Integral
 
-from numpy import nonzero, empty, asarray, float64, real
-import scipy.sparse as sp
-import scipy.sparse.linalg as sla
+from numpy import nonzero, empty, asarray
 from uncertainties import ufloat
 
 from openmc.data import DataLibrary, JOULE_PER_EV
@@ -894,79 +892,3 @@ class DepSystemSolver(ABC):
             to ``n0``.
 
         """
-
-
-class IPFCramSolver(DepSystemSolver):
-    r"""Abstract class that implements the IPF form of CRAM
-
-    Provides a :meth:`__call__` that utilizes an incomplete
-    partial factorization (IPF) for the Chebyshev Rational Approximation
-    Method (CRAM) [Pusa16]_
-
-    Concrete subclasses must provide two complex vectors :attr:`alpha`
-    and :attr:`theta` that make up the coefficients of the decompostion.
-    Vectors are expected to be of equal length ``N``.
-    Subclases are also expected to provide a coefficient :attr:`alpha0`
-    used in the final scaling step.
-
-    Attributes
-    ----------
-    alpha : numpy.ndarray
-        Complex residues of poles :attr:`theta` in the incomplete partial
-        factorization. Denoted as :math:`\tilde{\alpha}`
-    theta : numpy.ndarray
-        Complex poles :math:`\theta` of the rational approximation
-    alpha0 : float
-        Limit of the approximation at infinity
-
-    References
-    ----------
-
-    .. [Pusa16] M. Pusa, "Higher-Order Chebyshev Rational Approximation
-       Method and Application to Burnup Equations," Nuclear Science And
-       Engineering, 182:3,297-318
-       `DOI: 10.13182/NSE15-26 <https://doi.org/10.13182/NSE15-26>`_
-
-    """
-
-    @property
-    @abstractmethod
-    def alpha(self):
-        pass
-
-    @property
-    @abstractmethod
-    def theta(self):
-        pass
-
-    @property
-    @abstractmethod
-    def alpha0(self):
-        pass
-
-    def __call__(self, A, n0, dt):
-        """Solve depletion equations using IPF CRAM
-
-        Parameters
-        ----------
-        A : scipy.sparse.csr_matrix
-            Sparse transmutation matrix ``A[j, i]`` desribing rates at
-            which isotope ``i`` transmutes to isotope ``j``
-        n0 : numpy.ndarray
-            Initial compositions, typically given in number of atoms in some
-            material or an atom density
-        dt : float
-            Time [s] of the specific interval to be solved
-
-        Returns
-        -------
-        numpy.ndarray
-            Final compositions after ``dt``
-
-        """
-        A = sp.csr_matrix(A * dt, dtype=float64)
-        y = asarray(n0, dtype=float64)
-        ident = sp.eye(A.shape[0])
-        for alpha, theta in zip(self.alpha, self.theta):
-            y += 2*real(alpha*sla.spsolve(A - theta*ident, y))
-        return y * self.alpha0
diff --git a/openmc/deplete/cram.py b/openmc/deplete/cram.py
index 4fdba87bf..f4a847907 100644
--- a/openmc/deplete/cram.py
+++ b/openmc/deplete/cram.py
@@ -3,13 +3,17 @@
 Implements two different forms of CRAM for use in openmc.deplete.
 """
 
+import numbers
 from itertools import repeat
 from multiprocessing import Pool
 import time
 
 import numpy as np
+import scipy.sparse as sp
+import scipy.sparse.linalg as sla
 
-from .abc import IPFCramSolver
+from openmc.checkvalue import check_type, check_length
+from .abc import DepSystemSolver
 
 __all__ = [
     "deplete", "timed_deplete", "CRAM16", "CRAM48",
@@ -83,50 +87,51 @@ def timed_deplete(*args, **kwargs):
     return time.time() - start, results
 
 
-class Cram16Solver(IPFCramSolver):
-    r"""Solver implementing the 16th order IPF CRAM
+class IPFCramSolver(DepSystemSolver):
+    r"""Class for solving depletion systems with IPF Cram
 
-    Coefficients to :attr:`alpha`, :attr:`theta`, and :attr:`alpha0`
-    are from Table A.IV in [Pusa16]_.
+    Provides a :meth:`__call__` that utilizes an incomplete
+    partial factorization (IPF) for the Chebyshev Rational Approximation
+    Method (CRAM)
+
+    M. Pusa, "Higher-Order Chebyshev Rational Approximation
+    Method and Application to Burnup Equations," Nuclear Science And
+    Engineering, 182:3,297-318
+    `DOI: 10.13182/NSE15-26 <https://doi.org/10.13182/NSE15-26>`_
+
+    Parameters
+    ----------
+    alpha : numpy.ndarray
+        Complex residues of poles used in the factorization. Must be a
+        vector with even number of items.
+    theta : numpy.ndarray
+        Complex poles. Must have an equal size as ``alpha``.
+    alpha0 : float
+        Limit of the approximation at infinity
 
     Attributes
     ----------
     alpha : numpy.ndarray
         Complex residues of poles :attr:`theta` in the incomplete partial
         factorization. Denoted as :math:`\tilde{\alpha}`
-        in Algorithm 1 of [Pusa16]_
     theta : numpy.ndarray
         Complex poles :math:`\theta` of the rational approximation
     alpha0 : float
         Limit of the approximation at infinity
 
     """
-    alpha = np.array([
-        +5.464930576870210e+3 - 3.797983575308356e+4j,
-        +9.045112476907548e+1 - 1.115537522430261e+3j,
-        +2.344818070467641e+2 - 4.228020157070496e+2j,
-        +9.453304067358312e+1 - 2.951294291446048e+2j,
-        +7.283792954673409e+2 - 1.205646080220011e+5j,
-        +3.648229059594851e+1 - 1.155509621409682e+2j,
-        +2.547321630156819e+1 - 2.639500283021502e+1j,
-        +2.394538338734709e+1 - 5.650522971778156e+0j],
-        dtype=np.complex128)
 
-    theta = np.array([
-        +3.509103608414918 + 8.436198985884374j,
-        +5.948152268951177 + 3.587457362018322j,
-        -5.264971343442647 + 16.22022147316793j,
-        +1.419375897185666 + 10.92536348449672j,
-        +6.416177699099435 + 1.194122393370139j,
-        +4.993174737717997 + 5.996881713603942j,
-        -1.413928462488886 + 13.49772569889275j,
-        -10.84391707869699 + 19.27744616718165j],
-        dtype=np.complex128)
-
-    alpha0 = 2.124853710495224e-16
+    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)
+        check_type("alpha0", alpha0, numbers.Real)
+        self.alpha = alpha
+        self.theta = theta
+        self.alpha0 = alpha0
 
     def __call__(self, A, n0, dt):
-        """Solve using 16th order IPF CRAM [Pusa16]_
+        """Solve depletion equations using IPF CRAM
 
         Parameters
         ----------
@@ -145,114 +150,110 @@ class Cram16Solver(IPFCramSolver):
             Final compositions after ``dt``
 
         """
-        return super().__call__(A, n0, dt)
+        A = sp.csr_matrix(A * dt, dtype=np.float64)
+        y = np.asarray(n0, dtype=np.float64)
+        ident = sp.eye(A.shape[0])
+        for alpha, theta in zip(self.alpha, self.theta):
+            y += 2*np.real(alpha*sla.spsolve(A - theta*ident, y))
+        return y * self.alpha0
 
 
-class Cram48Solver(IPFCramSolver):
-    r"""Solver implementing the 48th order IPF CRAM
+# Coefficients for IPF Cram 16
+c16_alpha = np.array([
+    +5.464930576870210e+3 - 3.797983575308356e+4j,
+    +9.045112476907548e+1 - 1.115537522430261e+3j,
+    +2.344818070467641e+2 - 4.228020157070496e+2j,
+    +9.453304067358312e+1 - 2.951294291446048e+2j,
+    +7.283792954673409e+2 - 1.205646080220011e+5j,
+    +3.648229059594851e+1 - 1.155509621409682e+2j,
+    +2.547321630156819e+1 - 2.639500283021502e+1j,
+    +2.394538338734709e+1 - 5.650522971778156e+0j],
+    dtype=np.complex128)
 
-    Coefficients to :attr:`alpha`, :attr:`theta`, and :attr:`alpha0`
-    are from Table A.XII in [Pusa16]_.
+c16_theta = np.array([
+    +3.509103608414918 + 8.436198985884374j,
+    +5.948152268951177 + 3.587457362018322j,
+    -5.264971343442647 + 16.22022147316793j,
+    +1.419375897185666 + 10.92536348449672j,
+    +6.416177699099435 + 1.194122393370139j,
+    +4.993174737717997 + 5.996881713603942j,
+    -1.413928462488886 + 13.49772569889275j,
+    -10.84391707869699 + 19.27744616718165j],
+    dtype=np.complex128)
 
-    Attributes
-    ----------
-    alpha : numpy.ndarray
-        Complex residues of poles :attr:`theta` in the incomplete partial
-        factorization. Denoted as :math:`\tilde{\alpha}`
-        in Algorithm 1 of [Pusa16]_
-    theta : numpy.ndarray
-        Complex poles :math:`\theta` of the rational approximation
-    alpha0 : float
-        Limit of the approximation at infinity
+c16_alpha0 = 2.124853710495224e-16
+Cram16Solver = IPFCramSolver(c16_alpha, c16_theta, c16_alpha0)
+CRAM16 = Cram16Solver.__call__
 
-    """
+del c16_alpha, c16_alpha0, c16_theta
 
-    theta_r = np.array([
-        -4.465731934165702e+1, -5.284616241568964e+0,
-        -8.867715667624458e+0, +3.493013124279215e+0,
-        +1.564102508858634e+1, +1.742097597385893e+1,
-        -2.834466755180654e+1, +1.661569367939544e+1,
-        +8.011836167974721e+0, -2.056267541998229e+0,
-        +1.449208170441839e+1, +1.853807176907916e+1,
-        +9.932562704505182e+0, -2.244223871767187e+1,
-        +8.590014121680897e-1, -1.286192925744479e+1,
-        +1.164596909542055e+1, +1.806076684783089e+1,
-        +5.870672154659249e+0, -3.542938819659747e+1,
-        +1.901323489060250e+1, +1.885508331552577e+1,
-        -1.734689708174982e+1, +1.316284237125190e+1])
+# Coefficients for 48th order IPF Cram
 
-    theta_i = np.array([
-        +6.233225190695437e+1, +4.057499381311059e+1,
-        +4.325515754166724e+1, +3.281615453173585e+1,
-        +1.558061616372237e+1, +1.076629305714420e+1,
-        +5.492841024648724e+1, +1.316994930024688e+1,
-        +2.780232111309410e+1, +3.794824788914354e+1,
-        +1.799988210051809e+1, +5.974332563100539e+0,
-        +2.532823409972962e+1, +5.179633600312162e+1,
-        +3.536456194294350e+1, +4.600304902833652e+1,
-        +2.287153304140217e+1, +8.368200580099821e+0,
-        +3.029700159040121e+1, +5.834381701800013e+1,
-        +1.194282058271408e+0, +3.583428564427879e+0,
-        +4.883941101108207e+1, +2.042951874827759e+1])
+theta_r = np.array([
+    -4.465731934165702e+1, -5.284616241568964e+0,
+    -8.867715667624458e+0, +3.493013124279215e+0,
+    +1.564102508858634e+1, +1.742097597385893e+1,
+    -2.834466755180654e+1, +1.661569367939544e+1,
+    +8.011836167974721e+0, -2.056267541998229e+0,
+    +1.449208170441839e+1, +1.853807176907916e+1,
+    +9.932562704505182e+0, -2.244223871767187e+1,
+    +8.590014121680897e-1, -1.286192925744479e+1,
+    +1.164596909542055e+1, +1.806076684783089e+1,
+    +5.870672154659249e+0, -3.542938819659747e+1,
+    +1.901323489060250e+1, +1.885508331552577e+1,
+    -1.734689708174982e+1, +1.316284237125190e+1])
 
-    theta = np.array(theta_r + theta_i * 1j, dtype=np.complex128)
+theta_i = np.array([
+    +6.233225190695437e+1, +4.057499381311059e+1,
+    +4.325515754166724e+1, +3.281615453173585e+1,
+    +1.558061616372237e+1, +1.076629305714420e+1,
+    +5.492841024648724e+1, +1.316994930024688e+1,
+    +2.780232111309410e+1, +3.794824788914354e+1,
+    +1.799988210051809e+1, +5.974332563100539e+0,
+    +2.532823409972962e+1, +5.179633600312162e+1,
+    +3.536456194294350e+1, +4.600304902833652e+1,
+    +2.287153304140217e+1, +8.368200580099821e+0,
+    +3.029700159040121e+1, +5.834381701800013e+1,
+    +1.194282058271408e+0, +3.583428564427879e+0,
+    +4.883941101108207e+1, +2.042951874827759e+1])
 
-    alpha_r = np.array([
-        +6.387380733878774e+2, +1.909896179065730e+2,
-        +4.236195226571914e+2, +4.645770595258726e+2,
-        +7.765163276752433e+2, +1.907115136768522e+3,
-        +2.909892685603256e+3, +1.944772206620450e+2,
-        +1.382799786972332e+5, +5.628442079602433e+3,
-        +2.151681283794220e+2, +1.324720240514420e+3,
-        +1.617548476343347e+4, +1.112729040439685e+2,
-        +1.074624783191125e+2, +8.835727765158191e+1,
-        +9.354078136054179e+1, +9.418142823531573e+1,
-        +1.040012390717851e+2, +6.861882624343235e+1,
-        +8.766654491283722e+1, +1.056007619389650e+2,
-        +7.738987569039419e+1, +1.041366366475571e+2])
+c48_theta = np.array(theta_r + theta_i * 1j, dtype=np.complex128)
 
-    alpha_i = np.array([
-        -6.743912502859256e+2, -3.973203432721332e+2,
-        -2.041233768918671e+3, -1.652917287299683e+3,
-        -1.783617639907328e+4, -5.887068595142284e+4,
-        -9.953255345514560e+3, -1.427131226068449e+3,
-        -3.256885197214938e+6, -2.924284515884309e+4,
-        -1.121774011188224e+3, -6.370088443140973e+4,
-        -1.008798413156542e+6, -8.837109731680418e+1,
-        -1.457246116408180e+2, -6.388286188419360e+1,
-        -2.195424319460237e+2, -6.719055740098035e+2,
-        -1.693747595553868e+2, -1.177598523430493e+1,
-        -4.596464999363902e+3, -1.738294585524067e+3,
-        -4.311715386228984e+1, -2.777743732451969e+2])
+alpha_r = np.array([
+    +6.387380733878774e+2, +1.909896179065730e+2,
+    +4.236195226571914e+2, +4.645770595258726e+2,
+    +7.765163276752433e+2, +1.907115136768522e+3,
+    +2.909892685603256e+3, +1.944772206620450e+2,
+    +1.382799786972332e+5, +5.628442079602433e+3,
+    +2.151681283794220e+2, +1.324720240514420e+3,
+    +1.617548476343347e+4, +1.112729040439685e+2,
+    +1.074624783191125e+2, +8.835727765158191e+1,
+    +9.354078136054179e+1, +9.418142823531573e+1,
+    +1.040012390717851e+2, +6.861882624343235e+1,
+    +8.766654491283722e+1, +1.056007619389650e+2,
+    +7.738987569039419e+1, +1.041366366475571e+2])
 
-    alpha = np.array(alpha_r + alpha_i * 1j, dtype=np.complex128)
+alpha_i = np.array([
+    -6.743912502859256e+2, -3.973203432721332e+2,
+    -2.041233768918671e+3, -1.652917287299683e+3,
+    -1.783617639907328e+4, -5.887068595142284e+4,
+    -9.953255345514560e+3, -1.427131226068449e+3,
+    -3.256885197214938e+6, -2.924284515884309e+4,
+    -1.121774011188224e+3, -6.370088443140973e+4,
+    -1.008798413156542e+6, -8.837109731680418e+1,
+    -1.457246116408180e+2, -6.388286188419360e+1,
+    -2.195424319460237e+2, -6.719055740098035e+2,
+    -1.693747595553868e+2, -1.177598523430493e+1,
+    -4.596464999363902e+3, -1.738294585524067e+3,
+    -4.311715386228984e+1, -2.777743732451969e+2])
 
-    del theta_i, theta_r, alpha_r, alpha_i
+c48_alpha = np.array(alpha_r + alpha_i * 1j, dtype=np.complex128)
 
-    alpha0 = 2.258038182743983e-47
+c48_alpha0 = 2.258038182743983e-47
 
-    def __call__(self, A, n0, dt):
-        """Solve using 48th order IPF CRAM [Pusa16]_
+Cram48Solver = IPFCramSolver(c48_alpha, c48_theta, c48_alpha0)
 
-        Parameters
-        ----------
-        A : scipy.sparse.csr_matrix
-            Sparse transmutation matrix ``A[j, i]`` desribing rates at
-            which isotope ``i`` transmutes to isotope ``j``
-        n0 : numpy.ndarray
-            Initial compositions, typically given in number of atoms in some
-            material or an atom density
-        dt : float
-            Time [s] of the specific interval to be solved
+del c48_alpha, c48_alpha0, c48_theta, alpha_r, alpha_i, theta_r, theta_i
 
-        Returns
-        -------
-        numpy.ndarray
-            Final compositions after ``dt``
+CRAM48 = Cram48Solver.__call__
 
-        """
-        return super().__call__(A, n0, dt)
-
-
-CRAM16 = Cram16Solver().__call__
-CRAM48 = Cram48Solver().__call__