crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 05:35:49 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge a349be698e into 05d01274a7

Browse source
This commit is contained in:
Ethan Peterson 2026-07-20 13:34:32 +01:00 • committed by GitHub
commit 2d779d14a4
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
15 changed files with 1145 additions and 261 deletions
Download patch file Download diff file Expand all files Collapse all files

2
CMakeLists.txt
View file

@ -424,6 +424,8 @@ list(APPEND libopenmc_SOURCES
src/settings.cpp
src/simulation.cpp
src/source.cpp
src/sparse_matrix.cpp
src/bateman_solvers.cpp
src/state_point.cpp
src/string_utils.cpp
src/summary.cpp

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

@ -153,13 +153,6 @@ data, such as number densities and reaction rates for each material.
The following class and functions are used to solve the depletion equations,
with :func:`cram.CRAM48` being the default.
.. autosummary::
:toctree: generated
:nosignatures:
:template: myintegrator.rst
cram.IPFCramSolver
.. autosummary::
:toctree: generated
:nosignatures:
@ -275,8 +268,8 @@ transport-independent depletion) back on to the :class:`abc.TransportOperator`
abc.FissionYieldHelper
abc.ReactionRateHelper
Custom integrators or depletion solvers can be developed by subclassing from
the following abstract base classes:
Custom integrators can be developed by subclassing from the following abstract
base classes:
.. autosummary::
:toctree: generated
@ -285,7 +278,6 @@ the following abstract base classes:
abc.Integrator
abc.SIIntegrator
abc.DepSystemSolver
R2S Automation
--------------

102
include/openmc/bateman_solvers.h Normal file
View file

@ -0,0 +1,102 @@
//! \file bateman_solvers.h
//! \brief Solvers for the Bateman depletion equations
#ifndef OPENMC_BATEMAN_SOLVERS_H
#define OPENMC_BATEMAN_SOLVERS_H
#include <complex>
#include "openmc/sparse_matrix.h"
#include "openmc/vector.h"
namespace openmc {
//==============================================================================
//! Numeric LU factorization values for a CRAM shifted linear system
//!
//! Stores the complex-valued L/U entries for (dt*A - theta*I) corresponding
//! to a previously computed SymbolicLUFactorization. The U diagonal is stored
//! both in the CSC values array and as precomputed reciprocals for faster back
//! substitution.
//==============================================================================
struct NumericLUFactorization {
vector<std::complex<double>> l_data;
vector<std::complex<double>> u_data;
vector<std::complex<double>> u_diag_inv;
};
//! Numeric left-looking LU factorization of the CRAM shifted operator
//! (dt*A - theta*I) without explicitly forming it. Uses no pivoting; this is
//! safe because A is Metzler (non-negative off-diagonal) and |Im(theta)| >=
//! 1.194 for every CRAM pole, so every pivot satisfies |u_jj| >= |Im(theta)|.
void numeric_factorize_cram(const CSCMatrix& A, double dt,
std::complex<double> theta, const SymbolicLUFactorization& symbolic,
NumericLUFactorization& numeric, vector<std::complex<double>>& work);
//! Solve LUx = b using a symbolic LU pattern and matching numeric values.
void triangular_solve_lu(const vector<double>& b,
const SymbolicLUFactorization& symbolic,
const NumericLUFactorization& numeric, vector<std::complex<double>>& x);
//==============================================================================
//! Abstract base class for Bateman equation solvers
//!
//! Solves dN/dt = A*N over an interval dt, given initial composition N(0).
//==============================================================================
class BatemanSolver {
public:
virtual ~BatemanSolver() = default;
//! Solve the Bateman equations
//! \param A Sparse transmutation matrix (n x n)
//! \param n0 Initial atom densities [n]
//! \param dt Time interval [s]
//! \param substeps Number of substeps to use within dt
//! \return Final atom densities [n]
virtual vector<double> solve(const CSCMatrix& A, const vector<double>& n0,
double dt, int substeps = 1) = 0;
};
//==============================================================================
//! IPF CRAM solver for general transmutation matrices
//!
//! Implements the Incomplete Partial Fraction form of the Chebyshev
//! Rational Approximation Method (CRAM), as described in:
//! M. Pusa, "Higher-Order Chebyshev Rational Approximation Method and
//! Application to Burnup Equations," Nucl. Sci. Eng., 182:3, 297-318 (2016).
//!
//! The numeric factorization uses left-looking column LU without
//! pivoting. Each call to solve() performs a symbolic factorization
//! to compute L/U sparsity patterns, then reuses those patterns for
//! all pole solves within that call. Pivoting is unnecessary
//! because the transmutation matrix is Metzler (non-negative off-diagonal)
//! and the CRAM poles have nonzero imaginary parts (|Im(theta)| >= 1.194).
//! For any real Metzler matrix R and complex shift theta with Im(theta) != 0,
//! unpivoted Gaussian elimination on (R - theta*I) produces pivots u_jj
//! satisfying |u_jj| >= |Im(theta)|, guaranteeing non-singular factorization.
//! Since pivoting is not needed, the L/U sparsity patterns are deterministic
//! and identical across all poles, allowing a single symbolic phase to serve
//! all 24 (CRAM48) or 8 (CRAM16) linear solves.
//==============================================================================
class IPFCramSolver : public BatemanSolver {
public:
explicit IPFCramSolver(int order = 48);
//! Solve using full LU factorization (general transmutation matrix).
vector<double> solve(const CSCMatrix& A, const vector<double>& n0, double dt,
int substeps = 1) override;
private:
// --- CRAM coefficients (non-owning views of the static tables) ---
int n_poles_; //!< Number of poles (k/2)
const std::complex<double>* alpha_; //!< Residues [n_poles]
const std::complex<double>* theta_; //!< Poles [n_poles]
double alpha0_; //!< Limit at infinity
};
} // namespace openmc
#endif // OPENMC_BATEMAN_SOLVERS_H

22
include/openmc/capi.h
View file

