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

Merge 9fbf7aca06 into 61c8a59cff

Browse source
This commit is contained in:
Lorenzo Chierici 2026-07-21 07:53:51 +00:00 • committed by GitHub
commit b4eff7cbfb
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
6 changed files with 213 additions and 32 deletions
Download patch file Download diff file Expand all files Collapse all files

54
docs/source/methods/depletion.rst
View file

@ -339,3 +339,57 @@ where:
Note that mass conservation is guaranteed by transferring the number
of atoms directly.
--------------
External Source Rates
--------------
OpenMC allows the addition of external source rates to the depletion matrix.
This is useful for modeling the feed or removal of fixed amounts of nuclides
to or from a depletable material. Rates are specified as a mass flow (default
units of g/s) and converted to atom/s source terms for the nuclides in the
composition vector. A positive rate corresponds to feed and a negative rate to
removal.
Mathematically, this adds an external source term :math:`S_i` to the depletion
equation:
.. math::
\frac{dN_i}{dt} = \sum_j A_{ij} N_j + S_i
The resulting linear system is non-homogeneous but can be recast in homogeneous
form by augmenting the nuclide vector with a constant component equal to unity:
.. math::
\frac{d}{dt}\begin{bmatrix}\mathbf{N}\\ 1\end{bmatrix} =
\begin{bmatrix}
\mathbf{A} & \mathbf{S}\\
\mathbf{0} & 0
\end{bmatrix}
\begin{bmatrix}
\mathbf{N}\\
1
\end{bmatrix}
where :math:`\mathbf{S}` is the vector of external source rates in atom/s.
The resulting system can be solved with the same integration algorithms that are
used in the absence of the external source term.
External source rates with transfer rates coupling materials
------------------------------------------------------------
In the presence of external source rates and coupled transfer rates between
materials, the off-diagonal transfer blocks must be augmented to the same
dimensions as the corresponding diagonal blocks. Using the transfer example
above, if an external source rate is applied to material 1 (the losing
material), the coupling matrix :math:`\mathbf{T_{21}}` is extended with an
additional column of zeroes so that its column count matches the augmented
:math:`\mathbf{A_{11}}`. If the external source is applied to material 2 (the
receiving material), :math:`\mathbf{T_{21}}` instead receives an additional row
of zeroes so that its row count matches the augmented :math:`\mathbf{A_{22}}`.

67
docs/source/usersguide/depletion.rst
View file

@ -450,8 +450,73 @@ to transfer xenon from one material to another, you'd use::
integrator.add_transfer_rate(mat1, ['Xe'], 0.1, destination_material=mat2)
External Source Rates
External source rates define a fixed mass feed or removal of nuclides to or
from a depletable material. Unlike transfer rates, which are proportional to
the instantaneous nuclide inventory, external source rates add a constant
source term to the depletion equations. This can be useful to model batch
refueling, makeup fuel addition, or fixed-rate off-gas removal.
External source rates are defined by calling the
:meth:`~openmc.deplete.abc.Integrator.add_external_source_rate()` method
directly from one of the Integrator classes::
...
integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power)
integrator.add_external_source_rate(...)
Defining external source rates
------------------------------
The :meth:`~openmc.deplete.abc.Integrator.add_external_source_rate()` method
requires a :class:`~openmc.Material` instance (alternatively, a material id or
the name) as the depletable material, a composition dictionary giving the
relative weight fractions of elements and/or nuclides in the feed or removal
stream, and a mass flow rate.
.. caution::
Make sure you set the rate value with the right sign. A positive rate
corresponds to feed, while a negative rate corresponds to removal. This is
the opposite convention used for transfer rates.
The ``rate_units`` argument specifies the units for the mass flow rate. The
default is ``g/s``, but ``g/min``, ``g/h``, ``g/d``, and ``g/a`` are also valid
options.
For example, to feed uranium-235 into a material at 10 g/day, you'd use::
mat1 = openmc.Material(material_id=1, name='fuel')
...
integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power)
composition = {'U235': 1.0}
# by openmc.Material object
integrator.add_external_source_rate(mat1, composition, 10, rate_units='g/d')
# or by material id
integrator.add_external_source_rate(1, composition, 10, rate_units='g/d')
# or by material name
integrator.add_external_source_rate('fuel', composition, 10, rate_units='g/d')
Composition keys may be nuclides (e.g., ``'U235'``) or naturally abundant
elements (e.g., ``'U'``). When an element is specified, the mass flow is
distributed across its naturally occurring isotopes according to their natural
abundances.
The optional ``timesteps`` argument restricts the external source rate to
specific depletion step indices. If omitted, the rate is applied at every step.
Combining with transfer rates
-----------------------------
External source rates can be used together with transfer rates, including
transfers between materials via ``destination_material``. See
:ref:`methods_depletion` for the augmented-matrix formulation used when both
features are active.
Comparing to Other Codes
========================
Comparing depletion results from OpenMC with those from another code, such as
MCNP or Serpent, requires more than constructing equivalent transport models.

