From 52b49631faf60e71494c1718734f89108b535705 Mon Sep 17 00:00:00 2001 From: Perry Date: Wed, 11 Feb 2026 22:07:02 -0800 Subject: [PATCH 1/3] Clip negative atom densities after CRAM solves CRAM rational approximation can produce tiny negative atom densities as numerical artifacts. These propagate to depletion_results.h5 and to subsequent timesteps. Clip all CRAM outputs to non-negative in _timed_deplete(), the single funnel point all integrators pass through. --- openmc/deplete/abc.py | 2 ++ tests/unit_tests/test_deplete_cram.py | 32 +++++++++++++++++++++++++++ 2 files changed, 34 insertions(+) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 056f7c2737a..3c0e8e4c4d3 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -736,6 +736,8 @@ def _timed_deplete(self, n, rates, dt, i=None, matrix_func=None): results = deplete( self._solver, self.chain, n, rates, dt, i, matrix_func, self.transfer_rates, self.external_source_rates) + # filter-out negative atom densities (non-physical) + results = [np.maximum(r, 0.0) for r in results] return time.time() - start, results @abstractmethod diff --git a/tests/unit_tests/test_deplete_cram.py b/tests/unit_tests/test_deplete_cram.py index 8987fbd7aa6..e5051492770 100644 --- a/tests/unit_tests/test_deplete_cram.py +++ b/tests/unit_tests/test_deplete_cram.py @@ -3,10 +3,13 @@ Compares a few Mathematica matrix exponentials to CRAM16/CRAM48. """ +from unittest.mock import patch, MagicMock + from pytest import approx import numpy as np import scipy.sparse as sp from openmc.deplete.cram import CRAM16, CRAM48 +from openmc.deplete.abc import Integrator def test_CRAM16(): @@ -35,3 +38,32 @@ def test_CRAM48(): z0 = np.array((0.904837418035960, 0.576799023327476)) assert z == approx(z0) + + +def test_timed_deplete_clips_negatives(): + """Verify _timed_deplete clips negative atom densities to zero.""" + integrator = MagicMock(spec=Integrator) + integrator._timed_deplete = Integrator._timed_deplete.__get__(integrator) + integrator._solver = MagicMock() + integrator.chain = MagicMock() + integrator.transfer_rates = None + integrator.external_source_rates = None + + fake_results = [ + np.array([1.0e10, -1.0e-35, 0.0, -1.0e-47, 5.0e8]), + np.array([-1.0e-20, 3.0e12, -1.0e-40]), + ] + with patch('openmc.deplete.abc.deplete', return_value=fake_results): + _, results = integrator._timed_deplete(None, None, 1.0) + + for r in results: + assert np.all(r >= 0.0), f"Negative values found: {r[r < 0]}" + # Positive values unchanged + assert results[0][0] == 1.0e10 + assert results[0][4] == 5.0e8 + assert results[1][1] == 3.0e12 + # Negatives clipped to zero + assert results[0][1] == 0.0 + assert results[0][3] == 0.0 + assert results[1][0] == 0.0 + assert results[1][2] == 0.0 From c9f45cb1e85a6d50fb29b7bea26495e9b64831a2 Mon Sep 17 00:00:00 2001 From: Perry Date: Thu, 12 Mar 2026 22:29:52 -0700 Subject: [PATCH 2/3] =?UTF-8?q?=E2=97=8F=20Add=20opt-in=20clipping=20of=20?= =?UTF-8?q?negative=20atom=20densities=20after=20CRAM=20solve?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit CRAM rational approximations can undershoot zero for nearly-depleted nuclides. Add clip_negative_densities parameter (default False) to Integrator and SIIntegrator that clips results to zero with a one-time warning reporting the worst nuclide. - New bool param on Integrator.__init__ and SIIntegrator.__init__ - Warning identifies nuclide with most negative density --- openmc/deplete/abc.py | 31 ++++++++++++++++++++++++--- tests/unit_tests/test_deplete_cram.py | 2 ++ 2 files changed, 30 insertions(+), 3 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 3c0e8e4c4d3..75837792d9b 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -583,6 +583,11 @@ class Integrator(ABC): `source_rates` should be the same as the initial run. .. versionadded:: 0.15.1 + clip_negative_densities : bool, optional + Whether to clip negative atom densities to zero after each CRAM + solve. CRAM can produce small negative values. Defaults to False. + + .. versionadded:: 0.15.4 Attributes ---------- @@ -632,6 +637,7 @@ def __init__( timestep_units: str = 's', solver: str = "cram48", continue_timesteps: bool = False, + clip_negative_densities: bool = False, ): if continue_timesteps and operator.prev_res is None: raise ValueError("Continuation run requires passing prev_results.") @@ -684,6 +690,8 @@ def __init__( self.timesteps = np.asarray(seconds) self.source_rates = np.asarray(source_rates) + self.clip_negative_densities = clip_negative_densities + self._warned_negative_density = False self.transfer_rates = None self.external_source_rates = None @@ -736,8 +744,17 @@ def _timed_deplete(self, n, rates, dt, i=None, matrix_func=None): results = deplete( self._solver, self.chain, n, rates, dt, i, matrix_func, self.transfer_rates, self.external_source_rates) - # filter-out negative atom densities (non-physical) - results = [np.maximum(r, 0.0) for r in results] + if self.clip_negative_densities: + for j, r in enumerate(results): + neg_mask = r < 0 + if neg_mask.any(): + if not self._warned_negative_density: + idx = np.argmin(r) + nuc = self.chain.nuclides[idx].name + warn(f"Clipping negative atom densities " + f"(worst: {nuc} = {r[idx]:.4e})") + self._warned_negative_density = True + results[j] = np.maximum(r, 0.0) return time.time() - start, results @abstractmethod @@ -1118,6 +1135,11 @@ class SIIntegrator(Integrator): `source_rates` should be the same as the initial run. .. versionadded:: 0.15.1 + clip_negative_densities : bool, optional + Whether to clip negative atom densities to zero after each CRAM + solve. CRAM can produce small negative values. Defaults to False. + + .. versionadded:: 0.15.4 Attributes ---------- @@ -1161,12 +1183,15 @@ def __init__( n_steps: int = 10, solver: str = "cram48", continue_timesteps: bool = False, + clip_negative_densities: bool = False, ): check_type("n_steps", n_steps, Integral) check_greater_than("n_steps", n_steps, 0) super().__init__( operator, timesteps, power, power_density, source_rates, - timestep_units=timestep_units, solver=solver, continue_timesteps=continue_timesteps) + timestep_units=timestep_units, solver=solver, + continue_timesteps=continue_timesteps, + clip_negative_densities=clip_negative_densities) self.n_steps = n_steps def _get_bos_data_from_operator(self, step_index, step_power, n_bos): diff --git a/tests/unit_tests/test_deplete_cram.py b/tests/unit_tests/test_deplete_cram.py index e5051492770..502db73c95c 100644 --- a/tests/unit_tests/test_deplete_cram.py +++ b/tests/unit_tests/test_deplete_cram.py @@ -48,6 +48,8 @@ def test_timed_deplete_clips_negatives(): integrator.chain = MagicMock() integrator.transfer_rates = None integrator.external_source_rates = None + integrator.clip_negative_densities = True + integrator._warned_negative_density = False fake_results = [ np.array([1.0e10, -1.0e-35, 0.0, -1.0e-47, 5.0e8]), From 23746f1c40a9aac7530d63dfc5e3974da45c16ff Mon Sep 17 00:00:00 2001 From: Perry Date: Fri, 13 Mar 2026 22:48:03 -0700 Subject: [PATCH 3/3] Remove Test --- tests/unit_tests/test_deplete_cram.py | 34 --------------------------- 1 file changed, 34 deletions(-) diff --git a/tests/unit_tests/test_deplete_cram.py b/tests/unit_tests/test_deplete_cram.py index 502db73c95c..8987fbd7aa6 100644 --- a/tests/unit_tests/test_deplete_cram.py +++ b/tests/unit_tests/test_deplete_cram.py @@ -3,13 +3,10 @@ Compares a few Mathematica matrix exponentials to CRAM16/CRAM48. """ -from unittest.mock import patch, MagicMock - from pytest import approx import numpy as np import scipy.sparse as sp from openmc.deplete.cram import CRAM16, CRAM48 -from openmc.deplete.abc import Integrator def test_CRAM16(): @@ -38,34 +35,3 @@ def test_CRAM48(): z0 = np.array((0.904837418035960, 0.576799023327476)) assert z == approx(z0) - - -def test_timed_deplete_clips_negatives(): - """Verify _timed_deplete clips negative atom densities to zero.""" - integrator = MagicMock(spec=Integrator) - integrator._timed_deplete = Integrator._timed_deplete.__get__(integrator) - integrator._solver = MagicMock() - integrator.chain = MagicMock() - integrator.transfer_rates = None - integrator.external_source_rates = None - integrator.clip_negative_densities = True - integrator._warned_negative_density = False - - fake_results = [ - np.array([1.0e10, -1.0e-35, 0.0, -1.0e-47, 5.0e8]), - np.array([-1.0e-20, 3.0e12, -1.0e-40]), - ] - with patch('openmc.deplete.abc.deplete', return_value=fake_results): - _, results = integrator._timed_deplete(None, None, 1.0) - - for r in results: - assert np.all(r >= 0.0), f"Negative values found: {r[r < 0]}" - # Positive values unchanged - assert results[0][0] == 1.0e10 - assert results[0][4] == 5.0e8 - assert results[1][1] == 3.0e12 - # Negatives clipped to zero - assert results[0][1] == 0.0 - assert results[0][3] == 0.0 - assert results[1][0] == 0.0 - assert results[1][2] == 0.0