crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-26 13:15:39 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge pull request #1356 from drewejohnson/dep-solver-class

Browse source 
Provide class-based depletion solvers
This commit is contained in:
Paul Romano 2019-09-30 13:23:29 -05:00 • committed by GitHub
commit ed7123f4e7
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
3 changed files with 204 additions and 133 deletions
Download patch file Download diff file Expand all files Collapse all files

6
docs/source/pythonapi/deplete.rst
View file

@ -185,8 +185,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 +195,5 @@ base classes:
abc.Integrator
abc.SIIntegrator
abc.DepSystemSolver
abc.IPFCramSolver

37
openmc/deplete/abc.py
View file

@ -855,3 +855,40 @@ 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``.
"""

294
openmc/deplete/cram.py
View file

@ -3,6 +3,7 @@
Implements two different forms of CRAM for use in openmc.deplete.
"""
import numbers
from itertools import repeat
from multiprocessing import Pool
import time
@ -11,7 +12,12 @@ import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as sla
from . import comm
from openmc.checkvalue import check_type, check_length
from .abc import DepSystemSolver
__all__ = [
"deplete", "timed_deplete", "CRAM16", "CRAM48",
"Cram16Solver", "Cram48Solver"]
def deplete(chain, x, rates, dt, matrix_func=None):
@ -81,147 +87,173 @@ def timed_deplete(*args, **kwargs):
return time.time() - start, results
def CRAM16(A, n0, dt):
"""Chebyshev Rational Approximation Method, order 16
class IPFCramSolver(DepSystemSolver):
r"""Class for solving depletion systems with 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.
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
----------
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.
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
Returns
-------
numpy.array
Results of the matrix exponent.
"""
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]
alpha0 = 2.124853710495224e-16
k = 8
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
y *= alpha0
return y
def CRAM48(A, n0, dt):
"""Chebyshev Rational Approximation Method, order 48
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.
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}`
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 = np.array(theta_r + theta_i * 1j, dtype=np.complex128)
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
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]
def __call__(self, A, n0, dt):
"""Solve depletion equations using IPF CRAM
alpha0 = 2.258038182743983e-47
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
k = 24
Returns
-------
numpy.ndarray
Final compositions after ``dt``
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
"""
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
# 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)
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)
c16_alpha0 = 2.124853710495224e-16
Cram16Solver = IPFCramSolver(c16_alpha, c16_theta, c16_alpha0)
CRAM16 = Cram16Solver.__call__
del c16_alpha, c16_alpha0, c16_theta
# Coefficients for 48th order IPF Cram
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])
c48_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])
c48_alpha = np.array(alpha_r + alpha_i * 1j, dtype=np.complex128)
c48_alpha0 = 2.258038182743983e-47
Cram48Solver = IPFCramSolver(c48_alpha, c48_theta, c48_alpha0)
del c48_alpha, c48_alpha0, c48_theta, alpha_r, alpha_i, theta_r, theta_i
CRAM48 = Cram48Solver.__call__
y *= alpha0
return y