10
openmc/deplete/abc.py
View file

@ -1028,11 +1028,6 @@ class Integrator(ABC):
self.transfer_rates = TransferRates(
self.operator, materials, len(self.timesteps))
if self.external_source_rates is not None and destination_material:
raise ValueError('Currently is not possible to set a transfer rate '
'with destination matrial in combination with '
'external source rates.')
self.transfer_rates.set_transfer_rate(
material, components, transfer_rate, transfer_rate_units,
timesteps, destination_material)
@ -1074,11 +1069,6 @@ class Integrator(ABC):
self.external_source_rates = ExternalSourceRates(
self.operator, materials, len(self.timesteps))
if self.transfer_rates is not None and self.transfer_rates.index_transfer:
raise ValueError('Currently is not possible to set an external '
'source rate in combination with transfer rates '
'with destination matrial.')
self.external_source_rates.set_external_source_rate(
material, composition, rate, rate_units, timesteps)

79
openmc/deplete/pool.py
View file

@ -6,10 +6,10 @@ from itertools import repeat, starmap
from multiprocessing import Pool
import numpy as np
from scipy.sparse import hstack
from scipy.sparse import hstack, vstack
from openmc.mpi import comm
from .._sparse_compat import block_array
from .._sparse_compat import block_array, csc_array
# Configurable switch that enables / disables the use of
# multiprocessing routines during depletion
@ -62,7 +62,7 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
Time in [s] to deplete for
current_timestep : int
Current timestep index
maxtrix_func : callable, optional
matrix_func : callable, optional
Function to form the depletion matrix after calling ``matrix_func(chain,
rates, fission_yields)``, where ``fission_yields = {parent: {product:
yield_frac}}`` Expected to return the depletion matrix required by
@ -87,7 +87,6 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
list contains the number of [atom] of each nuclide.
"""
fission_yields = chain.fission_yields
if len(fission_yields) == 1:
fission_yields = repeat(fission_yields[0])
@ -103,6 +102,7 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
matrices = map(matrix_func, repeat(chain), rates, fission_yields,
*matrix_args)
n_solve = n
if (transfer_rates is not None and
current_timestep in transfer_rates.external_timesteps):
# Calculate transfer rate terms as diagonal matrices
@ -119,11 +119,32 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
matrices[mat_idx] = chain.add_redox_term(matrices[mat_idx],
transfer_rates.redox[mat_id][0],
transfer_rates.redox[mat_id][1])
# Add external sources if present
if (external_source_rates is not None and
current_timestep in external_source_rates.external_timesteps):
n_solve = [arr.copy() for arr in n]
sources = map(chain.form_ext_source_term, repeat(external_source_rates),
repeat(current_timestep), external_source_rates.local_mats)
# stack vector column at the end of the matrix
matrices = [
hstack([matrix, source])
for matrix, source in zip(matrices, sources)
]
# Homogenize diagonal matrices to be square
for i, matrix in enumerate(matrices):
if matrix.shape[0] + 1 == matrix.shape[1]:
# Add a row of zeroes to the matrix and append 1 to the last row of the nuclide vector
matrices[i] = vstack([matrix, csc_array((1, matrix.shape[1]))])
n_solve[i] = np.append(n_solve[i], 1.0)
# Set transfer rate terms with destination material if present
if current_timestep in transfer_rates.index_transfer:
# Gather all on comm.rank 0
matrices = comm.gather(matrices)
n = comm.gather(n)
n = comm.gather(n_solve)
if comm.rank == 0:
# Expand lists
@ -132,20 +153,27 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
# Calculate transfer rate terms as diagonal matrices
transfer_pair = {}
for mat_pair in transfer_rates.index_transfer[current_timestep]:
for mat_pair in dict.fromkeys(transfer_rates.index_transfer[current_timestep]):
transfer_matrix = chain.form_rr_term(transfer_rates,
current_timestep,
mat_pair)
# check if destination material has a redox control
if mat_pair[0] in transfer_rates.redox:
transfer_matrix = chain.add_redox_term(transfer_matrix,
transfer_rates.redox[mat_pair[0]][0],
transfer_rates.redox[mat_pair[0]][1])
# Add external source rates if present
if (external_source_rates is not None and
current_timestep in external_source_rates.external_timesteps):
if len(external_source_rates.get_components(mat_pair[0], current_timestep)) > 0:
transfer_matrix = vstack([transfer_matrix,
csc_array((1, transfer_matrix.shape[1]))])
if len(external_source_rates.get_components(mat_pair[1], current_timestep)) > 0:
transfer_matrix = hstack([transfer_matrix,
csc_array((transfer_matrix.shape[0], 1))])
transfer_pair[mat_pair] = transfer_matrix
# Combine all matrices together in a single matrix of matrices
# to be solved in one go
# Combine all matrices together in a single block matrix of matrices
# to be solved on one rank
n_rows = n_cols = len(transfer_rates.burnable_mats)
rows = []
for row in range(n_rows):
@ -172,18 +200,26 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
# Split back the nuclide vector result into the original form
n_result = np.split(n_result, np.cumsum([len(i) for i in n])[:-1])
# Remove extra value at the end of the nuclide vectors if external source rates are present
if (external_source_rates is not None and
current_timestep in external_source_rates.external_timesteps):
external_source_rates.reformat_nuclide_vectors(n_result)
#Clamp negative values to 0
for n_material in n_result:
np.maximum(n_material, 0.0, out=n_material)
else:
n_result = None
# Braodcast result to other ranks
n_result = comm.bcast(n_result)
# Distribute results across MPI
n_result = _distribute(n_result)
return n_result
if (external_source_rates is not None and
# If only external source rates are present
elif (external_source_rates is not None and
current_timestep in external_source_rates.external_timesteps):
n_solve = [arr.copy() for arr in n]
# Calculate external source term vectors
sources = map(chain.form_ext_source_term, repeat(external_source_rates),
repeat(current_timestep), external_source_rates.local_mats)
@ -193,15 +229,14 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
hstack([matrix, source])
for matrix, source in zip(matrices, sources)
]
# Add a last row of zeroes to the matrices and append 1 to the last row
# of the nuclide vectors
for i, matrix in enumerate(matrices):
if not np.equal(*matrix.shape):
matrix.resize(matrix.shape[1], matrix.shape[1])
n[i] = np.append(n[i], 1.0)
inputs = zip(matrices, n, repeat(dt), repeat(substeps))
if matrix.shape[0] + 1 == matrix.shape[1]:
matrices[i] = vstack([matrix, csc_array((1, matrix.shape[1]))])
n_solve[i] = np.append(n_solve[i], 1.0)
inputs = zip(matrices, n_solve, repeat(dt), repeat(substeps))
if USE_MULTIPROCESSING:
with Pool(NUM_PROCESSES) as pool:
@ -209,10 +244,12 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
else:
n_result = list(starmap(func, inputs))
# Remove extra value at the end of the nuclide vectors
# Remove extra value at the end of the nuclide vectors if external source rates are present
if (external_source_rates is not None and
current_timestep in external_source_rates.external_timesteps):
external_source_rates.reformat_nuclide_vectors(n)
external_source_rates.reformat_nuclide_vectors(n_result)
#Clamp negative values to 0
for n_material in n_result:
np.maximum(n_material, 0.0, out=n_material)
return n_result

BIN
tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_ext_source_and_transfer.h5 Normal file
View file

Binary file not shown.

35
tests/regression_tests/deplete_with_transfer_rates/test.py
View file

@ -126,3 +126,38 @@ def test_external_source_rates(run_in_tmpdir, model, rate, power, ref_result):
assert_atoms_equal(res_ref, res_test, tol=1e-3)
assert_reaction_rates_equal(res_ref, res_test, tol=1e-3)
@pytest.mark.parametrize("external_source_rate, transfer_rate, power, ref_result", [
(1e-1, 1e-1, 174., 'depletion_with_ext_source_and_transfer'),
])
def test_external_source_rates_with_transfer_rates(run_in_tmpdir, model, external_source_rate,
transfer_rate, power, ref_result):
"""Tests external_rates depletion class with external source rates and transfer rates"""
chain_file = Path(__file__).parents[2] / 'chain_simple.xml'
external_source_vector = {'U': 1}
op = CoupledOperator(model, chain_file)
op.round_number = True
integrator = openmc.deplete.PredictorIntegrator(
op, [1], power, timestep_units='d')
integrator.add_external_source_rate('f', external_source_vector, external_source_rate)
integrator.add_transfer_rate('f', ['U235'], transfer_rate, destination_material='w')
integrator.integrate()
# Get path to test and reference results
path_test = op.output_dir / 'depletion_results.h5'
path_reference = Path(__file__).with_name(f'ref_{ref_result}.h5')
# If updating results, do so and return
if config['update']:
shutil.copyfile(str(path_test), str(path_reference))
return
# Load the reference/test results
res_ref = openmc.deplete.Results(path_reference)
res_test = openmc.deplete.Results(path_test)
assert_atoms_equal(res_ref, res_test, tol=1e-3)
assert_reaction_rates_equal(res_ref, res_test, tol=1e-3)