@ -352,6 +352,28 @@ int openmc_properties_export(const char* filename);
// \return Error code
int openmc_properties_import(const char* filename);
//! Solve a Bateman system dN/dt = A*N over interval dt using CRAM.
//!
//! The transmutation matrix A is provided in CSC (Compressed Sparse Column)
//! format. The result is written to a caller-allocated output buffer.
//!
//! All pointer arguments must be non-null. CSC row indices within each
//! column must be sorted in ascending order.
//!
//! \param[in] n Matrix dimension (number of nuclides)
//! \param[in] indptr CSC column pointers [n+1]
//! \param[in] indices CSC row indices [nnz]
//! \param[in] data CSC nonzero values [nnz]
//! \param[in] n0 Initial atom densities [n]
//! \param[in] dt Time interval in seconds
//! \param[in] order CRAM approximation order (16 or 48)
//! \param[in] substeps Number of substeps to use within dt
//! \param[out] result Final atom densities [n]
//! \return Error code
int openmc_cram_solve(int n, const int* indptr, const int* indices,
const double* data, const double* n0, double dt, int order, int substeps,
double* result);
// Error codes
extern int OPENMC_E_UNASSIGNED;
extern int OPENMC_E_ALLOCATE;

94
include/openmc/sparse_matrix.h Normal file
View file

@ -0,0 +1,94 @@
//! \file sparse_matrix.h
//! \brief Compressed Sparse Column (CSC) sparsity pattern and matrix classes
#ifndef OPENMC_SPARSE_MATRIX_H
#define OPENMC_SPARSE_MATRIX_H
#include <utility> // for move
#include "openmc/vector.h"
namespace openmc {
//==============================================================================
//! CSC sparsity pattern (structure only, no values)
//!
//! Stores the compressed column structure of a square sparse matrix:
//! column pointers and row indices. Rows within each column are sorted in
//! ascending order.
//==============================================================================
class CSCPattern {
public:
// Constructors
CSCPattern() = default;
CSCPattern(int n, vector<int> indptr, vector<int> indices);
// Accessors
int n() const { return n_; }
int nnz() const { return static_cast<int>(indices_.size()); }
const vector<int>& indptr() const { return indptr_; }
const vector<int>& indices() const { return indices_; }
//! Return a new pattern with all diagonal entries forced present.
//! Existing entries (including any diagonals already present) are preserved.
CSCPattern with_diagonal() const;
//! Structural equality check. Two patterns are equal iff they have the same
//! dimension, identical column pointers, and identical row indices.
bool operator==(const CSCPattern& other) const;
bool operator!=(const CSCPattern& other) const { return !(*this == other); }
private:
int n_ {0}; //!< Matrix dimension
vector<int> indptr_; //!< Column pointers [n+1]
vector<int> indices_; //!< Row indices [nnz], sorted within each column
};
//==============================================================================
//! CSC sparse matrix with real (double) values
//!
//! Associates a double-precision value with each structural nonzero defined
//! by the underlying CSCPattern.
//==============================================================================
class CSCMatrix {
public:
// Constructors
CSCMatrix() = default;
CSCMatrix(
int n, vector<int> indptr, vector<int> indices, vector<double> data);
// Accessors
int n() const { return pattern_.n(); }
int nnz() const { return pattern_.nnz(); }
const CSCPattern& pattern() const { return pattern_; }
const vector<int>& indptr() const { return pattern_.indptr(); }
const vector<int>& indices() const { return pattern_.indices(); }
const vector<double>& data() const { return data_; }
private:
CSCPattern pattern_; //!< Structural pattern
vector<double> data_; //!< Values [nnz]
};
//==============================================================================
//! Symbolic LU factorization for CSC matrices
//!
//! Stores the shared structural information for left-looking column LU
//! factorization without pivoting.
//==============================================================================
struct SymbolicLUFactorization {
CSCPattern pattern;
CSCPattern l_pattern;
CSCPattern u_pattern;
};
//! Compute symbolic LU fill patterns for left-looking column LU without
//! pivoting.
SymbolicLUFactorization symbolic_factorize(CSCPattern pattern);
} // namespace openmc
#endif // OPENMC_SPARSE_MATRIX_H

43
openmc/deplete/abc.py
View file

@ -1,7 +1,7 @@
"""abc module.
This module contains Abstract Base Classes for implementing operator,
integrator, depletion system solver, and operator helper classes
integrator, and operator helper classes
"""
from __future__ import annotations
@ -37,7 +37,7 @@ from .keff_search_control import _KeffSearchControl
__all__ = [
"OperatorResult", "TransportOperator",
"ReactionRateHelper", "NormalizationHelper", "FissionYieldHelper",
"Integrator", "SIIntegrator", "DepSystemSolver", "add_params"]
"Integrator", "SIIntegrator", "add_params"]
def _normalize_timesteps(
@ -1448,42 +1448,3 @@ class SIIntegrator(Integrator):
self.operator.write_bos_data(self._i_res + len(self))
self.operator.finalize()
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, substeps=1):
"""Solve the linear system of equations for depletion
Parameters
----------
A : scipy.sparse.csc_array
Sparse transmutation matrix ``A[j, i]`` describing 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
substeps : int, optional
Number of substeps to use when the solver supports substepping.
Returns
-------
numpy.ndarray
Final compositions after ``dt``. Should be of identical shape
to ``n0``.
"""

241
openmc/deplete/cram.py
View file

