From d112eddbc39470b9df8d2d5b8797d69e7d1002c1 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Thu, 26 Sep 2019 20:26:13 -0400 Subject: [PATCH] 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 `_ - - """ - - @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 `_ + + 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__