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

NJOY ERRORR module not working with energy grid that start at 0 ev and adding tests to check the safeguards for self-covariance matrices

Browse source
This commit is contained in:
Grego01 2026-04-16 11:29:31 -04:00
parent e2c28ebb3f
commit f7d58fb315
2 changed files with 390 additions and 91 deletions
Download patch file Download diff file Expand all files Collapse all files

9
openmc/data/mf33_njoy.py
View file

@ -51,14 +51,17 @@ def _validate_energy_grid_ev(ek: Sequence[float]) -> List[float]:
ek = [float(x) for x in ek]
if len(ek) < 2:
raise ValueError("Energy grid must have at least 2 boundaries (G+1).")
if ek[0] <= 0.0:
raise ValueError(
f"Energy grid lower boundary must be positive (got {ek[0]:g} eV). "
f"ERRORR cannot integrate from zero energy. Use a small positive "
f"value like 1e-5 eV instead."
)
for i in range(1, len(ek)):
if not (ek[i] > ek[i-1]):
raise ValueError("Energy grid boundaries must be strictly increasing (in eV).")
return ek
# ------------------------- NJOY deck builder -------------------------
def _moder_input(nin: int, nout: int) -> str:

472
tests/unit_tests/test_xs_covariance.py
View file

@ -16,33 +16,40 @@ import pytest
from openmc.data.xs_covariance_njoy import (
CovFactorResult,
NeutronXSCovariances,
RawBlockDiagnostics,
compute_covariance_factor,
validate_raw_mf33_reactions,
)
# ---------------------------------------------------------------------------
# Helpers
# ---------------------------------------------------------------------------
def _make_posdef_matrix(n, seed=42):
# Default group count and energy grid used by validate_raw_mf33_reactions tests
G = 4
EK = np.logspace(-5, 7, G + 1).tolist() # G+1 boundaries in eV
def _posdef(n, seed=42):
"""Return a random (n x n) symmetric positive-definite matrix."""
rng = np.random.default_rng(seed)
A = rng.standard_normal((n, n))
return A @ A.T + 0.1 * np.eye(n)
def _make_semidef_matrix(n, rank, seed=42):
def _semidef(n, rank, seed=42):
"""Return a random (n x n) symmetric PSD matrix with given rank."""
rng = np.random.default_rng(seed)
B = rng.standard_normal((n, rank))
return B @ B.T
def _make_mock_covariances(G=10, seed=42):
"""Build a NeutronXSCovariances object with synthetic data."""
energy = np.logspace(-5, 7, G + 1)
cov_mat = _make_posdef_matrix(G, seed=seed)
cov_mat = _posdef(G, seed=seed)
result = compute_covariance_factor(cov_mat)
reactions = {
2: {
"MAT": 2631, "MF": 33, "MT": 2,
@ -51,7 +58,7 @@ def _make_mock_covariances(G=10, seed=42):
},
}
factors = {2: {2: result}}
return NeutronXSCovariances(
name="Fe56",
energy_grid_ev=energy,
@ -60,144 +67,151 @@ def _make_mock_covariances(G=10, seed=42):
temperature_k=293.6,
factor_results=factors,
)
# ---------------------------------------------------------------------------
# Tests — compute_covariance_factor
# ---------------------------------------------------------------------------
def _wrap_self(mt, M):
"""Wrap a single self-block into the reactions dict structure."""
return {mt: {"COVS": {mt: M}}}
def _wrap_pair(mt, mt1, M_self, M_cross,
M_self1=None, M_partner=None):
"""Wrap a (mt, mt1) pair with cross-block(s)."""
reactions = {
mt: {"COVS": {mt: M_self, mt1: M_cross}},
mt1: {"COVS": {}},
}
if M_self1 is not None:
reactions[mt1]["COVS"][mt1] = M_self1
if M_partner is not None:
reactions[mt1]["COVS"][mt] = M_partner
return reactions
# ===========================================================================
# compute_covariance_factor
# ===========================================================================
def test_factor_posdef_cholesky():
"""Positive-definite input should use the Cholesky path."""
Sigma = _make_posdef_matrix(5)
Sigma = _posdef(5)
result = compute_covariance_factor(Sigma)
assert result.method == "cholesky"
assert result.effective_rank == 5
assert result.full_size == 5
assert result.L.shape == (5, 5)
np.testing.assert_allclose(result.L @ result.L.T, Sigma, atol=1e-10)
def test_factor_semidef_eigen_qr():
"""Rank-deficient input should fall back to eigen_qr."""
Sigma = _make_semidef_matrix(8, rank=3)
Sigma = _semidef(8, rank=3)
result = compute_covariance_factor(Sigma)
assert result.method == "eigen_qr"
assert result.effective_rank == 3
assert result.full_size == 8
assert result.L.shape == (8, 3)
np.testing.assert_allclose(result.L @ result.L.T, Sigma, atol=1e-8)
def test_factor_zero_matrix():
"""All-zero input should return rank-0 with the zero_matrix method."""
result = compute_covariance_factor(np.zeros((6, 6)))
assert result.method == "zero_matrix"
assert result.effective_rank == 0
assert result.full_size == 6
assert result.L.shape == (6, 0)
def test_factor_partial_zero_variance():
"""Matrix with some zero-diagonal rows should still factorise the
non-zero sub-block correctly."""
Sigma = np.zeros((5, 5))
small = _make_posdef_matrix(3, seed=99)
small = _posdef(3, seed=99)
Sigma[1:4, 1:4] = small
result = compute_covariance_factor(Sigma)
assert result.effective_rank == 3
assert result.L.shape[0] == 5
np.testing.assert_allclose(result.L @ result.L.T, Sigma, atol=1e-10)
def test_factor_symmetry_enforced():
"""Slightly asymmetric input should not crash — symmetry is forced."""
Sigma = _make_posdef_matrix(4)
Sigma = _posdef(4)
Sigma[0, 1] += 1e-12
result = compute_covariance_factor(Sigma)
assert result.effective_rank == 4
def test_cov_factor_result_fields():
"""CovFactorResult should have exactly the four documented fields."""
r = CovFactorResult(
L=np.eye(2), effective_rank=2, full_size=2, method="test",
)
assert r.L.shape == (2, 2)
assert r.effective_rank == 2
assert r.full_size == 2
assert r.method == "test"
# ---------------------------------------------------------------------------
# Tests — HDF5 round-trip (standalone file)
# ---------------------------------------------------------------------------
# ===========================================================================
# HDF5 round-trip (standalone file)
# ===========================================================================
def test_hdf5_roundtrip_with_raw_covariance(tmp_path):
"""Full round-trip storing both raw covariance and L factors."""
cov = _make_mock_covariances(G=10)
h5_path = tmp_path / "cov_test.h5"
cov.to_hdf5(h5_path, store_raw_covariance=True)
cov_read = NeutronXSCovariances.from_hdf5(h5_path)
# Metadata
assert cov_read.mat == cov.mat
assert cov_read.temperature_k == cov.temperature_k
np.testing.assert_allclose(cov_read.energy_grid_ev, cov.energy_grid_ev)
# Reaction keys preserved
assert set(cov_read.reactions.keys()) == set(cov.reactions.keys())
# Raw covariance matrices match
for mt in cov.reactions:
for mt1 in cov.reactions[mt]["COVS"]:
original = cov.reactions[mt]["COVS"][mt1]
loaded = cov_read.reactions[mt]["COVS"][mt1]
np.testing.assert_allclose(loaded, original, atol=1e-12)
# L factors match
for mt in cov.factor_results:
for mt1 in cov.factor_results[mt]:
L_orig = cov.factor_results[mt][mt1].L
L_read = cov_read.factor_results[mt][mt1].L
np.testing.assert_allclose(L_read, L_orig, atol=1e-12)
def test_hdf5_roundtrip_factors_only(tmp_path):
"""Round-trip with store_raw_covariance=False.
The raw covariance should be reconstructed as L @ L.T on read.
"""
cov = _make_mock_covariances(G=8)
h5_path = tmp_path / "cov_factors_only.h5"
cov.to_hdf5(h5_path, store_raw_covariance=False)
cov_read = NeutronXSCovariances.from_hdf5(h5_path)
# Verify "reactions" group is absent in the HDF5
with h5py.File(h5_path, "r") as f:
assert "reactions" not in f["mf33"]
# COVS should be reconstructed from L factors
for mt in cov.reactions:
for mt1 in cov.reactions[mt]["COVS"]:
original = cov.reactions[mt]["COVS"][mt1]
reconstructed = cov_read.reactions[mt]["COVS"][mt1]
np.testing.assert_allclose(reconstructed, original, atol=1e-8)
def test_hdf5_schema_attributes(tmp_path):
"""Verify the expected HDF5 attributes and datasets exist."""
cov = _make_mock_covariances(G=5)
h5_path = tmp_path / "schema_check.h5"
cov.to_hdf5(h5_path)
with h5py.File(h5_path, "r") as f:
mf33 = f["mf33"]
assert mf33.attrs["format"] == b"openmc.mf33.v1"
@ -206,18 +220,18 @@ def test_hdf5_schema_attributes(tmp_path):
assert "energy_grid_ev" in mf33
assert "mts" in mf33
assert "factors" in mf33
# Energy grid has G+1 entries
assert mf33["energy_grid_ev"].shape == (6,)
def test_hdf5_multiple_reactions(tmp_path):
"""Round-trip with multiple MT values."""
G = 6
energy = np.logspace(-5, 7, G + 1)
cov_2 = _make_posdef_matrix(G, seed=10)
cov_102 = _make_posdef_matrix(G, seed=20)
G_local = 6
energy = np.logspace(-5, 7, G_local + 1)
cov_2 = _posdef(G_local, seed=10)
cov_102 = _posdef(G_local, seed=20)
reactions = {
2: {"MAT": 2631, "MF": 33, "MT": 2,
"ZA": 26056.0, "AWR": 55.45,
@ -233,7 +247,7 @@ def test_hdf5_multiple_reactions(tmp_path):
h5_path = tmp_path / "multi_mt.h5"
cov_obj.to_hdf5(h5_path)
cov_read = NeutronXSCovariances.from_hdf5(h5_path)
assert set(cov_read.reactions.keys()) == {2, 102}
np.testing.assert_allclose(
cov_read.reactions[2]["COVS"][2], cov_2, atol=1e-12,
@ -241,42 +255,324 @@ def test_hdf5_multiple_reactions(tmp_path):
np.testing.assert_allclose(
cov_read.reactions[102]["COVS"][102], cov_102, atol=1e-12,
)
# ---------------------------------------------------------------------------
# Tests — write_mf33_group / _read_mf33_group (embedded in another HDF5)
# ---------------------------------------------------------------------------
# ===========================================================================
# write_mf33_group / _read_mf33_group (embedded in another HDF5)
# ===========================================================================
def test_write_read_mf33_into_existing_group(tmp_path):
"""Write mf33 into a pre-existing HDF5 group, then read it back."""
cov = _make_mock_covariances(G=5)
h5_path = tmp_path / "embedded.h5"
with h5py.File(h5_path, "w") as f:
nuc = f.create_group("Fe56")
cov_root = nuc.create_group("covariance")
cov.write_mf33_group(cov_root, store_raw_covariance=True)
with h5py.File(h5_path, "r") as f:
mf33_group = f["Fe56"]["covariance"]["mf33"]
cov_read = NeutronXSCovariances._read_mf33_group(
mf33_group, name="Fe56",
)
assert cov_read.name == "Fe56"
assert cov_read.mat == 2631
assert 2 in cov_read.reactions
def test_overwrite_existing_mf33(tmp_path):
"""Calling write_mf33_group twice should replace, not duplicate."""
cov = _make_mock_covariances(G=4)
h5_path = tmp_path / "overwrite.h5"
with h5py.File(h5_path, "w") as f:
root = f.create_group("root")
cov.write_mf33_group(root)
cov.write_mf33_group(root)
with h5py.File(h5_path, "r") as f:
assert "mf33" in f["root"]
assert "mf33" in f["root"]
# ===========================================================================
# validate_raw_mf33_reactions — hard failures (shape & finiteness)
# ===========================================================================
def test_validate_wrong_shape_fails():
M = np.eye(G + 1) # deliberately wrong
diags = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
dx = diags[2][2]
assert dx.status == "fail"
assert any("Wrong shape" in m for m in dx.messages)
# Downstream diagnostics untouched (None) because the check `continue`s
assert dx.n_negative_eigenvalues is None
assert dx.max_abs_correlation is None
def test_validate_nan_in_matrix_fails():
M = np.eye(G)
M[1, 1] = np.nan
diags = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
dx = diags[2][2]
assert dx.status == "fail"
assert any("NaN or Inf" in m for m in dx.messages)
assert dx.n_negative_eigenvalues is None
# ===========================================================================
# validate_raw_mf33_reactions — clean self-block baseline
# ===========================================================================
def test_validate_clean_selfblock_passes():
# Small variances (~0.01) so rel_std ~ 0.1 → well below warn threshold
M = 0.01 * _posdef(G, seed=1)
diags = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
dx = diags[2][2]
assert dx.status == "pass"
assert dx.mt == 2 and dx.mt1 == 2
assert dx.is_self_covariance is True
# Diagnostic fields should be populated for self-blocks
assert dx.n_negative_eigenvalues is not None
assert dx.max_abs_correlation is not None
assert dx.max_rel_std is not None
assert dx.n_negative_eigenvalues == 0
assert dx.max_abs_correlation <= 1.0 + 1e-10
# ===========================================================================
# validate_raw_mf33_reactions — symmetry
# ===========================================================================
def test_validate_asymmetric_selfblock_fails():
M = _posdef(G, seed=2)
M[0, 1] += 0.5 # break symmetry noticeably
diags = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
dx = diags[2][2]
assert dx.status == "fail"
assert any("Not symmetric" in m for m in dx.messages)
# ===========================================================================
# validate_raw_mf33_reactions — diagonal: negative / zero
# ===========================================================================
def test_validate_negative_diagonal_fails():
M = np.eye(G) * 0.01
M[2, 2] = -0.01
diags = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
dx = diags[2][2]
assert dx.status == "fail"
assert any("negative diagonal" in m for m in dx.messages)
def test_validate_zero_diagonal_warns():
# One group has zero variance; rest are fine.
M = np.eye(G) * 0.01
M[1, 1] = 0.0
diags = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
dx = diags[2][2]
assert dx.status == "warn"
assert any("zero-variance" in m for m in dx.messages)
def test_validate_negative_and_zero_diagonal_is_fail_not_warn():
"""If both negative and zero diagonals are present, fail wins."""
M = np.eye(G) * 0.01
M[0, 0] = 0.0
M[1, 1] = -0.01
diags = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
assert diags[2][2].status == "fail"
# ===========================================================================
# validate_raw_mf33_reactions — eigenvalue spectrum
# ===========================================================================
def test_validate_small_negative_eigenvalue_is_informational_only():
"""neg_ratio well below warn threshold (1e-3) → status stays pass, but
an informational message is appended."""
rng = np.random.default_rng(5)
Q, _ = np.linalg.qr(rng.standard_normal((G, G)))
eigs = np.array([1.0, 0.5, 0.2, -1e-6])
M = (Q * eigs) @ Q.T
M = 0.5 * (M + M.T) # numerical symmetry
diags = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
dx = diags[2][2]
# Diagonal might not all be positive due to random rotation;
# check only the eigenvalue-specific behavior.
if dx.status == "pass":
assert dx.n_negative_eigenvalues >= 1
assert dx.negativity_ratio is not None
assert dx.negativity_ratio < 1e-3
assert any("small negative" in m for m in dx.messages)
def test_validate_moderate_negative_eigenvalue_warns():
"""neg_ratio between warn (1e-3) and fail (1e-1) → warn."""
eigs = np.array([1.0, 0.5, 0.2, -1e-2]) # ratio = 1e-2
rng = np.random.default_rng(7)
Q, _ = np.linalg.qr(rng.standard_normal((G, G)))
M = (Q * eigs) @ Q.T
M = 0.5 * (M + M.T)
diags = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
dx = diags[2][2]
# Filter out the (possibly flaky) diagonal/correlation side-effects:
assert dx.n_negative_eigenvalues == 1
assert 1e-3 <= dx.negativity_ratio < 1e-1
assert any("negativity ratio" in m and "threshold" in m for m in dx.messages)
def test_validate_large_negative_eigenvalue_fails():
"""neg_ratio >= fail threshold (1e-1) → fail."""
eigs = np.array([1.0, 0.5, 0.2, -0.5]) # ratio = 0.5
rng = np.random.default_rng(11)
Q, _ = np.linalg.qr(rng.standard_normal((G, G)))
M = (Q * eigs) @ Q.T
M = 0.5 * (M + M.T)
diags = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
dx = diags[2][2]
assert dx.status == "fail"
assert dx.n_negative_eigenvalues == 1
assert dx.negativity_ratio >= 1e-1
assert any("repaired matrix will differ" in m for m in dx.messages)
# ===========================================================================
# validate_raw_mf33_reactions — correlation bound |rho| <= 1
# ===========================================================================
def test_validate_correlation_overshoot_warns():
"""Construct a symmetric matrix with an off-diagonal that exceeds the
geometric-mean bound sqrt(d_ii * d_jj)."""
M = np.eye(G) * 1.0
M[0, 1] = M[1, 0] = 1.5 # rho = 1.5 > 1
diags = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
dx = diags[2][2]
# Status is at least warn; could be fail because this also produces a
# negative eigenvalue (2x2 block has det < 0).
assert dx.status in ("warn", "fail")
assert dx.max_abs_correlation > 1.0
assert any("Correlation matrix" in m and "rho" in m for m in dx.messages)
# ===========================================================================
# validate_raw_mf33_reactions — relative uncertainty magnitude
# ===========================================================================
def test_validate_rel_std_above_warn_warns():
"""sqrt(diag) > rel_std_warn (1.0) but < rel_std_fail (10.0) → warn."""
M = np.eye(G) * 0.01
M[0, 0] = 4.0 # sigma_rel = 2.0 → between 1 and 10
diags = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
dx = diags[2][2]
assert dx.status == "warn"
assert dx.n_groups_rel_std_above_1 == 1
assert dx.n_groups_rel_std_above_10 == 0
assert any("Relative uncertainty" in m for m in dx.messages)
def test_validate_rel_std_above_fail_threshold_flagged_as_artifact():
"""sqrt(diag) > rel_std_fail (10.0) → warn with 'artifact' message."""
M = np.eye(G) * 0.01
M[0, 0] = 200.0 # sigma_rel ~ 14 > 10
diags = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
dx = diags[2][2]
assert dx.status == "warn"
assert dx.n_groups_rel_std_above_10 >= 1
assert any("processing artifact" in m for m in dx.messages)
# ===========================================================================
# validate_raw_mf33_reactions — threshold customization
# ===========================================================================
def test_validate_custom_negativity_thresholds():
"""Verify that relaxing thresholds turns a fail into a non-fail."""
eigs = np.array([1.0, 0.5, 0.2, -0.5])
rng = np.random.default_rng(11)
Q, _ = np.linalg.qr(rng.standard_normal((G, G)))
M = (Q * eigs) @ Q.T
M = 0.5 * (M + M.T)
# With defaults: fail.
strict = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
assert strict[2][2].status == "fail"
# With extremely loose thresholds, the eigenvalue-spectrum branch no
# longer fires. Diagonal / correlation checks may still warn depending
# on the random rotation, but no "repaired matrix" message.
loose = validate_raw_mf33_reactions(
_wrap_self(2, M), EK,
negativity_ratio_warn=10.0,
negativity_ratio_fail=100.0,
)
assert not any(
"repaired matrix will differ" in m
for m in loose[2][2].messages
)
# ===========================================================================
# validate_raw_mf33_reactions — cross-block partner consistency
# ===========================================================================
def test_validate_cross_block_partner_absent_passes():
"""ERRORR stores only the upper triangle; missing partner is expected."""
Mself = 0.01 * _posdef(G, seed=21)
Mcross = np.full((G, G), 0.001)
reactions = _wrap_pair(2, 102, Mself, Mcross) # no partner at (102, 2)
diags = validate_raw_mf33_reactions(reactions, EK)
dx = diags[2][102]
assert dx.is_self_covariance is False
assert dx.status == "pass"
# Cross-blocks shouldn't populate self-only fields
assert dx.n_negative_eigenvalues is None
assert dx.max_abs_correlation is None
def test_validate_cross_block_partner_consistent_passes():
Mself = 0.01 * _posdef(G, seed=22)
Mself1 = 0.01 * _posdef(G, seed=23)
Mcross = np.full((G, G), 0.001)
reactions = _wrap_pair(
2, 102, Mself, Mcross,
M_self1=Mself1,
M_partner=Mcross.T, # correct partner
)
diags = validate_raw_mf33_reactions(reactions, EK)
assert diags[2][102].status == "pass"
assert diags[102][2].status == "pass"
def test_validate_cross_block_partner_shape_mismatch_fails():
Mself = 0.01 * _posdef(G, seed=24)
Mcross = np.full((G, G), 0.001)
bad_partner = np.full((G, G + 1), 0.001) # wrong shape
reactions = _wrap_pair(
2, 102, Mself, Mcross,
M_self1=0.01 * _posdef(G, seed=25),
M_partner=bad_partner,
)
diags = validate_raw_mf33_reactions(reactions, EK)
assert diags[2][102].status == "fail"
assert any("Partner wrong shape" in m for m in diags[2][102].messages)
def test_validate_cross_block_partner_not_transpose_fails():
Mself = 0.01 * _posdef(G, seed=28)
Mcross = np.full((G, G), 0.001)
wrong_partner = np.full((G, G), 0.002) # finite, right shape, wrong values
reactions = _wrap_pair(
2, 102, Mself, Mcross,
M_self1=0.01 * _posdef(G, seed=29),
M_partner=wrong_partner,
)
diags = validate_raw_mf33_reactions(reactions, EK)
assert diags[2][102].status == "fail"
msg = " ".join(diags[2][102].messages)
assert "!=" in msg or "C(" in msg