@ -1,205 +1,82 @@
"""Chebyshev Rational Approximation Method module
"""Chebyshev Rational Approximation Method module.
Implements two different forms of CRAM for use in openmc.deplete.
Provides :func:`CRAM16` and :func:`CRAM48` solvers for use in
:mod:`openmc.deplete`, backed by the C++ ``IPFCramSolver`` implementation in
:mod:`openmc.lib.deplete`.
"""
from functools import partial
import numbers
import numpy as np
from scipy.sparse.linalg import spsolve, splu
from openmc.checkvalue import check_type, check_length, check_greater_than
from .abc import DepSystemSolver
from .._sparse_compat import csc_array, eye_array
from openmc.lib.deplete import cram_solve
from .._sparse_compat import csc_array
__all__ = ["CRAM16", "CRAM48", "Cram16Solver", "Cram48Solver", "IPFCramSolver"]
__all__ = ["CRAM16", "CRAM48"]
class IPFCramSolver(DepSystemSolver):
r"""CRAM depletion solver that uses incomplete partial factorization
def _cram_solve(A, n0, dt, order, substeps=1):
"""Single-material CRAM solve via the C++ ``IPFCramSolver`` backend."""
return cram_solve(
csc_array(A, dtype=np.float64), np.asarray(n0, dtype=np.float64),
float(dt), order=order, substeps=substeps)
Provides a :meth:`__call__` that utilizes an incomplete
partial factorization (IPF) for the Chebyshev Rational Approximation
Method (CRAM), as described in the following paper: M. Pusa, "`Higher-Order
Chebyshev Rational Approximation Method and Application to Burnup Equations
<https://doi.org/10.13182/NSE15-26>`_," Nucl. Sci. Eng., 182:3, 297-318.
When `substeps` > 1, the time interval is split into `substeps` identical
sub-intervals and LU factorizations are reused across them, as described
in: A. Isotalo and M. Pusa, "`Improving the Accuracy of the Chebyshev
Rational Approximation Method Using Substeps
def CRAM48(A, n0, dt, substeps=1):
r"""Solve depletion equations using 48th order IPF CRAM.
Implements the Incomplete Partial Fraction form of the Chebyshev
Rational Approximation Method (CRAM) as described in M. Pusa,
"`Higher-Order Chebyshev Rational Approximation Method and Application to
Burnup Equations <https://doi.org/10.13182/NSE15-26>`_," Nucl. Sci. Eng.,
182:3, 297-318. When ``substeps > 1`` the time interval is split into
equal sub-intervals and LU factorizations are reused across them as
described in A. Isotalo and M. Pusa, "`Improving the Accuracy of the
Chebyshev Rational Approximation Method Using Substeps
<https://doi.org/10.13182/NSE15-67>`_," Nucl. Sci. Eng., 183:1, 65-77.
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
A : scipy.sparse.csc_array
Sparse transmutation matrix ``A[j, i]`` describing 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 interval to be solved.
substeps : int, optional
Number of equal substeps to take within ``dt``. Defaults to 1.
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
Returns
-------
numpy.ndarray
Final compositions after ``dt``.
"""
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, substeps=1):
"""Solve depletion equations using IPF CRAM
Parameters
----------
A : scipy.sparse.csc_array
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
substeps : int, optional
Number of substeps per depletion interval.
Returns
-------
numpy.ndarray
Final compositions after ``dt``
"""
check_type("substeps", substeps, numbers.Integral)
check_greater_than("substeps", substeps, 0)
step_dt = dt if substeps == 1 else dt / substeps
A = step_dt * csc_array(A, dtype=np.float64)
ident = eye_array(A.shape[0], format='csc')
if substeps == 1:
solvers = [partial(spsolve, A - theta * ident) for theta in self.theta]
else:
# Pre-compute LU factorizations and reuse them across substeps.
solvers = [splu(A - theta * ident).solve for theta in self.theta]
y = n0.copy()
for _ in range(substeps):
for alpha, solve in zip(self.alpha, solvers):
y += 2 * np.real(alpha * solve(y))
y *= self.alpha0
return y
return _cram_solve(A, n0, dt, order=48, substeps=substeps)
# 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)
def CRAM16(A, n0, dt, substeps=1):
r"""Solve depletion equations using 16th order IPF CRAM.
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)
See :func:`CRAM48` for a description of the method and substep behavior.
c16_alpha0 = 2.124853710495224e-16
Cram16Solver = IPFCramSolver(c16_alpha, c16_theta, c16_alpha0)
CRAM16 = Cram16Solver.__call__
Parameters
----------
A : scipy.sparse.csc_array
Sparse transmutation matrix ``A[j, i]`` describing 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 interval to be solved.
substeps : int, optional
Number of equal substeps to take within ``dt``. Defaults to 1.
del c16_alpha, c16_alpha0, c16_theta
Returns
-------
numpy.ndarray
Final compositions after ``dt``.
# 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__
"""
return _cram_solve(A, n0, dt, order=16, substeps=substeps)

1
openmc/lib/__init__.py
View file

@ -66,6 +66,7 @@ from .math import *
from .plot import *
from .weight_windows import *
from .dagmc import *
from . import deplete
# Flag to denote whether or not openmc.lib.init has been called
# TODO: Establish and use a flag in the C++ code to represent the status of the

71
openmc/lib/deplete.py Normal file
View file

@ -0,0 +1,71 @@
"""Ctypes bindings for C++ depletion solvers."""
from ctypes import c_int, c_double
import numbers
import numpy as np
from numpy.ctypeslib import ndpointer
from .error import _error_handler
from . import _dll
_array_1d_int = ndpointer(dtype=np.int32, ndim=1, flags='CONTIGUOUS')
_array_1d_dbl = ndpointer(dtype=np.float64, ndim=1, flags='CONTIGUOUS')
# --- CRAM single-material solve ---
_dll.openmc_cram_solve.restype = c_int
_dll.openmc_cram_solve.errcheck = _error_handler
_dll.openmc_cram_solve.argtypes = [
c_int, # n
_array_1d_int, # indptr
_array_1d_int, # indices
_array_1d_dbl, # data
_array_1d_dbl, # n0
c_double, # dt
c_int, # order
c_int, # substeps
_array_1d_dbl, # result
]
def cram_solve(A, n0, dt, order=48, substeps=1):
"""Solve a single Bateman system using C++ CRAM.
Parameters
----------
A : scipy.sparse.csc_array or scipy.sparse.csc_matrix
Sparse transmutation matrix in CSC format.
n0 : numpy.ndarray
Initial atom number vector.
dt : float
Time step in seconds.
order : int
CRAM approximation order (16 or 48).
substeps : int
Number of equal substeps to use within ``dt``.
Returns
-------
numpy.ndarray
Final atom numbers.
"""
if order not in (16, 48):
raise ValueError(f"CRAM order must be 16 or 48, got {order}")
if not isinstance(substeps, numbers.Integral):
raise TypeError(f"substeps must be an integer, got {type(substeps)}")
if substeps <= 0:
raise ValueError(f"substeps must be positive, got {substeps}")
n = A.shape[0]
indptr = np.asarray(A.indptr, dtype=np.int32)
indices = np.asarray(A.indices, dtype=np.int32)
data = np.asarray(A.data, dtype=np.float64)
n0 = np.asarray(n0, dtype=np.float64)
result = np.empty(n, dtype=np.float64)
_dll.openmc_cram_solve(
n, indptr, indices, data, n0, dt, order, int(substeps), result)
return result

377
src/bateman_solvers.cpp Normal file
View file

