diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 565b9e2dcf..b42aa3e745 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 `_ + + """ + + @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 31049bb724..4fdba87bfb 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) - `_ 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) - `_ 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__