@ -0,0 +1,377 @@
//! \file bateman_solvers.cpp
//! \brief Implementation of Bateman equation solvers
#include "openmc/bateman_solvers.h"
#include <complex>
#include "openmc/capi.h"
#include "openmc/error.h"
namespace openmc {
namespace {
// Fast complex multiply: (a+bi)(c+di) = (ac-bd) + (ad+bc)i
// Avoids GCC's __muldc3 which includes NaN recovery we don't need.
inline std::complex<double> fast_cmul(
std::complex<double> x, std::complex<double> y)
{
double a = x.real(), b = x.imag();
double c = y.real(), d = y.imag();
return {a * c - b * d, a * d + b * c};
}
// Fast complex reciprocal: 1/(a+bi) = (a-bi) / (a^2 + b^2)
inline std::complex<double> fast_crecip(std::complex<double> z)
{
double a = z.real();
double b = z.imag();
double denom = a * a + b * b;
return {a / denom, -b / denom};
}
} // anonymous namespace
//==============================================================================
// LU factorization for the CRAM shifted operator
//==============================================================================
void numeric_factorize_cram(const CSCMatrix& A, double dt,
std::complex<double> theta, const SymbolicLUFactorization& symbolic,
NumericLUFactorization& numeric, vector<std::complex<double>>& work)
{
int n = A.n();
const auto& a_indptr = A.indptr();
const auto& a_indices = A.indices();
const auto& a_data = A.data();
const auto& sp_indptr = symbolic.pattern.indptr();
const auto& sp_indices = symbolic.pattern.indices();
const auto& l_indptr = symbolic.l_pattern.indptr();
const auto& l_rowidx = symbolic.l_pattern.indices();
const auto& u_indptr = symbolic.u_pattern.indptr();
const auto& u_rowidx = symbolic.u_pattern.indices();
if (static_cast<int>(work.size()) != n) {
work.resize(n);
}
numeric.l_data.resize(symbolic.l_pattern.nnz());
numeric.u_data.resize(symbolic.u_pattern.nnz());
numeric.u_diag_inv.resize(n);
for (int j = 0; j < n; ++j) {
// Scatter the shifted operator column: M[:, j] = dt * A[:, j] - theta *
// I[:, j]
{
int a_pos = a_indptr[j];
int a_end = a_indptr[j + 1];
for (int sp_pos = sp_indptr[j]; sp_pos < sp_indptr[j + 1]; ++sp_pos) {
int row = sp_indices[sp_pos];
std::complex<double> val(0.0, 0.0);
if (a_pos < a_end && a_indices[a_pos] == row) {
val = dt * a_data[a_pos];
++a_pos;
}
if (row == j) {
val -= theta;
}
work[row] = val;
}
}
// Left-looking updates.
for (int up = u_indptr[j]; up < u_indptr[j + 1] - 1; ++up) {
int k = u_rowidx[up];
std::complex<double> ukj = work[k];
numeric.u_data[up] = ukj;
for (int lp = l_indptr[k]; lp < l_indptr[k + 1]; ++lp) {
work[l_rowidx[lp]] -= fast_cmul(numeric.l_data[lp], ukj);
}
}
std::complex<double> inv_ujj = fast_crecip(work[j]);
numeric.u_diag_inv[j] = inv_ujj;
numeric.u_data[u_indptr[j + 1] - 1] = work[j];
for (int lp = l_indptr[j]; lp < l_indptr[j + 1]; ++lp) {
numeric.l_data[lp] = fast_cmul(work[l_rowidx[lp]], inv_ujj);
}
// Reset the scatter buffer to zero for the next column. The union of
// U[:, j] and L[:, j] row indices is a superset of every row written by
// the scatter and left-looking updates above, so these two loops suffice.
for (int up = u_indptr[j]; up < u_indptr[j + 1]; ++up) {
work[u_rowidx[up]] = {0.0, 0.0};
}
for (int lp = l_indptr[j]; lp < l_indptr[j + 1]; ++lp) {
work[l_rowidx[lp]] = {0.0, 0.0};
}
// Invariant: the union of U[:, j] rows and L[:, j] rows is a superset of
// every row touched by the scatter and left-looking updates above, so
// these two loops fully reset `work` to zero before the next column.
}
}
void triangular_solve_lu(const vector<double>& b,
const SymbolicLUFactorization& symbolic,
const NumericLUFactorization& numeric, vector<std::complex<double>>& x)
{
int n = symbolic.pattern.n();
const auto& l_indptr = symbolic.l_pattern.indptr();
const auto& l_rowidx = symbolic.l_pattern.indices();
const auto& u_indptr = symbolic.u_pattern.indptr();
const auto& u_rowidx = symbolic.u_pattern.indices();
if (static_cast<int>(x.size()) != n) {
x.resize(n);
}
for (int j = 0; j < n; ++j) {
x[j] = std::complex<double>(b[j], 0.0);
}
for (int j = 0; j < n; ++j) {
for (int lp = l_indptr[j]; lp < l_indptr[j + 1]; ++lp) {
x[l_rowidx[lp]] -= fast_cmul(numeric.l_data[lp], x[j]);
}
}
for (int j = n - 1; j >= 0; --j) {
x[j] = fast_cmul(x[j], numeric.u_diag_inv[j]);
for (int up = u_indptr[j]; up < u_indptr[j + 1] - 1; ++up) {
x[u_rowidx[up]] -= fast_cmul(numeric.u_data[up], x[j]);
}
}
}
//==============================================================================
// CRAM coefficient tables
//
// Coefficients from M. Pusa, "Higher-Order Chebyshev Rational Approximation
// Method and Application to Burnup Equations," Nucl. Sci. Eng., 182:3,
// 297-318 (2016).
//==============================================================================
namespace {
// --- CRAM16 coefficients (8 poles) ---
constexpr std::complex<double> cram16_theta[] = {
{+3.509103608414918e+0, +8.436198985884374e+0},
{+5.948152268951177e+0, +3.587457362018322e+0},
{-5.264971343442647e+0, +1.622022147316793e+1},
{+1.419375897185666e+0, +1.092536348449672e+1},
{+6.416177699099435e+0, +1.194122393370139e+0},
{+4.993174737717997e+0, +5.996881713603942e+0},
{-1.413928462488886e+0, +1.349772569889275e+1},
{-1.084391707869699e+1, +1.927744616718165e+1},
};
constexpr std::complex<double> cram16_alpha[] = {
{+5.464930576870210e+3, -3.797983575308356e+4},
{+9.045112476907548e+1, -1.115537522430261e+3},
{+2.344818070467641e+2, -4.228020157070496e+2},
{+9.453304067358312e+1, -2.951294291446048e+2},
{+7.283792954673409e+2, -1.205646080220011e+5},
{+3.648229059594851e+1, -1.155509621409682e+2},
{+2.547321630156819e+1, -2.639500283021502e+1},
{+2.394538338734709e+1, -5.650522971778156e+0},
};
constexpr double cram16_alpha0 = 2.124853710495224e-16;
// --- CRAM48 coefficients (24 poles) ---
constexpr std::complex<double> cram48_theta[] = {
{-4.465731934165702e+1, +6.233225190695437e+1},
{-5.284616241568964e+0, +4.057499381311059e+1},
{-8.867715667624458e+0, +4.325515754166724e+1},
{+3.493013124279215e+0, +3.281615453173585e+1},
{+1.564102508858634e+1, +1.558061616372237e+1},
{+1.742097597385893e+1, +1.076629305714420e+1},
{-2.834466755180654e+1, +5.492841024648724e+1},
{+1.661569367939544e+1, +1.316994930024688e+1},
{+8.011836167974721e+0, +2.780232111309410e+1},
{-2.056267541998229e+0, +3.794824788914354e+1},
{+1.449208170441839e+1, +1.799988210051809e+1},
{+1.853807176907916e+1, +5.974332563100539e+0},
{+9.932562704505182e+0, +2.532823409972962e+1},
{-2.244223871767187e+1, +5.179633600312162e+1},
{+8.590014121680897e-1, +3.536456194294350e+1},
{-1.286192925744479e+1, +4.600304902833652e+1},
{+1.164596909542055e+1, +2.287153304140217e+1},
{+1.806076684783089e+1, +8.368200580099821e+0},
{+5.870672154659249e+0, +3.029700159040121e+1},
{-3.542938819659747e+1, +5.834381701800013e+1},
{+1.901323489060250e+1, +1.194282058271408e+0},
{+1.885508331552577e+1, +3.583428564427879e+0},
{-1.734689708174982e+1, +4.883941101108207e+1},
{+1.316284237125190e+1, +2.042951874827759e+1},
};
constexpr std::complex<double> cram48_alpha[] = {
{+6.387380733878774e+2, -6.743912502859256e+2},
{+1.909896179065730e+2, -3.973203432721332e+2},
{+4.236195226571914e+2, -2.041233768918671e+3},
{+4.645770595258726e+2, -1.652917287299683e+3},
{+7.765163276752433e+2, -1.783617639907328e+4},
{+1.907115136768522e+3, -5.887068595142284e+4},
{+2.909892685603256e+3, -9.953255345514560e+3},
{+1.944772206620450e+2, -1.427131226068449e+3},
{+1.382799786972332e+5, -3.256885197214938e+6},
{+5.628442079602433e+3, -2.924284515884309e+4},
{+2.151681283794220e+2, -1.121774011188224e+3},
{+1.324720240514420e+3, -6.370088443140973e+4},
{+1.617548476343347e+4, -1.008798413156542e+6},
{+1.112729040439685e+2, -8.837109731680418e+1},
{+1.074624783191125e+2, -1.457246116408180e+2},
{+8.835727765158191e+1, -6.388286188419360e+1},
{+9.354078136054179e+1, -2.195424319460237e+2},
{+9.418142823531573e+1, -6.719055740098035e+2},
{+1.040012390717851e+2, -1.693747595553868e+2},
{+6.861882624343235e+1, -1.177598523430493e+1},
{+8.766654491283722e+1, -4.596464999363902e+3},
{+1.056007619389650e+2, -1.738294585524067e+3},
{+7.738987569039419e+1, -4.311715386228984e+1},
{+1.041366366475571e+2, -2.777743732451969e+2},
};
constexpr double cram48_alpha0 = 2.258038182743983e-47;
} // anonymous namespace
//==============================================================================
// IPFCramSolver implementation
//==============================================================================
IPFCramSolver::IPFCramSolver(int order)
{
if (order == 16) {
n_poles_ = 8;
alpha_ = cram16_alpha;
theta_ = cram16_theta;
alpha0_ = cram16_alpha0;
} else if (order == 48) {
n_poles_ = 24;
alpha_ = cram48_alpha;
theta_ = cram48_theta;
alpha0_ = cram48_alpha0;
} else {
throw std::invalid_argument {
fmt::format("CRAM order must be 16 or 48, got {}.", order)};
}
}
vector<double> IPFCramSolver::solve(
const CSCMatrix& A, const vector<double>& n0, double dt, int substeps)
{
if (substeps <= 0) {
throw std::invalid_argument {
fmt::format("substeps must be positive, got {}.", substeps)};
}
int n = A.n();
double step_dt = dt / substeps;
// Symbolic factorization: compute L/U sparsity patterns for this matrix.
// Reused across all pole solves below.
auto symbolic = symbolic_factorize(A.pattern().with_diagonal());
vector<std::complex<double>> work(n);
vector<std::complex<double>> x(n);
// IPF CRAM iteration:
// y_0 = n0
// y_{k+1} = y_k + 2*Re(alpha_k * (A*dt - theta_k*I)^{-1} * y_k)
// result = alpha0 * y_final
vector<double> y(n0.begin(), n0.end());
if (substeps > 1) {
vector<NumericLUFactorization> cached_pole_factorizations(n_poles_);
for (int p = 0; p < n_poles_; ++p) {
numeric_factorize_cram(
A, step_dt, theta_[p], symbolic, cached_pole_factorizations[p], work);
}
for (int step = 0; step < substeps; ++step) {
for (int p = 0; p < n_poles_; ++p) {
triangular_solve_lu(y, symbolic, cached_pole_factorizations[p], x);
// y += 2 * Re(alpha_p * x_p)
for (int i = 0; i < n; ++i) {
auto ax = fast_cmul(alpha_[p], x[i]);
y[i] += 2.0 * ax.real();
}
}
for (int i = 0; i < n; ++i) {
y[i] *= alpha0_;
}
}
} else {
// Reuse a single numeric container while recomputing each pole on demand.
NumericLUFactorization current_pole_factorization;
for (int p = 0; p < n_poles_; ++p) {
numeric_factorize_cram(
A, step_dt, theta_[p], symbolic, current_pole_factorization, work);
triangular_solve_lu(y, symbolic, current_pole_factorization, x);
// y += 2 * Re(alpha_p * x_p)
for (int i = 0; i < n; ++i) {
auto ax = fast_cmul(alpha_[p], x[i]);
y[i] += 2.0 * ax.real();
}
}
for (int i = 0; i < n; ++i) {
y[i] *= alpha0_;
}
}
return y;
}
} // namespace openmc
//==============================================================================
// C API
//==============================================================================
using namespace openmc;
extern "C" int openmc_cram_solve(int n, const int* indptr, const int* indices,
const double* data, const double* n0, double dt, int order, int substeps,
double* result)
{
if (!indptr || !indices || !data || !n0 || !result) {
set_errmsg("openmc_cram_solve: null pointer argument");
return OPENMC_E_INVALID_ARGUMENT;
}
if (n < 0) {
set_errmsg(fmt::format("matrix dimension must be non-negative, got {}", n));
return OPENMC_E_INVALID_ARGUMENT;
}
try {
int nnz = indptr[n];
CSCMatrix A(n, vector<int>(indptr, indptr + n + 1),
vector<int>(indices, indices + nnz), vector<double>(data, data + nnz));
vector<double> n0_vec(n0, n0 + n);
IPFCramSolver solver(order);
vector<double> y = solver.solve(A, n0_vec, dt, substeps);
std::copy(y.begin(), y.end(), result);
} catch (const std::invalid_argument& e) {
set_errmsg(e.what());
return OPENMC_E_INVALID_ARGUMENT;
} catch (const std::exception& e) {
set_errmsg(e.what());
return OPENMC_E_DATA;
}
return 0;
}

193
src/sparse_matrix.cpp Normal file
View file

@ -0,0 +1,193 @@
//! \file sparse_matrix.cpp
//! \brief Implementation of CSCPattern and CSCMatrix
#include "openmc/sparse_matrix.h"
#include <algorithm>
#include <stdexcept>
#include <fmt/core.h>
namespace openmc {
//==============================================================================
// CSCPattern implementation
//==============================================================================
CSCPattern::CSCPattern(int n, vector<int> indptr, vector<int> indices)
: n_(n), indptr_(std::move(indptr)), indices_(std::move(indices))
{
if (static_cast<int>(indptr_.size()) != n_ + 1) {
throw std::invalid_argument {fmt::format(
"CSCPattern: indptr size ({}) != n + 1 ({})", indptr_.size(), n_ + 1)};
}
if (indptr_[0] != 0) {
throw std::invalid_argument {
fmt::format("CSCPattern: indptr[0] ({}) != 0", indptr_[0])};
}
if (indptr_[n_] != static_cast<int>(indices_.size())) {
throw std::invalid_argument {
fmt::format("CSCPattern: indptr[n] ({}) != indices size ({})",
indptr_[n_], indices_.size())};
}
for (int j = 0; j < n_; ++j) {
for (int k = indptr_[j]; k < indptr_[j + 1]; ++k) {
if (indices_[k] < 0 || indices_[k] >= n_) {
throw std::invalid_argument {fmt::format(
"CSCPattern: row index {} out of bounds [0, {}) in column {}",
indices_[k], n_, j)};
}
if (k > indptr_[j] && indices_[k - 1] >= indices_[k]) {
throw std::invalid_argument {
fmt::format("CSCPattern: row indices not sorted in column {} "
"(indices[{}]={} >= indices[{}]={})",
j, k - 1, indices_[k - 1], k, indices_[k])};
}
}
}
}
//==============================================================================
// CSCMatrix implementation
//==============================================================================
CSCMatrix::CSCMatrix(
int n, vector<int> indptr, vector<int> indices, vector<double> data)
: pattern_(n, std::move(indptr), std::move(indices)), data_(std::move(data))
{
if (static_cast<int>(data_.size()) != pattern_.nnz()) {
throw std::invalid_argument {
fmt::format("CSCMatrix: data size ({}) != pattern nnz ({})", data_.size(),
pattern_.nnz())};
}
}
bool CSCPattern::operator==(const CSCPattern& other) const
{
return n_ == other.n_ && indptr_ == other.indptr_ &&
indices_ == other.indices_;
}
CSCPattern CSCPattern::with_diagonal() const
{
// Single-pass merge: for each column, copy the existing sorted row indices
// while inserting `col` in its sorted position if absent. The final nnz is
// at most nnz() + n_, so we reserve that upper bound.
vector<int> new_indptr(n_ + 1);
vector<int> new_indices;
new_indices.reserve(indices_.size() + n_);
for (int col = 0; col < n_; ++col) {
new_indptr[col] = static_cast<int>(new_indices.size());
bool diag_written = false;
for (int idx = indptr_[col]; idx < indptr_[col + 1]; ++idx) {
int row = indices_[idx];
if (!diag_written && row >= col) {
new_indices.push_back(col);
diag_written = true;
if (row == col)
continue; // don't duplicate an existing diagonal
}
new_indices.push_back(row);
}
if (!diag_written) {
new_indices.push_back(col);
}
}
new_indptr[n_] = static_cast<int>(new_indices.size());
return CSCPattern(n_, std::move(new_indptr), std::move(new_indices));
}
//==============================================================================
// LU factorization helpers
//==============================================================================
SymbolicLUFactorization symbolic_factorize(CSCPattern pattern)
{
int n = pattern.n();
const auto& indptr = pattern.indptr();
const auto& indices = pattern.indices();
// Build L and U fill patterns column by column.
vector<vector<int>> l_cols(n);
vector<bool> marked(n, false);
vector<int> u_work, l_work;
// Temporary storage for U column patterns before CSC assembly.
vector<vector<int>> u_cols(n);
for (int j = 0; j < n; ++j) {
u_work.clear();
l_work.clear();
// Scatter: mark off-diagonal structural entries in column j.
for (int p = indptr[j]; p < indptr[j + 1]; ++p) {
int i = indices[p];
if (i == j)
continue;
if (!marked[i]) {
marked[i] = true;
if (i < j) {
u_work.push_back(i);
} else {
l_work.push_back(i);
}
}
}
// Propagate fill through previously discovered L columns.
for (size_t idx = 0; idx < u_work.size(); ++idx) {
int k = u_work[idx];
for (int row : l_cols[k]) {
if (row == j)
continue;
if (!marked[row]) {
marked[row] = true;
if (row < j) {
u_work.push_back(row);
} else {
l_work.push_back(row);
}
}
}
}
std::sort(u_work.begin(), u_work.end());
std::sort(l_work.begin(), l_work.end());
u_cols[j] = u_work;
l_cols[j] = l_work;
for (int k : u_work)
marked[k] = false;
for (int i : l_work)
marked[i] = false;
}
vector<int> l_indptr(n + 1);
vector<int> l_rowidx;
for (int j = 0; j < n; ++j) {
l_indptr[j] = static_cast<int>(l_rowidx.size());
for (int r : l_cols[j])
l_rowidx.push_back(r);
}
l_indptr[n] = static_cast<int>(l_rowidx.size());
// Store the diagonal as the last U entry in each column.
vector<int> u_indptr(n + 1);
vector<int> u_rowidx;
for (int j = 0; j < n; ++j) {
u_indptr[j] = static_cast<int>(u_rowidx.size());
for (int r : u_cols[j])
u_rowidx.push_back(r);
u_rowidx.push_back(j);
}
u_indptr[n] = static_cast<int>(u_rowidx.size());
return {std::move(pattern),
CSCPattern(n, std::move(l_indptr), std::move(l_rowidx)),
CSCPattern(n, std::move(u_indptr), std::move(u_rowidx))};
}
} // namespace openmc

1
tests/cpp_unit_tests/CMakeLists.txt
View file

@ -1,6 +1,7 @@
set(TEST_NAMES
test_distribution
test_file_utils
test_sparse_matrix
test_tally
test_interpolate
test_math

166
tests/cpp_unit_tests/test_sparse_matrix.cpp Normal file
View file

@ -0,0 +1,166 @@
#include <stdexcept>
#include <catch2/catch_test_macros.hpp>
#include "openmc/sparse_matrix.h"
using namespace openmc;
TEST_CASE("CSCPattern construction and accessors")
{
// 3x3 pattern:
// col 0: rows 0, 2
// col 1: row 1
// col 2: rows 0, 2
vector<int> indptr = {0, 2, 3, 5};
vector<int> indices = {0, 2, 1, 0, 2};
CSCPattern pat(3, indptr, indices);
CHECK(pat.n() == 3);
CHECK(pat.nnz() == 5);
CHECK(pat.indptr() == indptr);
CHECK(pat.indices() == indices);
}
TEST_CASE("CSCPattern default constructor")
{
CSCPattern pat;
CHECK(pat.n() == 0);
CHECK(pat.nnz() == 0);
}
TEST_CASE("CSCPattern operator==")
{
vector<int> indptr = {0, 2, 3, 5};
vector<int> indices = {0, 2, 1, 0, 2};
CSCPattern a(3, indptr, indices);
CSCPattern b(3, indptr, indices);
CHECK(a == b);
CHECK_FALSE(a != b);
// Different dimension
CSCPattern c(4, {0, 0, 0, 0, 0}, {});
CHECK(a != c);
// Same dimension but different pattern
CSCPattern d(3, {0, 1, 2, 3}, {0, 1, 2});
CHECK(a != d);
}
TEST_CASE("CSCPattern with_diagonal inserts missing diagonals")
{
// Pattern with no diagonal entries:
// col 0: row 2
// col 1: (empty)
// col 2: row 0
CSCPattern pat(3, {0, 1, 1, 2}, {2, 0});
auto wd = pat.with_diagonal();
CHECK(wd.n() == 3);
// Original 2 entries + 3 missing diagonals = 5
CHECK(wd.nnz() == 5);
// Verify diagonals are present and rows stay sorted
const auto& ip = wd.indptr();
const auto& ix = wd.indices();
// col 0: rows 0(new diag), 2
CHECK(ip[0] == 0);
CHECK(ix[0] == 0);
CHECK(ix[1] == 2);
// col 1: row 1(new diag)
CHECK(ip[1] == 2);
CHECK(ix[2] == 1);
// col 2: rows 0, 2(new diag)
CHECK(ip[2] == 3);
CHECK(ix[3] == 0);
CHECK(ix[4] == 2);
CHECK(ip[3] == 5);
}
TEST_CASE("CSCPattern with_diagonal preserves existing diagonals")
{
// Full diagonal already present
CSCPattern pat(3, {0, 2, 4, 6}, {0, 1, 1, 2, 0, 2});
auto wd = pat.with_diagonal();
CHECK(wd.n() == 3);
CHECK(wd.nnz() == 6); // unchanged
CHECK(wd == pat);
}
TEST_CASE("CSCPattern with_diagonal partial diagonals")
{
// 3x3 pattern: col 0 has diag, col 1 missing, col 2 has diag
// col 0: rows 0, 2
// col 1: row 0
// col 2: rows 1, 2
CSCPattern pat(3, {0, 2, 3, 5}, {0, 2, 0, 1, 2});
auto wd = pat.with_diagonal();
CHECK(wd.nnz() == 6); // 5 + 1 missing diagonal for col 1
// col 1 should now have rows 0, 1 (diagonal inserted in sorted order)
const auto& ip = wd.indptr();
const auto& ix = wd.indices();
CHECK(ix[ip[1]] == 0);
CHECK(ix[ip[1] + 1] == 1);
}
TEST_CASE("CSCMatrix construction and accessors")
{
// 3x3 matrix:
// [1 0 2]
// [0 3 0]
// [4 0 5]
vector<int> indptr = {0, 2, 3, 5};
vector<int> indices = {0, 2, 1, 0, 2};
vector<double> data = {1.0, 4.0, 3.0, 2.0, 5.0};
CSCMatrix mat(3, indptr, indices, data);
CHECK(mat.n() == 3);
CHECK(mat.nnz() == 5);
CHECK(mat.indptr() == indptr);
CHECK(mat.indices() == indices);
CHECK(mat.data() == data);
// Verify pattern matches equivalent standalone CSCPattern
CSCPattern pat(3, indptr, indices);
CHECK(mat.pattern() == pat);
}
TEST_CASE("CSCMatrix default constructor")
{
CSCMatrix mat;
CHECK(mat.n() == 0);
CHECK(mat.nnz() == 0);
}
TEST_CASE("CSCPattern rejects malformed inputs")
{
// indptr size mismatch (n=3 requires 4 entries)
CHECK_THROWS_AS(CSCPattern(3, {0, 1, 2}, {0, 1}), std::invalid_argument);
// indptr[0] must be 0
CHECK_THROWS_AS(CSCPattern(2, {1, 1, 2}, {0, 1}), std::invalid_argument);
// indptr[n] must equal nnz
CHECK_THROWS_AS(CSCPattern(2, {0, 1, 3}, {0, 1}), std::invalid_argument);
// Row index out of bounds
CHECK_THROWS_AS(CSCPattern(2, {0, 1, 2}, {0, 5}), std::invalid_argument);
// Unsorted row indices within a column
CHECK_THROWS_AS(CSCPattern(3, {0, 2, 2, 2}, {2, 0}), std::invalid_argument);
// Duplicate row indices within a column (not strictly ascending)
CHECK_THROWS_AS(CSCPattern(3, {0, 2, 2, 2}, {1, 1}), std::invalid_argument);
}
TEST_CASE("CSCMatrix rejects data/nnz size mismatch")
{
// Pattern has 2 nonzeros, only 1 data value supplied
CHECK_THROWS_AS(
CSCMatrix(2, {0, 1, 2}, {0, 1}, {1.0}), std::invalid_argument);
}

77
tests/unit_tests/test_deplete_cram.py
View file

@ -8,93 +8,118 @@ import numpy as np
import pytest
import scipy.sparse as sp
from pytest import approx
from openmc.deplete.cram import (CRAM16, CRAM48, Cram16Solver, Cram48Solver,
IPFCramSolver)
from openmc.deplete.cram import CRAM16, CRAM48
from openmc.lib.deplete import cram_solve
def test_CRAM16():
"""Test 16-term CRAM."""
"""Test 16th order CRAM against Mathematica reference."""
x = np.array([1.0, 1.0])
mat = sp.csr_matrix([[-1.0, 0.0], [-2.0, -3.0]])
mat = sp.csc_array([[-1.0, 0.0], [-2.0, -3.0]])
dt = 0.1
z = CRAM16(mat, x, dt)
# Solution from mathematica
z0 = np.array((0.904837418035960, 0.576799023327476))
assert z == approx(z0)
def test_CRAM48():
"""Test 48-term CRAM."""
"""Test 48th order CRAM against Mathematica reference."""
x = np.array([1.0, 1.0])
mat = sp.csr_matrix([[-1.0, 0.0], [-2.0, -3.0]])
mat = sp.csc_array([[-1.0, 0.0], [-2.0, -3.0]])
dt = 0.1
z = CRAM48(mat, x, dt)
# Solution from mathematica
z0 = np.array((0.904837418035960, 0.576799023327476))
assert z == approx(z0)
def test_cram_solve_cpp():
"""Test C++ cram_solve binding directly against Mathematica reference."""
A = sp.csc_array([[-1.0, 0.0], [-2.0, -3.0]])
n0 = np.array([1.0, 1.0])
dt = 0.1
z48 = cram_solve(A, n0, dt, order=48)
z16 = cram_solve(A, n0, dt, order=16)
z_ref = np.array((0.904837418035960, 0.576799023327476))
assert z48 == approx(z_ref)
assert z16 == approx(z_ref)
def test_substeps1_matches_original():
"""substeps=1 must be bitwise identical to original spsolve path."""
"""substeps=1 must be bitwise identical to original result."""
x = np.array([1.0, 1.0])
mat = sp.csr_matrix([[-1.0, 0.0], [-2.0, -3.0]])
mat = sp.csc_array([[-1.0, 0.0], [-2.0, -3.0]])
dt = 0.1
z_orig = CRAM48(mat, x, dt)
z_sub1 = CRAM48(mat, x, dt, substeps=1)
z_cpp = cram_solve(mat, x, dt, order=48, substeps=1)
np.testing.assert_array_equal(z_sub1, z_orig)
np.testing.assert_array_equal(z_cpp, z_orig)
def test_substeps2_matches_two_half_steps():
"""substeps=2 must match two independent CRAM calls with dt/2."""
x = np.array([1.0, 1.0])
mat = sp.csr_matrix([[-1.0, 0.0], [-2.0, -3.0]])
mat = sp.csc_array([[-1.0, 0.0], [-2.0, -3.0]])
dt = 1.0
# Two manual half-steps using original spsolve path
z_half = CRAM48(mat, x, dt / 2)
z_two = CRAM48(mat, z_half, dt / 2)
# Single call with substeps=2
z_sub2 = CRAM48(mat, x, dt, substeps=2)
z_cpp = cram_solve(mat, x, dt, order=48, substeps=2)
assert z_sub2 == approx(z_two, rel=1e-12)
assert z_cpp == approx(z_two, rel=1e-12)
@pytest.mark.parametrize("substeps", [0, -1])
def test_invalid_substeps(substeps):
"""substeps must be a positive integer at call time."""
x = np.array([1.0, 1.0])
mat = sp.csr_matrix([[-1.0, 0.0], [-2.0, -3.0]])
mat = sp.csc_array([[-1.0, 0.0], [-2.0, -3.0]])
with pytest.raises(ValueError, match="substeps"):
CRAM48(mat, x, 0.1, substeps=substeps)
with pytest.raises(ValueError, match="substeps"):
cram_solve(mat, x, 0.1, order=48, substeps=substeps)
@pytest.mark.parametrize("order", [0, 1, 32, 64])
def test_invalid_order(order):
"""cram_solve must reject unsupported orders."""
mat = sp.csc_array([[-1.0, 0.0], [-2.0, -3.0]])
x = np.array([1.0, 1.0])
with pytest.raises(ValueError, match="order"):
cram_solve(mat, x, 0.1, order=order)
def test_substeps_self_convergence():
"""Increasing substeps converges toward reference solution.
Uses CRAM16 (alpha0 ~ 2e-16) where substep convergence is visible.
CRAM48 (alpha0 ~ 2e-47) is already near machine precision for small
systems; its correctness is verified by the other substep tests.
Uses CRAM16 where substep convergence is visible.
"""
mat = sp.csr_matrix([[-1.0, 0.0], [-2.0, -3.0]])
mat = sp.csc_array([[-1.0, 0.0], [-2.0, -3.0]])
x = np.array([1.0, 1.0])
dt = 50 # lambda*dt = 50 and 150, stresses CRAM16
dt = 50.0
n_ref = CRAM16(mat, x, dt, substeps=128)
prev_err = np.inf
for s in [1, 2, 4, 8, 16]:
n_s = CRAM16(mat, x, dt, substeps=s)
err = np.linalg.norm(n_s - n_ref) / np.linalg.norm(n_ref)
assert err < prev_err, \
f"substeps={s} error {err:.2e} not less than previous {prev_err:.2e}"
prev_err = err
for substeps in [1, 2, 4, 8, 16]:
n_sub = CRAM16(mat, x, dt, substeps=substeps)
err = np.linalg.norm(n_sub - n_ref) / np.linalg.norm(n_ref)
assert err < prev_err
prev_err = err

4
tests/unit_tests/test_deplete_integrator.py
View file

@ -250,9 +250,9 @@ def test_integrator(run_in_tmpdir, scheme):
integrator = bundle.solver(operator, [0.75], 1, solver="cram16")
assert integrator.solver is cram.CRAM16
integrator = bundle.solver(operator, [0.75], 1, solver=cram.Cram48Solver,
integrator = bundle.solver(operator, [0.75], 1, solver=cram.CRAM48,
substeps=2)
assert integrator.solver is cram.Cram48Solver
assert integrator.solver is cram.CRAM48
assert integrator.substeps == 2
integrator.solver = mock_good_solver