From e2b0b97df9cfd61a5ed4ac50f08d386e9ab4a88e Mon Sep 17 00:00:00 2001 From: Grego01 Date: Wed, 4 Mar 2026 15:21:20 -0500 Subject: [PATCH 01/22] New implementation to load covariance data from HDF5 files --- include/openmc/reaction.h | 7 + include/openmc/settings.h | 1 + openmc/data/mf33_njoy.py | 200 ++++++++++++++ openmc/data/neutron.py | 149 +++++++++++ openmc/data/xs_covariance_njoy.py | 423 ++++++++++++++++++++++++++++++ openmc/settings.py | 16 ++ src/nuclide.cpp | 100 +++++++ src/settings.cpp | 8 + 8 files changed, 904 insertions(+) create mode 100644 openmc/data/mf33_njoy.py create mode 100644 openmc/data/xs_covariance_njoy.py diff --git a/include/openmc/reaction.h b/include/openmc/reaction.h index 3314d18666f..30a25ac0e77 100644 --- a/include/openmc/reaction.h +++ b/include/openmc/reaction.h @@ -5,6 +5,7 @@ #define OPENMC_REACTION_H #include +#include #include "hdf5.h" @@ -64,6 +65,12 @@ class Reaction { bool redundant_; //!< redundant reaction? vector xs_; //!< Cross section at each temperature vector products_; //!< Reaction products + + struct CovData { + vector L; + int num_groups {0}; + }; + std::unordered_map cholesky_; // MT → Cholesky factor }; //============================================================================== diff --git a/include/openmc/settings.h b/include/openmc/settings.h index b369c99fef8..d0947000ca4 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -197,6 +197,7 @@ extern int trigger_batch_interval; //!< Batch interval for triggers extern "C" int verbosity; //!< How verbose to make output extern double weight_cutoff; //!< Weight cutoff for Russian roulette extern double weight_survive; //!< Survival weight after Russian roulette +extern bool load_covariance; //!< Load covariance data for sampling } // namespace settings diff --git a/openmc/data/mf33_njoy.py b/openmc/data/mf33_njoy.py new file mode 100644 index 00000000000..a105d46016a --- /dev/null +++ b/openmc/data/mf33_njoy.py @@ -0,0 +1,200 @@ +#!/usr/bin/env python3 +""" +Minimal NJOY/ERRORR driver to produce MF=33 (multigroup covariance) as text. +Run a short NJOY chain and return dict with keys: mat, ek, tape33. + +Assumptions +----------- +- `energy_grid_ev` is a strictly increasing sequence of group boundaries in **eV** (length G+1). +- Output is relative covariances (IRELCO=1) on an explicit group structure (IGN=1). +""" + +from __future__ import annotations + +import os +import shutil +import subprocess as sp +from pathlib import Path +from tempfile import TemporaryDirectory +from typing import Dict, Optional, Sequence, Any, Tuple, List + + +# ------------------------- small helpers ------------------------- + +def _read_mat_from_endf(endf_path: Path) -> int: + """Infer MAT from an ENDF-6 evaluation (text).""" + first_nonzero = None + with endf_path.open("r", errors="ignore") as f: + for line in f: + if len(line) < 75: + continue + try: + mat = int(line[66:70]) + mf = int(line[70:72]) + mt = int(line[72:75]) + except ValueError: + continue + if mat != 0 and first_nonzero is None: + first_nonzero = mat + if mat != 0 and mf == 1 and mt == 451: + return mat + if first_nonzero is None: + raise ValueError(f"Could not infer MAT from ENDF file: {endf_path}") + return first_nonzero + + +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).") + + 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: + return f"moder\n{nin:d} {nout:d} /\n" + + +def _reconr_input(endfin: int, pendfout: int, mat: int, tol: float, header: str) -> str: + return ( + "reconr\n" + f"{endfin:d} {pendfout:d} /\n" + f"'{header}'/\n" + f"{mat:d} 0 0 /\n" + f"{tol:g} 0. /\n" + "0/\n" + ) + + +def _broadr_input(endfin: int, pendfin: int, pendfout: int, mat: int, tol: float, temperature: float) -> str: + return ( + "broadr\n" + f"{endfin:d} {pendfin:d} {pendfout:d} /\n" + f"{mat:d} 1 0 0 0. /\n" + f"{tol:g} /\n" + f"{temperature:.1f} /\n" + "0 /\n" + ) + + +def _errorr_mf33_input( + *, + endfin: int, + pendfin: int, + errorrout: int, + mat: int, + temperature: float, + ek: Sequence[float], +) -> str: + # Fixed / stable defaults + ign = 1 # explicit group structure + iwt = 2 # constant weight + iprint = 0 # quiet + irelco = 1 # relative covariances + irespr = 1 # recommended resonance parameter processing + iread = 0 # don't restrict MT list (process all available in MF=33) + + nk = len(ek) - 1 + return ( + "errorr\n" + f"{endfin:d} {pendfin:d} 0 {errorrout:d} 0 /\n" + f"{mat:d} {ign:d} {iwt:d} {iprint:d} {irelco:d} /\n" + f"{iprint:d} {temperature:.1f} /\n" + f"{iread:d} 33 {irespr:d}/\n" + f"{nk:d} /\n" + + " ".join(f"{x:.5e}" for x in ek) + " /\n" + ) + + +def _build_deck(mat: int, ek: Sequence[float], temperature: float) -> str: + # Common convention: negative units are formatted. + # Use a minimal chain: MODER -> RECONR -> BROADR -> MODER -> ERRORR + tol = 0.001 + e = 21 # internal endf copy + p = 22 # pendf after reconr + b = 23 # pendf after broadr + out_errorr = 33 + + deck = "" + deck += _moder_input(20, -e) + deck += _reconr_input(-e, -p, mat=mat, tol=tol, header="openmc_mf33") + deck += _broadr_input(-e, -p, -b, mat=mat, tol=tol, temperature=temperature) + deck += _errorr_mf33_input(endfin=-e, pendfin=-b, errorrout=out_errorr, mat=mat, temperature=temperature, ek=ek) + deck += "stop\n" + return deck + + +def _run_njoy(deck: str, endf: Path, exe: Optional[str]) -> Dict[int, str]: + if exe is None: + exe = os.environ.get("NJOY") + if not exe: + raise ValueError("NJOY executable not provided and NJOY env var is not set.") + + with TemporaryDirectory() as td: + tmpdir = Path(td) + shutil.copy(endf, tmpdir / "tape20") + + proc = sp.Popen( + exe, + shell=True, + cwd=str(tmpdir), + stdin=sp.PIPE, + stdout=True, # sp.DEVNULL + stderr=True, #sp.DEVNULL + ) + proc.communicate(input=deck.encode()) + if proc.returncode != 0: + raise RuntimeError(f"NJOY failed with return code {proc.returncode} (run with a debug wrapper to capture tapes).") + + tp33 = tmpdir / "tape33" + if not tp33.exists(): + raise RuntimeError("NJOY ran but did not produce tape33 (ERRORR output).") + + return {33: tp33.read_text()} + +def generate_errorr_mf33( + endf_path: str | Path, + energy_grid_ev: Sequence[float], + *, + njoy_exec: Optional[str] = None, + mat: Optional[int] = None, + temperature: float = 293.6, +) -> Dict[str, Any]: + """Run NJOY/ERRORR and return MF=33 tape33 as text. + + Parameters + ---------- + endf_path + Path to ENDF-6 evaluation (text). + energy_grid_ev + Explicit group boundaries (G+1 values) in **eV**. + njoy_exec + NJOY executable/command. If omitted, uses $NJOY. + mat + MAT number. If omitted, inferred from the evaluation. + temperature + Processing temperature in K (single temperature). + + Returns + ------- + dict + Keys: "mat", "ek" (validated eV grid), "tape33" (text). + """ + endf_path = Path(endf_path).resolve() + if not endf_path.exists(): + raise FileNotFoundError(endf_path) + + ek = _validate_energy_grid_ev(energy_grid_ev) + if mat is None: + mat = _read_mat_from_endf(endf_path) + + deck = _build_deck(mat=mat, ek=ek, temperature=float(temperature)) + tapes = _run_njoy(deck=deck, endf=endf_path, exe=njoy_exec) + + return {"mat": int(mat), "ek": ek, "tape33": tapes[33]} diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index 71927cbed6c..fe2dd56fa7f 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -108,6 +108,7 @@ def __init__(self, name, atomic_number, mass_number, metastable, self.reactions = {} self._urr = {} self._resonances = None + self._mg_covariance = None def __contains__(self, mt): return mt in self.reactions @@ -228,6 +229,38 @@ def urr(self, urr): cv.check_type('probability tables', value, ProbabilityTables) self._urr = urr + @property + def mg_covariance(self): + """Multigroup cross-section covariance data (your MF=33-derived product). + + Expected schema (dict-like): + - energy_bounds: (G+1,) float + - mts: (N,) int + - matrix: (...) float + - representation: optional str, e.g. "relative" + """ + return self._mg_covariance + + @mg_covariance.setter + def mg_covariance(self, cov): + if cov is None: + self._mg_covariance = None + return + cv.check_type('mg_covariance', cov, Mapping) + for k in ('energy_bounds', 'mts', 'matrix'): + if k not in cov: + raise KeyError(f"mg_covariance missing required key '{k}'") + self._mg_covariance = cov + + # Optional alias if you already use `data.covariance = ...` somewhere + @property + def covariance(self): + return self._mg_covariance + + @covariance.setter + def covariance(self, cov): + self.mg_covariance = cov + @property def temperatures(self): return [f"{int(round(kT / K_BOLTZMANN))}K" for kT in self.kTs] @@ -428,6 +461,79 @@ def export_to_hdf5(self, path, mode='a', libver='earliest'): if self.fission_energy is not None: fer_group = g.create_group('fission_energy_release') self.fission_energy.to_hdf5(fer_group) + + # Write covariance data (per-reaction MF=33 schema) + if self._mg_covariance is not None: + cov = self._mg_covariance + cov_root = g.require_group("covariance") + + # Prefer the shared writer if available (NeutronXSCovariances) + if hasattr(cov, "write_mf33_group"): + cov.write_mf33_group(cov_root) + else: + # Fallback: write directly from duck-typed attributes + if "mf33" in cov_root: + del cov_root["mf33"] + mf33 = cov_root.create_group("mf33") + + mf33.attrs["format"] = np.bytes_("openmc.mf33.v1") + mf33.attrs["source"] = np.bytes_("njoy errorr module") + mf33.attrs["relative"] = 1 + if getattr(cov, "mat", None) is not None: + mf33.attrs["mat"] = int(cov.mat) + if getattr(cov, "temperature_k", None) is not None: + mf33.attrs["temperature_k"] = float(cov.temperature_k) + + mf33.create_dataset( + "energy_grid_ev", + data=np.asarray(cov.energy_grid_ev, dtype=np.float64), + ) + + greact = mf33.create_group("reactions") + gchol = mf33.create_group("cholesky") + + # Grab cached Cholesky factors if available + chol = getattr(cov, "cholesky_factors", None) or {} + + for mt, sec in cov.reactions.items(): + gmt = greact.create_group(str(int(mt))) + gmt.attrs["ZA"] = float(sec.get("ZA", 0.0)) + gmt.attrs["AWR"] = float(sec.get("AWR", 0.0)) + + gmt_chol = gchol.create_group(str(int(mt))) + gmt_chol.attrs["ZA"] = float(sec.get("ZA", 0.0)) + gmt_chol.attrs["AWR"] = float(sec.get("AWR", 0.0)) + + covs = sec.get("COVS", {}) + for mt1, M in covs.items(): + M_arr = np.asarray(M, dtype=np.float64) + gmt.create_dataset( + str(int(mt1)), + data=M_arr, + compression="gzip", + shuffle=True, + ) + + # Use cached L or compute on the fly + L = chol.get(mt, {}).get(mt1, None) + if L is None: + import scipy.linalg as _la + try: + L = _la.cholesky(M_arr, lower=True) + except _la.LinAlgError: + evals, V = _la.eigh(M_arr) + max_e = evals.max() if evals.max() > 0 else 1.0 + evals[evals / max_e < 1e-10] = 0.0 + evals[evals < 0.0] = 0.0 + A_tmp = V * np.sqrt(evals)[np.newaxis, :] + _, R_tmp = _la.qr(A_tmp.T, mode='economic') + L = R_tmp.T + gmt_chol.create_dataset( + str(int(mt1)), + data=np.asarray(L, dtype=np.float64), + compression="gzip", + shuffle=True, + ) @classmethod def from_hdf5(cls, group_or_filename): @@ -503,6 +609,49 @@ def from_hdf5(cls, group_or_filename): if 'fission_energy_release' in group: fer_group = group['fission_energy_release'] data.fission_energy = FissionEnergyRelease.from_hdf5(fer_group) + + # Read covariance data (per-reaction MF=33 schema) + if 'covariance' in group and 'mf33' in group['covariance']: + try: + from .xs_covariance import NeutronXSCovariances + data._mg_covariance = NeutronXSCovariances._read_mf33_group( + group['covariance']['mf33'], name=name, + ) + except ImportError: + # Fallback: read into a lightweight SimpleNamespace + import types + mf33 = group['covariance']['mf33'] + ek = np.asarray(mf33['energy_grid_ev'][...], dtype=np.float64) + reactions = {} + for mt_str, gmt in mf33['reactions'].items(): + mt = int(mt_str) + covs = {} + for mt1_str, ds in gmt.items(): + covs[int(mt1_str)] = np.asarray(ds[...], dtype=np.float64) + reactions[mt] = { + 'ZA': float(gmt.attrs.get('ZA', 0.0)), + 'AWR': float(gmt.attrs.get('AWR', 0.0)), + 'COVS': covs, + } + + # Read pre-computed Cholesky factors if present + chol = None + if 'cholesky' in mf33: + chol = {} + for mt_str, gmt_chol in mf33['cholesky'].items(): + mt = int(mt_str) + chol[mt] = {} + for mt1_str, ds in gmt_chol.items(): + chol[mt][int(mt1_str)] = np.asarray(ds[...], dtype=np.float64) + + obj = types.SimpleNamespace( + energy_grid_ev=ek, + reactions=reactions, + mat=int(mf33.attrs['mat']) if 'mat' in mf33.attrs else None, + temperature_k=float(mf33.attrs['temperature_k']) if 'temperature_k' in mf33.attrs else None, + cholesky_factors=chol, + ) + data._mg_covariance = obj return data diff --git a/openmc/data/xs_covariance_njoy.py b/openmc/data/xs_covariance_njoy.py new file mode 100644 index 00000000000..eaa7fd9b0d8 --- /dev/null +++ b/openmc/data/xs_covariance_njoy.py @@ -0,0 +1,423 @@ +"""openmc.data.xs_covariance + +Tiny model + HDF5 I/O for multigroup **cross section** covariances (ERRORR MF=33). + +This is intentionally narrow: +- neutron MF=33 only +- explicit energy grid (group edges in eV) +- **relative** covariances (what ERRORR typically produces with IRELCO=1) + +Public API +---------- +NeutronXSCovariances.from_endf(...) +NeutronXSCovariances.to_hdf5(...) +NeutronXSCovariances.from_hdf5(...) + +Everything else is internal parsing. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from pathlib import Path +from typing import Any, Dict, List, Optional, Sequence, Tuple + +import numpy as np +import scipy.linalg as la +from .mf33_njoy import generate_errorr_mf33 + + +# ----------------------------------------------------------------------------- +# ERRORR tape33 (ENDF-6 text) parser for MF=33 +# ----------------------------------------------------------------------------- + +from collections import namedtuple + +_CONT = namedtuple("CONT", "C1 C2 L1 L2 N1 N2") +_LIST = namedtuple("LIST", "C1 C2 L1 L2 NPL N2 B") + + +def _endf_int(field: str) -> int: + field = field.strip() + if not field: + return 0 + try: + return int(field) + except ValueError: + return int(float(field)) + + +def _endf_float(field: str) -> float: + s = field.strip() + if not s: + return 0.0 + if "e" in s.lower(): + return float(s) + # ENDF-style: mantissa + exponent sign+digits, no 'e' + for i in range(len(s) - 1, 0, -1): + if s[i] in "+-" and s[i - 1].isdigit(): + mant = s[:i] + exp = s[i:] + return float(f"{mant}e{exp}") + return float(s) + + +def _parse_endf_ids(line: str) -> Tuple[int, int, int]: + if len(line) < 75: + return (0, 0, 0) + return (int(line[66:70]), int(line[70:72]), int(line[72:75])) + + +def _read_cont_line(line: str) -> _CONT: + return _CONT( + _endf_float(line[0:11]), + _endf_float(line[11:22]), + _endf_int(line[22:33]), + _endf_int(line[33:44]), + _endf_int(line[44:55]), + _endf_int(line[55:66]), + ) + + +def _read_list_record(lines: List[str], ipos: int) -> Tuple[_LIST, int]: + C = _read_cont_line(lines[ipos]) + npl = int(C.N1) + ipos += 1 + vals: List[float] = [] + nlines = int(np.ceil(npl / 6.0)) + for _ in range(nlines): + ln = lines[ipos] + for j in range(6): + a = 11 * j + b = a + 11 + if a >= 66: + break + vals.append(_endf_float(ln[a:b])) + ipos += 1 + vals = vals[:npl] + return _LIST(C.C1, C.C2, C.L1, C.L2, npl, C.N2, vals), ipos + + +def _parse_mf33_section_lines(section_lines: List[str], mat: int, mt: int) -> Dict[str, Any]: + i = 0 + C0 = _read_cont_line(section_lines[i]) + i += 1 + out: Dict[str, Any] = {"MAT": mat, "MF": 33, "MT": mt, "ZA": C0.C1, "AWR": C0.C2} + + reaction_pairs: Dict[int, np.ndarray] = {} + for _ in range(int(C0.N2)): # number of reaction pairs + C = _read_cont_line(section_lines[i]) + i += 1 + mt1 = int(C.L2) + ng = int(C.N2) + M = np.zeros((ng, ng)) + while True: + L, i = _read_list_record(section_lines, i) + ngcol = int(L.L1) + grow = int(L.N2) + gcol = int(L.L2) + # LIST rows are 1-based + M[grow - 1, gcol - 1 : gcol - 1 + ngcol] = np.asarray(L.B, dtype=float) + if grow >= ng: + break + reaction_pairs[mt1] = M + + out["COVS"] = reaction_pairs + return out + + +def parse_errorr_mf33_text(tape33_text: str, mat: Optional[int] = None) -> Dict[int, Dict[str, Any]]: + lines = tape33_text.splitlines(True) + reactions: Dict[int, Dict[str, Any]] = {} + + current: List[str] = [] + current_key: Optional[Tuple[int, int, int]] = None + + def flush(): + nonlocal current, current_key + if current_key is None or not current: + current = [] + current_key = None + return + m, mf, mt = current_key + if mf == 33 and mt > 0 and (mat is None or m == mat): + reactions[mt] = _parse_mf33_section_lines(current, m, mt) + current = [] + current_key = None + + for ln in lines: + m, mf, mt = _parse_endf_ids(ln) + + if m == 0 and mf == 0 and mt == 0: + flush() + continue + if mt == 0: + flush() + continue + + key = (m, mf, mt) + if current_key is None: + current_key = key + elif key != current_key: + flush() + current_key = key + + current.append(ln) + + flush() + return reactions + + +# ----------------------------------------------------------------------------- +# Cholesky factor computation (with eigendecomposition fallback) +# ----------------------------------------------------------------------------- + +def _compute_cholesky_factor(cov_matrix: np.ndarray, tol: float = 1e-10) -> np.ndarray: + """Compute the lower-triangular Cholesky factor L such that Σ ≈ L L^T. + + If the matrix is symmetric positive-definite, standard Cholesky is used. + Otherwise, falls back to eigendecomposition with negative-eigenvalue + clamping (matching the logic in sampling_algorithm.cpp): + 1. Eigendecompose: Σ = V diag(λ) V^T + 2. Zero out eigenvalues with λ/λ_max < tol + 3. Form A = V diag(√λ) + 4. QR-decompose A^T = Q R → L = R^T + + Parameters + ---------- + cov_matrix : np.ndarray + Square covariance matrix (G x G). + tol : float + Eigenvalue tolerance: eigenvalues with λ/λ_max < tol are zeroed. + + Returns + ------- + np.ndarray + Lower-triangular factor L (G x G). + """ + try: + return la.cholesky(cov_matrix, lower=True) + except la.LinAlgError: + pass + + eigenvalues, V = la.eigh(cov_matrix) + max_eig = eigenvalues.max() if eigenvalues.max() > 0 else 1.0 + eigenvalues[eigenvalues / max_eig < tol] = 0.0 + eigenvalues[eigenvalues < 0.0] = 0.0 + + A = V * np.sqrt(eigenvalues)[np.newaxis, :] + _, R = la.qr(A.T, mode='economic') + return R.T + + +# ----------------------------------------------------------------------------- +# Data model + HDF5 I/O +# ----------------------------------------------------------------------------- + +@dataclass +class NeutronXSCovariances: + """Neutron multigroup **relative** covariance data for cross sections (MF=33).""" + + name: str + energy_grid_ev: np.ndarray # shape (G+1,) + reactions: Dict[int, Dict[str, Any]] # MT -> parsed dict from ERRORR + mat: Optional[int] = None + temperature_k: Optional[float] = None # what we processed at (single T) + cholesky_factors: Optional[Dict[int, Dict[int, np.ndarray]]] = None # MT -> MT1 -> L + + @classmethod + def from_endf( + cls, + endf_path: str | Path, + energy_grid_ev: Sequence[float], + *, + njoy_exec: Optional[str] = None, + mat: Optional[int] = None, + temperature: float = 293.6, + name: Optional[str] = None, + ) -> "NeutronXSCovariances": + res = generate_errorr_mf33( + endf_path, + energy_grid_ev, + njoy_exec=njoy_exec, + mat=mat, + temperature=temperature, + ) + mat_used = int(res["mat"]) + ek = np.asarray(res["ek"], dtype=float) + + reactions = parse_errorr_mf33_text(res["tape33"], mat=mat_used) + + # Pre-compute Cholesky factors for all covariance sub-blocks + chol: Dict[int, Dict[int, np.ndarray]] = {} + for mt, sec in reactions.items(): + chol[mt] = {} + for mt1, M in sec.get("COVS", {}).items(): + chol[mt][mt1] = _compute_cholesky_factor( + np.asarray(M, dtype=np.float64) + ) + + if name is None: + name = Path(endf_path).stem + + return cls( + name=str(name), + energy_grid_ev=ek, + reactions=reactions, + mat=mat_used, + temperature_k=float(temperature), + cholesky_factors=chol, + ) + + def to_hdf5(self, filename: str | Path) -> None: + """Write covariances to a standalone HDF5 file. + + Layout matches the "covariance/mf33" group that + :meth:'IncidentNeutron.export_to_hdf5' produces, so that both + paths share one schema and one reader. + """ + import h5py + + filename = Path(filename) + filename.parent.mkdir(parents=True, exist_ok=True) + + with h5py.File(filename, "w") as f: + self.write_mf33_group(f) + + def write_mf33_group(self, h5_group) -> None: + """Write the "mf33" sub-tree into an already-open HDF5 group. + + Parameters + ---------- + h5_group : h5py.Group + Parent group. A child group called "mf33" will be created + (or replaced) directly under it. + """ + if "mf33" in h5_group: + del h5_group["mf33"] + mf33 = h5_group.create_group("mf33") + + mf33.attrs["format"] = np.bytes_("openmc.mf33.v1") + mf33.attrs["source"] = np.bytes_("njoy errorr") + mf33.attrs["relative"] = 1 # int flag – portable across languages + if self.mat is not None: + mf33.attrs["mat"] = int(self.mat) + if self.temperature_k is not None: + mf33.attrs["temperature_k"] = float(self.temperature_k) + + mf33.create_dataset( + "energy_grid_ev", + data=np.asarray(self.energy_grid_ev, dtype=np.float64), + ) + + greact = mf33.create_group("reactions") + gchol = mf33.create_group("cholesky") + + # Compute Cholesky factors on-the-fly if not already cached + chol = self.cholesky_factors or {} + + for mt, sec in self.reactions.items(): + gmt = greact.create_group(str(int(mt))) + gmt.attrs["ZA"] = float(sec.get("ZA", 0.0)) + gmt.attrs["AWR"] = float(sec.get("AWR", 0.0)) + + gmt_chol = gchol.create_group(str(int(mt))) + gmt_chol.attrs["ZA"] = float(sec.get("ZA", 0.0)) + gmt_chol.attrs["AWR"] = float(sec.get("AWR", 0.0)) + + covs: Dict[int, np.ndarray] = sec.get("COVS", {}) + for mt1, M in covs.items(): + M_arr = np.asarray(M, dtype=np.float64) + gmt.create_dataset( + str(int(mt1)), + data=M_arr, + compression="gzip", + shuffle=True, + ) + + # Use cached L if available, otherwise compute + L = chol.get(mt, {}).get(mt1, None) + if L is None: + L = _compute_cholesky_factor(M_arr) + gmt_chol.create_dataset( + str(int(mt1)), + data=np.asarray(L, dtype=np.float64), + compression="gzip", + shuffle=True, + ) + + @classmethod + def from_hdf5(cls, filename: str | Path, name: Optional[str] = None) -> "NeutronXSCovariances": + """Read covariances back from an HDF5 file. + + Supports two layouts: + - **standalone** file written by :meth:`to_hdf5` (``mf33`` group + at root) + - **embedded** inside an OpenMC incident-neutron file + (``/covariance/mf33``) + """ + import h5py + + filename = Path(filename) + with h5py.File(filename, "r") as f: + mf33 = cls._find_mf33_group(f) + return cls._read_mf33_group(mf33, name=name or filename.stem) + + @staticmethod + def _find_mf33_group(f) -> "h5py.Group": + """Locate the ``mf33`` group regardless of nesting depth.""" + # Standalone file: /mf33 + if "mf33" in f: + return f["mf33"] + # Embedded: //covariance/mf33 + for key in f: + g = f[key] + if isinstance(g, type(f)) and "covariance" in g: # h5py.Group + cov = g["covariance"] + if "mf33" in cov: + return cov["mf33"] + raise KeyError("No mf33 group found in HDF5 file.") + + @classmethod + def _read_mf33_group(cls, mf33, *, name: str = "") -> "NeutronXSCovariances": + """Build an instance from an ``mf33`` h5py.Group.""" + mat_val = mf33.attrs.get("mat", None) + temp_val = mf33.attrs.get("temperature_k", None) + + ek = np.asarray(mf33["energy_grid_ev"][...], dtype=np.float64) + + reactions: Dict[int, Dict[str, Any]] = {} + greact = mf33["reactions"] + for mt_str, gmt in greact.items(): + mt = int(mt_str) + covs: Dict[int, np.ndarray] = {} + for mt1_str, ds in gmt.items(): + covs[int(mt1_str)] = np.asarray(ds[...], dtype=np.float64) + + reactions[mt] = { + "MAT": int(mat_val) if mat_val is not None else 0, + "MF": 33, + "MT": mt, + "ZA": float(gmt.attrs.get("ZA", 0.0)), + "AWR": float(gmt.attrs.get("AWR", 0.0)), + "COVS": covs, + } + + # Read pre-computed Cholesky factors if present + chol: Optional[Dict[int, Dict[int, np.ndarray]]] = None + if "cholesky" in mf33: + chol = {} + for mt_str, gmt_chol in mf33["cholesky"].items(): + mt = int(mt_str) + chol[mt] = {} + for mt1_str, ds in gmt_chol.items(): + chol[mt][int(mt1_str)] = np.asarray(ds[...], dtype=np.float64) + + return cls( + name=str(name), + energy_grid_ev=ek, + reactions=reactions, + mat=int(mat_val) if mat_val is not None else None, + temperature_k=float(temp_val) if temp_val is not None else None, + cholesky_factors=chol, + ) \ No newline at end of file diff --git a/openmc/settings.py b/openmc/settings.py index 3cf662e311e..06bed94afff 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -462,6 +462,7 @@ def __init__(self, **kwargs): self._max_tracks = None self._max_secondaries = None self._use_decay_photons = None + self._load_covariance = None self._random_ray = {} @@ -1396,6 +1397,15 @@ def free_gas_threshold(self, free_gas_threshold: float | None): cv.check_greater_than('free gas threshold', free_gas_threshold, 0.0) self._free_gas_threshold = free_gas_threshold + @property + def load_covariance(self): + return self._load_covariance + + @load_covariance.setter + def load_covariance(self, value): + cv.check_type('load_covariance', value, bool) + self._load_covariance = value + def _create_run_mode_subelement(self, root): elem = ET.SubElement(root, "run_mode") elem.text = self._run_mode.value @@ -1739,6 +1749,11 @@ def _create_use_decay_photons_subelement(self, root): if self._use_decay_photons is not None: element = ET.SubElement(root, "use_decay_photons") element.text = str(self._use_decay_photons).lower() + + def _create_load_covariance_subelement(self, root): + if self._load_covariance is not None: + element = ET.SubElement(root, "load_covariance") + element.text = str(self._load_covariance).lower() def _create_resonance_scattering_subelement(self, root): res = self.resonance_scattering @@ -2466,6 +2481,7 @@ def to_xml_element(self, mesh_memo=None): self._create_use_decay_photons_subelement(element) self._create_source_rejection_fraction_subelement(element) self._create_free_gas_threshold_subelement(element) + self._create_load_covariance_subelement(element) # Clean the indentation in the file to be user-readable clean_indentation(element) diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 69e603a7c66..85e4bb7b780 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -1107,6 +1107,102 @@ extern "C" size_t nuclides_size() return data::nuclides.size(); } +//============================================================================== +// Covariance loading +//============================================================================== + +void load_nuclide_covariance(const std::string& name) +{ + + // Make sure the nuclide is loaded + auto it_nuc = data::nuclide_map.find(name); + if (it_nuc == data::nuclide_map.end()) + return; + + int i_nuclide = it_nuc->second; + auto& nuc = data::nuclides[i_nuclide]; + + // Find the same neutron library file used for cross sections + LibraryKey key {Library::Type::neutron, name}; + const auto& it_lib = data::library_map.find(key); + if (it_lib == data::library_map.end()) + return; + + int idx = it_lib->second; + const auto& filename = data::libraries[idx].path_; + + // Open the standard neutron XS library file + hid_t file_id = file_open(filename, 'r'); + check_data_version(file_id); + + // Open nuclide group: "/Fe56" + if (!object_exists(file_id, name.c_str())) { + file_close(file_id); + return; + } + hid_t nuc_group = open_group(file_id, name.c_str()); + + // Walk: covariance/mf33/cholesky + if (!object_exists(nuc_group, "covariance")) { + close_group(nuc_group); + file_close(file_id); + return; + } + hid_t cov_group = open_group(nuc_group, "covariance"); + + if (!object_exists(cov_group, "mf33")) { + close_group(cov_group); + close_group(nuc_group); + file_close(file_id); + return; + } + hid_t mf33_group = open_group(cov_group, "mf33"); + + if (!object_exists(mf33_group, "cholesky")) { + close_group(mf33_group); + close_group(cov_group); + close_group(nuc_group); + file_close(file_id); + return; + } + hid_t chol_group = open_group(mf33_group, "cholesky"); + + write_message(6, "Reading covariance data for {} from {}", name, filename); + + // For each reaction, load diagonal block: cholesky/{MT}/{MT} + for (auto& rx : nuc->reactions_) { + + std::string mt_str = std::to_string(rx->mt_); + + if (!object_exists(chol_group, mt_str.c_str())) + continue; + + hid_t mt_group = open_group(chol_group, mt_str.c_str()); + + if (object_exists(mt_group, mt_str.c_str())) { + // The dataset is 2D (n x n) + xt::xtensor L_matrix; + read_dataset(mt_group, mt_str.c_str(), L_matrix); + + int n = static_cast(L_matrix.shape(0)); + + // Flatten the 2D matrix into the storage vector + Reaction::CovData& cov = rx->cholesky_[rx->mt_]; + cov.L.assign(L_matrix.data(), L_matrix.data() + L_matrix.size()); + cov.num_groups = n; + } + + close_group(mt_group); + } + + // Close groups/files + close_group(chol_group); + close_group(mf33_group); + close_group(cov_group); + close_group(nuc_group); + file_close(file_id); +} + //============================================================================== // C API //============================================================================== @@ -1145,6 +1241,10 @@ extern "C" int openmc_load_nuclide(const char* name, const double* temps, int n) if (settings::temperature_multipole) read_multipole_data(i_nuclide); + // Read covariance data (Cholesky factors) if requested + if (settings::load_covariance) + load_nuclide_covariance(name); + // Read elemental data, if necessary if (settings::photon_transport) { auto element = to_element(name); diff --git a/src/settings.cpp b/src/settings.cpp index 5b472468fc4..c5ef118a6c5 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -148,6 +148,7 @@ int trigger_batch_interval {1}; int verbosity {-1}; double weight_cutoff {0.25}; double weight_survive {1.0}; +bool load_covariance {false}; } // namespace settings @@ -1268,6 +1269,13 @@ void read_settings_xml(pugi::xml_node root) settings::use_decay_photons = get_node_value_bool(root, "use_decay_photons"); } + + if (check_for_node(root, "load_covariance")) { + settings::load_covariance = + get_node_value_bool(root, "load_covariance"); + } + + } void free_memory_settings() From aff35aed020ef0a7f5bf845f7a7d3da3c8a919b6 Mon Sep 17 00:00:00 2001 From: Grego01 Date: Wed, 4 Mar 2026 15:23:53 -0500 Subject: [PATCH 02/22] Merge upstream changes --- .pre-commit-config.yaml | 7 +++++++ openmc/data/mf33_enjoy - Copy.py:Zone.Identifier | 0 openmc/data/xs_covariance_njoy - Copy.py:Zone.Identifier | 0 vendor/gsl-lite | 1 + vendor/xtensor | 1 + vendor/xtl | 1 + 6 files changed, 10 insertions(+) create mode 100644 .pre-commit-config.yaml create mode 100644 openmc/data/mf33_enjoy - Copy.py:Zone.Identifier create mode 100644 openmc/data/xs_covariance_njoy - Copy.py:Zone.Identifier create mode 160000 vendor/gsl-lite create mode 160000 vendor/xtensor create mode 160000 vendor/xtl diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 00000000000..084dd40d3b3 --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,7 @@ +repos: + - repo: https://github.com/pre-commit/mirrors-clang-format + rev: v12.0.1 # Use the version of clang-format that your project uses + hooks: + - id: clang-format + files: \\.cpp$ + args: [--style=file] \ No newline at end of file diff --git a/openmc/data/mf33_enjoy - Copy.py:Zone.Identifier b/openmc/data/mf33_enjoy - Copy.py:Zone.Identifier new file mode 100644 index 00000000000..e69de29bb2d diff --git a/openmc/data/xs_covariance_njoy - Copy.py:Zone.Identifier b/openmc/data/xs_covariance_njoy - Copy.py:Zone.Identifier new file mode 100644 index 00000000000..e69de29bb2d diff --git a/vendor/gsl-lite b/vendor/gsl-lite new file mode 160000 index 00000000000..913e86d49c6 --- /dev/null +++ b/vendor/gsl-lite @@ -0,0 +1 @@ +Subproject commit 913e86d49c6a1acca980f4e325378f9dc393493a diff --git a/vendor/xtensor b/vendor/xtensor new file mode 160000 index 00000000000..3634f2ded19 --- /dev/null +++ b/vendor/xtensor @@ -0,0 +1 @@ +Subproject commit 3634f2ded19e0cf38208c8b86cea9e1d7c8e397d diff --git a/vendor/xtl b/vendor/xtl new file mode 160000 index 00000000000..a7c1c5444df --- /dev/null +++ b/vendor/xtl @@ -0,0 +1 @@ +Subproject commit a7c1c5444dfc57f76620391af4c94785ff82c8d6 From 832d801eb9112a48eb3a04f821bbeff87ce7aed5 Mon Sep 17 00:00:00 2001 From: Grego01 Date: Wed, 4 Mar 2026 16:50:47 -0500 Subject: [PATCH 03/22] Modify xtensor to Tensor due to PR #3805 --- src/nuclide.cpp | 3 +-- vendor/gsl-lite | 1 - vendor/xtensor | 1 - vendor/xtl | 1 - 4 files changed, 1 insertion(+), 5 deletions(-) delete mode 160000 vendor/gsl-lite delete mode 160000 vendor/xtensor delete mode 160000 vendor/xtl diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 4f02cef0f53..d379062cfd5 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -1104,7 +1104,6 @@ extern "C" size_t nuclides_size() //============================================================================== // Covariance loading //============================================================================== - void load_nuclide_covariance(const std::string& name) { @@ -1175,7 +1174,7 @@ void load_nuclide_covariance(const std::string& name) if (object_exists(mt_group, mt_str.c_str())) { // The dataset is 2D (n x n) - xt::xtensor L_matrix; + tensor::Tensor L_matrix; read_dataset(mt_group, mt_str.c_str(), L_matrix); int n = static_cast(L_matrix.shape(0)); diff --git a/vendor/gsl-lite b/vendor/gsl-lite deleted file mode 160000 index 913e86d49c6..00000000000 --- a/vendor/gsl-lite +++ /dev/null @@ -1 +0,0 @@ -Subproject commit 913e86d49c6a1acca980f4e325378f9dc393493a diff --git a/vendor/xtensor b/vendor/xtensor deleted file mode 160000 index 3634f2ded19..00000000000 --- a/vendor/xtensor +++ /dev/null @@ -1 +0,0 @@ -Subproject commit 3634f2ded19e0cf38208c8b86cea9e1d7c8e397d diff --git a/vendor/xtl b/vendor/xtl deleted file mode 160000 index a7c1c5444df..00000000000 --- a/vendor/xtl +++ /dev/null @@ -1 +0,0 @@ -Subproject commit a7c1c5444dfc57f76620391af4c94785ff82c8d6 From 7a6e0bfa8ef9aa45b6196af78a6983a2f300c910 Mon Sep 17 00:00:00 2001 From: Grego01 Date: Wed, 4 Mar 2026 17:27:43 -0500 Subject: [PATCH 04/22] First implementation of sampling cross sections in OpenMC --- include/openmc/reaction.h | 7 ++++++ src/cross_sections.cpp | 7 ++++++ src/reaction.cpp | 53 +++++++++++++++++++++++++++++++++++++++ 3 files changed, 67 insertions(+) diff --git a/include/openmc/reaction.h b/include/openmc/reaction.h index e370c2eacf0..4c4a2a97c03 100644 --- a/include/openmc/reaction.h +++ b/include/openmc/reaction.h @@ -68,9 +68,16 @@ class Reaction { struct CovData { vector L; + vector e_bounds; int num_groups {0}; }; std::unordered_map cholesky_; // MT → Cholesky factor + + //! Apply multigroup perturbation factors to CE cross sections + //! \param[in] energy_grid Nuclide energy grid (the shared CE grid) + //! \param[in,out] seed Random number seed + void perturb_xs(const vector& energy_grid, uint64_t* seed); + }; //============================================================================== diff --git a/src/cross_sections.cpp b/src/cross_sections.cpp index ec9da5b8abc..d8c89e35871 100644 --- a/src/cross_sections.cpp +++ b/src/cross_sections.cpp @@ -214,6 +214,13 @@ void read_ce_cross_sections(const vector>& nuc_temps, throw std::runtime_error {openmc_err_msg}; already_read.insert(name); + + /*uint64_t xs_perturbation_seed = 42; + + auto& nuc = *data::nuclides.back(); + for (auto& rx : nuc.reactions_) { + rx->perturb_xs(nuc.grid_[0].energy, &xs_perturbation_seed); + }*/ } } diff --git a/src/reaction.cpp b/src/reaction.cpp index c02e0cc407e..d387ddeba05 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -4,6 +4,7 @@ #include #include #include // for move +#include #include @@ -179,6 +180,58 @@ double Reaction::collapse_rate(int64_t i_temp, span energy, return xs_flux_sum; } +void Reaction::perturb_xs( + const vector& energy_grid, uint64_t* seed) +{ + // Look up Cholesky factors for this reaction's MT + auto it = cholesky_.find(mt_); + if (it == cholesky_.end()) return; + + const auto& cov = it->second; + int ng = cov.num_groups; + + // 1. Sample z from N(0, I) + vector z(ng); + for (int g = 0; g < ng; ++g) { + double u1 = prn(seed); + double u2 = prn(seed); + z[g] = std::sqrt(-2.0 * std::log(u1)) * std::cos(2.0 * M_PI * u2); + } + + // 2. Correlated perturbations: delta = L * z + // L is stored column-major (lower triangular, ng x ng) + vector delta(ng, 0.0); + for (int i = 0; i < ng; ++i) + for (int j = 0; j <= i; ++j) + delta[i] += cov.L[i * ng + j] * z[j]; + + // 3. Multiplicative factors per group + vector factor(ng); + for (int g = 0; g < ng; ++g) + factor[g] = 1.0 + delta[g]; + + // 4. Apply to every temperature's CE cross section + for (auto& txs : xs_) { + int n_points = txs.value.size(); + for (int k = 0; k < n_points; ++k) { + // Map CE grid index to energy + int i_grid = k + txs.threshold; + double E = energy_grid[i_grid]; + + // Find which covariance group this energy falls in + // e_bounds is sorted ascending, size ng+1 + int g = lower_bound_index( + cov.e_bounds.cbegin(), cov.e_bounds.cend(), E); + g = std::min(g, ng - 1); + + txs.value[k] *= factor[g]; + + // Keep cross sections non-negative + if (txs.value[k] < 0.0) txs.value[k] = 0.0; + } + } +} + //============================================================================== // Non-member functions //============================================================================== From 7d6a7c6558ce48f821ad149e56389e18fb2e7496 Mon Sep 17 00:00:00 2001 From: Grego01 Date: Fri, 6 Mar 2026 11:35:25 -0500 Subject: [PATCH 05/22] Fixing bugs and more efficient loading of covariance data in OpenMC through HDF5 files --- include/openmc/reaction.h | 7 +- openmc/data/mf33_njoy.py | 280 +++++++++++++++---- openmc/data/xs_covariance_njoy.py | 428 +++++++++++++++++++++--------- src/nuclide.cpp | 41 ++- src/reaction.cpp | 37 ++- 5 files changed, 582 insertions(+), 211 deletions(-) diff --git a/include/openmc/reaction.h b/include/openmc/reaction.h index 4c4a2a97c03..e9b8cb76de7 100644 --- a/include/openmc/reaction.h +++ b/include/openmc/reaction.h @@ -67,9 +67,10 @@ class Reaction { vector products_; //!< Reaction products struct CovData { - vector L; - vector e_bounds; - int num_groups {0}; + vector L; //!< Flattened (full_size x effective_rank) factor + vector e_bounds; //!< Energy group boundaries [eV], length = full_size+1 + int full_size {0}; //!< Original matrix dimension G (number of groups) + int effective_rank {0}; //!< Number of columns retained (r <= G) }; std::unordered_map cholesky_; // MT → Cholesky factor diff --git a/openmc/data/mf33_njoy.py b/openmc/data/mf33_njoy.py index a105d46016a..64a703742e6 100644 --- a/openmc/data/mf33_njoy.py +++ b/openmc/data/mf33_njoy.py @@ -16,9 +16,13 @@ import subprocess as sp from pathlib import Path from tempfile import TemporaryDirectory -from typing import Dict, Optional, Sequence, Any, Tuple, List +from typing import Any, Dict, List, Optional, Sequence, Union +NJOY_tolerance = 0.01 +NJOY_sig0 = [1e10] +NJOY_thermr_emax = 10.0 + # ------------------------- small helpers ------------------------- def _read_mat_from_endf(endf_path: Path) -> int: @@ -34,9 +38,9 @@ def _read_mat_from_endf(endf_path: Path) -> int: mt = int(line[72:75]) except ValueError: continue - if mat != 0 and first_nonzero is None: + if mat > 0 and first_nonzero is None: first_nonzero = mat - if mat != 0 and mf == 1 and mt == 451: + if mat > 0 and mf == 1 and mt == 451: return mat if first_nonzero is None: raise ValueError(f"Could not infer MAT from ENDF file: {endf_path}") @@ -60,25 +64,57 @@ def _validate_energy_grid_ev(ek: Sequence[float]) -> List[float]: def _moder_input(nin: int, nout: int) -> str: return f"moder\n{nin:d} {nout:d} /\n" - -def _reconr_input(endfin: int, pendfout: int, mat: int, tol: float, header: str) -> str: +def _reconr_input(endfin, pendfout, mat, err=NJOY_tolerance): return ( "reconr\n" f"{endfin:d} {pendfout:d} /\n" - f"'{header}'/\n" + f"'mf33'/\n" f"{mat:d} 0 0 /\n" - f"{tol:g} 0. /\n" + f"{err:g} 0. /\n" "0/\n" ) +def _thermr_input(endfin, pendfin, pendfout, mat, temperature, + err=NJOY_tolerance, emax=NJOY_thermr_emax): + return ( + "thermr\n" + f"{endfin:d} {pendfin:d} {pendfout:d} /\n" + f"0 {mat:d} 20 1 1 0 0 1 221 0 /\n" + f"{temperature:.1f} /\n" + f"{err:g} {emax:g} /\n" + ) -def _broadr_input(endfin: int, pendfin: int, pendfout: int, mat: int, tol: float, temperature: float) -> str: +def _broadr_input(endfin, pendfin, pendfout, mat, temperature, err=NJOY_tolerance): return ( "broadr\n" f"{endfin:d} {pendfin:d} {pendfout:d} /\n" f"{mat:d} 1 0 0 0. /\n" - f"{tol:g} /\n" + f"{err:g} /\n" + f"{temperature:.1f} /\n" + "0 /\n" + ) + +def _unresr_input(endfin, pendfin, pendfout, mat, temperature, + sig0=NJOY_sig0): + nsig0 = len(sig0) + return ( + "unresr\n" + f"{endfin:d} {pendfin:d} {pendfout:d} /\n" + f"{mat:d} 1 {nsig0:d} 0 /\n" + f"{temperature:.1f} /\n" + + " ".join(f"{s:.2E}" for s in sig0) + " /\n" + "0 /\n" + ) + +def _purr_input(endfin, pendfin, pendfout, mat, temperature, + sig0=NJOY_sig0, bins=20, ladders=32): + nsig0 = len(sig0) + return ( + "purr\n" + f"{endfin:d} {pendfin:d} {pendfout:d} /\n" + f"{mat:d} 1 {nsig0:d} {bins:d} {ladders:d} 0 /\n" f"{temperature:.1f} /\n" + + " ".join(f"{s:.2E}" for s in sig0) + " /\n" "0 /\n" ) @@ -91,73 +127,155 @@ def _errorr_mf33_input( mat: int, temperature: float, ek: Sequence[float], + iwt: int = 2, + spectrum: Optional[Sequence[float]] = None, + relative: bool = True, + irespr: int = 1, + mt: Optional[Union[int, Sequence[int]]] = None, + iprint: bool = False, ) -> str: - # Fixed / stable defaults - ign = 1 # explicit group structure - iwt = 2 # constant weight - iprint = 0 # quiet - irelco = 1 # relative covariances - irespr = 1 # recommended resonance parameter processing - iread = 0 # don't restrict MT list (process all available in MF=33) - - nk = len(ek) - 1 - return ( - "errorr\n" - f"{endfin:d} {pendfin:d} 0 {errorrout:d} 0 /\n" - f"{mat:d} {ign:d} {iwt:d} {iprint:d} {irelco:d} /\n" - f"{iprint:d} {temperature:.1f} /\n" - f"{iread:d} 33 {irespr:d}/\n" - f"{nk:d} /\n" - + " ".join(f"{x:.5e}" for x in ek) + " /\n" + irelco = 1 if relative else 0 + iread = 1 if mt is not None else 0 + iwt_ = 1 if spectrum is not None else iwt + pf = int(iprint) + nk = len(ek) - 1 + + lines = [ + "errorr", + f"{endfin:d} {pendfin:d} 0 {errorrout:d} 0 /", + f"{mat:d} 1 {iwt_:d} {pf:d} {irelco} /", # ign=1 (explicit grid) + f"{pf:d} {temperature:.1f} /", + f"{iread:d} 33 {irespr:d} /", + ] + + # MT list + if iread == 1: + mtlist = [mt] if isinstance(mt, int) else list(mt) + lines.append(f"{len(mtlist):d} 0 /") + lines.append(" ".join(str(m) for m in mtlist) + " /") + + # Explicit energy grid + lines.append(f"{nk} /") + lines.append(" ".join(f"{x:.5e}" for x in ek) + " /") + + # User weight spectrum (TAB1, iwt=1) + if iwt_ == 1 and spectrum is not None: + n_pairs = len(spectrum) // 2 + lines.append(f" 0.0 0.0 0 0 1 {n_pairs:d} /") + lines.append(f"{n_pairs:d} 1 /") + for k in range(n_pairs): + lines.append(f"{spectrum[2*k]:.6e} {spectrum[2*k+1]:.6e} /") + lines.append("/") + + return "\n".join(lines) + "\n" + + +# ---- deck assembly --------------------------------------------------------- + +def _build_deck( + mat: int, + ek: Sequence[float], + temperature: float, + *, + err: float = NJOY_tolerance, + thermr: bool = False, + unresr: bool = False, + purr: bool = False, + iwt: int = 2, + spectrum: Optional[Sequence[float]] = None, + relative: bool = True, + irespr: int = 1, + mt: Optional[Union[int, Sequence[int]]] = None, + iprint: bool = False, +) -> str: + e = 21 # internal formatted ENDF + p = e + 1 # running PENDF tape number + deck = _moder_input(20, -e) + + # RECONR + deck += _reconr_input(-e, -p, mat=mat, err=err) + + # BROADR + if temperature > 0: + o = p + 1 + deck += _broadr_input(-e, -p, -o, mat=mat, + temperature=temperature, err=err) + p = o + + # THERMR (free-gas) + if thermr: + o = p + 1 + deck += _thermr_input(0, -p, -o, mat=mat, + temperature=temperature, err=err) + p = o + + # UNRESR + if unresr: + o = p + 1 + deck += _unresr_input(-e, -p, -o, mat=mat, + temperature=temperature) + p = o + + # PURR + if purr: + o = p + 1 + deck += _purr_input(-e, -p, -o, mat=mat, + temperature=temperature) + p = o + + # ERRORR (MF=33) + deck += _errorr_mf33_input( + endfin=-e, pendfin=-p, errorrout=33, + mat=mat, temperature=temperature, ek=ek, + iwt=iwt, spectrum=spectrum, relative=relative, + irespr=irespr, mt=mt, iprint=iprint, ) - - -def _build_deck(mat: int, ek: Sequence[float], temperature: float) -> str: - # Common convention: negative units are formatted. - # Use a minimal chain: MODER -> RECONR -> BROADR -> MODER -> ERRORR - tol = 0.001 - e = 21 # internal endf copy - p = 22 # pendf after reconr - b = 23 # pendf after broadr - out_errorr = 33 - - deck = "" - deck += _moder_input(20, -e) - deck += _reconr_input(-e, -p, mat=mat, tol=tol, header="openmc_mf33") - deck += _broadr_input(-e, -p, -b, mat=mat, tol=tol, temperature=temperature) - deck += _errorr_mf33_input(endfin=-e, pendfin=-b, errorrout=out_errorr, mat=mat, temperature=temperature, ek=ek) deck += "stop\n" return deck +# ---- NJOY runner ----------------------------------------------------------- + def _run_njoy(deck: str, endf: Path, exe: Optional[str]) -> Dict[int, str]: if exe is None: exe = os.environ.get("NJOY") if not exe: - raise ValueError("NJOY executable not provided and NJOY env var is not set.") + raise ValueError( + "NJOY executable not provided and $NJOY env var is unset." + ) with TemporaryDirectory() as td: tmpdir = Path(td) shutil.copy(endf, tmpdir / "tape20") + (tmpdir / "input").write_text(deck) proc = sp.Popen( - exe, - shell=True, - cwd=str(tmpdir), - stdin=sp.PIPE, - stdout=True, # sp.DEVNULL - stderr=True, #sp.DEVNULL + exe, shell=True, cwd=str(tmpdir), + stdin=sp.PIPE, stdout=sp.PIPE, stderr=sp.PIPE, ) - proc.communicate(input=deck.encode()) + stdout, stderr = proc.communicate(input=deck.encode()) + if proc.returncode != 0: - raise RuntimeError(f"NJOY failed with return code {proc.returncode} (run with a debug wrapper to capture tapes).") + # Save work dir for debugging + dest = Path("njoy_failed_outputs") + if dest.exists(): + shutil.rmtree(dest) + shutil.copytree(tmpdir, dest) + raise RuntimeError( + f"NJOY failed (rc={proc.returncode}). " + f"Work dir saved to '{dest}'.\n" + f"stderr: {stderr.decode(errors='replace')[:2000]}" + ) tp33 = tmpdir / "tape33" if not tp33.exists(): - raise RuntimeError("NJOY ran but did not produce tape33 (ERRORR output).") + raise RuntimeError("NJOY ran but did not produce tape33.") return {33: tp33.read_text()} + +# ---- public API ------------------------------------------------------------ + def generate_errorr_mf33( endf_path: str | Path, energy_grid_ev: Sequence[float], @@ -165,6 +283,19 @@ def generate_errorr_mf33( njoy_exec: Optional[str] = None, mat: Optional[int] = None, temperature: float = 293.6, + # processing chain + err: float = NJOY_tolerance, + minimal_processing: bool = True, + thermr: bool = False, + unresr: bool = False, + purr: bool = False, + # ERRORR options + iwt: int = 2, + spectrum: Optional[Sequence[float]] = None, + relative: bool = True, + irespr: int = 1, + mt: Optional[Union[int, Sequence[int]]] = None, + iprint: bool = False, ) -> Dict[str, Any]: """Run NJOY/ERRORR and return MF=33 tape33 as text. @@ -173,13 +304,35 @@ def generate_errorr_mf33( endf_path Path to ENDF-6 evaluation (text). energy_grid_ev - Explicit group boundaries (G+1 values) in **eV**. + Explicit group boundaries (G+1 values) in eV. njoy_exec - NJOY executable/command. If omitted, uses $NJOY. + NJOY executable/command. Falls back to $NJOY. mat - MAT number. If omitted, inferred from the evaluation. + MAT number. Inferred from the file when omitted. temperature - Processing temperature in K (single temperature). + Processing temperature in K. + err + Reconstruction tolerance for RECONR / BROADR. + minimal_processing + If True (default), force thermr/unresr/purr off. + thermr, unresr, purr + Individual module toggles. All ignored when + ``minimal_processing=True``. + iwt + Weight function option (default 2 = constant). + Ignored when *spectrum* is given. + spectrum + User weight function as flat [E0, W0, E1, W1, …]. + Forces iwt=1. + relative + Produce relative covariances (default True). + irespr + Resonance-parameter covariance processing + (0 = area sensitivity, 1 = 1%% sensitivity). + mt + Restrict ERRORR to specific reaction numbers. + iprint + NJOY print flag. Returns ------- @@ -194,7 +347,16 @@ def generate_errorr_mf33( if mat is None: mat = _read_mat_from_endf(endf_path) - deck = _build_deck(mat=mat, ek=ek, temperature=float(temperature)) - tapes = _run_njoy(deck=deck, endf=endf_path, exe=njoy_exec) + if minimal_processing: + thermr = unresr = purr = False + + deck = _build_deck( + mat=mat, ek=ek, temperature=float(temperature), + err=err, thermr=thermr, unresr=unresr, purr=purr, + iwt=iwt, spectrum=spectrum, relative=relative, + irespr=irespr, mt=mt, iprint=iprint, + ) + print("Running NJOY to produce MF33 covariance data") + tapes = _run_njoy(deck, endf_path, exe=njoy_exec) - return {"mat": int(mat), "ek": ek, "tape33": tapes[33]} + return {"mat": int(mat), "ek": ek, "tape33": tapes[33]} \ No newline at end of file diff --git a/openmc/data/xs_covariance_njoy.py b/openmc/data/xs_covariance_njoy.py index eaa7fd9b0d8..9c8db342ae6 100644 --- a/openmc/data/xs_covariance_njoy.py +++ b/openmc/data/xs_covariance_njoy.py @@ -18,14 +18,16 @@ from __future__ import annotations -from dataclasses import dataclass +import logging +from dataclasses import dataclass, field from pathlib import Path -from typing import Any, Dict, List, Optional, Sequence, Tuple +from typing import Any, Dict, List, Optional, Sequence, Tuple, Union import numpy as np import scipy.linalg as la from .mf33_njoy import generate_errorr_mf33 +log = logging.getLogger(__name__) # ----------------------------------------------------------------------------- # ERRORR tape33 (ENDF-6 text) parser for MF=33 @@ -172,43 +174,223 @@ def flush(): # Cholesky factor computation (with eigendecomposition fallback) # ----------------------------------------------------------------------------- -def _compute_cholesky_factor(cov_matrix: np.ndarray, tol: float = 1e-10) -> np.ndarray: - """Compute the lower-triangular Cholesky factor L such that Σ ≈ L L^T. +@dataclass +class CholeskyResult: + """Container for a Cholesky-like factorization and its diagnostics. + + Attributes + ---------- + L : np.ndarray + Lower-triangular factor, shape (G, r) where r = effective_rank. + Satisfies L @ L.T ≈ original covariance (after any regularization). + effective_rank : int + Number of positive eigenvalues retained. + full_size : int + Original matrix dimension G. + method : str + Which path was taken: "cholesky", "eigen_qr", or "zero_matrix". + condition_number : float + Ratio of largest to smallest *positive* singular value. + np.inf when the matrix is singular. + negative_eig_mass : float + Sum of absolute values of negative eigenvalues (before clamping). + Zero means the raw matrix was already PSD. + negative_eig_mass_pct : float + negative_eig_mass as a percentage of the Frobenius norm. + regularization_applied : float + The diagonal-loading factor actually used (0.0 if none). + nonzero_variance_indices : np.ndarray + Integer indices of rows/columns that had non-zero variance. + """ + L: np.ndarray + effective_rank: int + full_size: int + method: str = "eigen_qr" + condition_number: float = np.inf + negative_eig_mass: float = 0.0 + negative_eig_mass_pct: float = 0.0 + regularization_applied: float = 0.0 + nonzero_variance_indices: np.ndarray = field(default_factory=lambda: np.array([], dtype=int)) + + +def _enforce_symmetry(A: np.ndarray) -> np.ndarray: + """Force exact symmetry: A_sym = (A + A^T) / 2.""" + return 0.5 * (A + A.T) + + +def _strip_zero_variance(A: np.ndarray): + """Remove rows/columns with zero diagonal (zero variance). + + Returns + ------- + A_reduced : np.ndarray + Square sub-matrix with non-zero-variance entries only. + nz_indices : np.ndarray + Integer indices into the original matrix for the kept rows/cols. + """ + diag = np.diag(A) + nz = np.flatnonzero(diag > 0.0) + if len(nz) == 0: + return np.empty((0, 0), dtype=A.dtype), nz + return A[np.ix_(nz, nz)], nz + + +def _regularize(A: np.ndarray, correction: float) -> np.ndarray: + """Add scaled diagonal loading: A_reg = A + correction * diag(A).""" + D = np.diag(A).copy() + return A + np.diag(D * correction) + + +def _compute_diagnostics(eigenvalues: np.ndarray, fro_norm: float) -> dict: + """Compute negative-eigenvalue mass and condition number.""" + eig_neg = eigenvalues[eigenvalues < 0.0] + neg_mass = float(-np.sum(eig_neg)) + + eig_pos = eigenvalues[eigenvalues > 0.0] + if len(eig_pos) >= 2: + cond = float(eig_pos.max() / eig_pos.min()) + elif len(eig_pos) == 1: + cond = 1.0 + else: + cond = np.inf + + neg_pct = 100.0 * neg_mass / fro_norm if fro_norm > 0 else 0.0 + return dict( + negative_eig_mass=neg_mass, + negative_eig_mass_pct=neg_pct, + condition_number=cond, + ) + - If the matrix is symmetric positive-definite, standard Cholesky is used. - Otherwise, falls back to eigendecomposition with negative-eigenvalue - clamping (matching the logic in sampling_algorithm.cpp): - 1. Eigendecompose: Σ = V diag(λ) V^T - 2. Zero out eigenvalues with λ/λ_max < tol - 3. Form A = V diag(√λ) - 4. QR-decompose A^T = Q R → L = R^T +def compute_cholesky_factor( + cov_matrix: np.ndarray, + tol: float = 1e-10, + regularization: float = 0.0, +) -> CholeskyResult: + """Compute a lower-triangular factor L such that L @ L.T ≈ Sigma. Parameters ---------- cov_matrix : np.ndarray - Square covariance matrix (G x G). + Square covariance matrix (G * G). tol : float - Eigenvalue tolerance: eigenvalues with λ/λ_max < tol are zeroed. + Eigenvalue tolerance: eigenvalues with λ / λ_max < tol are + zeroed. + regularization : float + Diagonal-loading fraction (0.0 = none, 0.005 = 0.5 %). Returns ------- - np.ndarray - Lower-triangular factor L (G x G). + CholeskyResult """ + G = cov_matrix.shape[0] + + # --- 0. Enforce exact symmetry --- + A = _enforce_symmetry(np.asarray(cov_matrix, dtype=np.float64)) + + # --- 1. Strip zero-variance rows/columns --- + Ar, nz = _strip_zero_variance(A) + + if Ar.size == 0: + return CholeskyResult( + L=np.zeros((G, 0), dtype=np.float64), + effective_rank=0, + full_size=G, + method="zero_matrix", + condition_number=np.inf, + nonzero_variance_indices=nz, + ) + + # --- 2. Optional regularization --- + reg_applied = 0.0 + if regularization > 0.0: + Ar = _regularize(Ar, regularization) + reg_applied = regularization + + Gr = Ar.shape[0] + fro_norm = float(np.linalg.norm(Ar, "fro")) + + # --- 3a. Fast path: standard Cholesky --- try: - return la.cholesky(cov_matrix, lower=True) + Lr = la.cholesky(Ar, lower=True) + eig_pos = np.diag(Lr) ** 2 + cond = float(eig_pos.max() / eig_pos.min()) if eig_pos.min() > 0 else np.inf + + L_full = np.zeros((G, Gr), dtype=np.float64) + L_full[nz, :] = Lr + + return CholeskyResult( + L=L_full, + effective_rank=Gr, + full_size=G, + method="cholesky", + condition_number=cond, + negative_eig_mass=0.0, + negative_eig_mass_pct=0.0, + regularization_applied=reg_applied, + nonzero_variance_indices=nz, + ) except la.LinAlgError: pass - eigenvalues, V = la.eigh(cov_matrix) + # --- 3b. Fallback: eigendecomposition + QR --- + eigenvalues, V = la.eigh(Ar) + + diag_info = _compute_diagnostics(eigenvalues, fro_norm) + + if diag_info["negative_eig_mass"] > 0: + log.warning( + "Non-PSD block (size %d): negative eigenvalue mass = %.3e " + "(%.2f%% of Frobenius norm). Clamping to zero.", + Gr, + diag_info["negative_eig_mass"], + diag_info["negative_eig_mass_pct"], + ) + max_eig = eigenvalues.max() if eigenvalues.max() > 0 else 1.0 eigenvalues[eigenvalues / max_eig < tol] = 0.0 eigenvalues[eigenvalues < 0.0] = 0.0 - A = V * np.sqrt(eigenvalues)[np.newaxis, :] - _, R = la.qr(A.T, mode='economic') - return R.T + pos_mask = eigenvalues > 0.0 + r = int(pos_mask.sum()) + + if r == 0: + return CholeskyResult( + L=np.zeros((G, 0), dtype=np.float64), + effective_rank=0, + full_size=G, + method="eigen_qr", + condition_number=np.inf, + negative_eig_mass=diag_info["negative_eig_mass"], + negative_eig_mass_pct=diag_info["negative_eig_mass_pct"], + regularization_applied=reg_applied, + nonzero_variance_indices=nz, + ) + V_pos = V[:, pos_mask] + sqrt_lam = np.sqrt(eigenvalues[pos_mask]) + Athin = V_pos * sqrt_lam[np.newaxis, :] + + _, R = la.qr(Athin.T, mode="economic") + Lr = R.T + + L_full = np.zeros((G, r), dtype=np.float64) + L_full[nz, :] = Lr + + eig_pos_vals = eigenvalues[pos_mask] + cond = float(eig_pos_vals.max() / eig_pos_vals.min()) if eig_pos_vals.min() > 0 else np.inf + + return CholeskyResult( + L=L_full, + effective_rank=r, + full_size=G, + method="eigen_qr", + condition_number=cond, + negative_eig_mass=diag_info["negative_eig_mass"], + negative_eig_mass_pct=diag_info["negative_eig_mass_pct"], + regularization_applied=reg_applied, + nonzero_variance_indices=nz, + ) # ----------------------------------------------------------------------------- # Data model + HDF5 I/O @@ -223,7 +405,7 @@ class NeutronXSCovariances: reactions: Dict[int, Dict[str, Any]] # MT -> parsed dict from ERRORR mat: Optional[int] = None temperature_k: Optional[float] = None # what we processed at (single T) - cholesky_factors: Optional[Dict[int, Dict[int, np.ndarray]]] = None # MT -> MT1 -> L + cholesky_results: Optional[Dict[int, Dict[int, CholeskyResult]]] = None @classmethod def from_endf( @@ -235,6 +417,9 @@ def from_endf( mat: Optional[int] = None, temperature: float = 293.6, name: Optional[str] = None, + compute_cholesky: bool = True, + eig_tol: float = 1e-10, + regularization: float = 0.0, ) -> "NeutronXSCovariances": res = generate_errorr_mf33( endf_path, @@ -249,13 +434,28 @@ def from_endf( reactions = parse_errorr_mf33_text(res["tape33"], mat=mat_used) # Pre-compute Cholesky factors for all covariance sub-blocks - chol: Dict[int, Dict[int, np.ndarray]] = {} - for mt, sec in reactions.items(): - chol[mt] = {} - for mt1, M in sec.get("COVS", {}).items(): - chol[mt][mt1] = _compute_cholesky_factor( - np.asarray(M, dtype=np.float64) - ) + chol: Optional[Dict[int, Dict[int, CholeskyResult]]] = None + if compute_cholesky: + chol = {} + for mt_key, sec in reactions.items(): + chol[mt_key] = {} + for mt1, M in sec.get("COVS", {}).items(): + result = compute_cholesky_factor( + np.asarray(M, dtype=np.float64), + tol=eig_tol, + regularization=regularization, + ) + chol[mt_key][mt1] = result + log.info( + " MT %s -> MT1 %s: method=%-9s rank=%d/%d " + "cond=%.2e neg_mass=%.2e%%", + mt_key, mt1, + result.method, + result.effective_rank, + result.full_size, + result.condition_number, + result.negative_eig_mass_pct, + ) if name is None: name = Path(endf_path).stem @@ -266,10 +466,14 @@ def from_endf( reactions=reactions, mat=mat_used, temperature_k=float(temperature), - cholesky_factors=chol, + cholesky_results=chol, ) + + # ----------------------------------------------------------------- + # HDF5 writing + # ----------------------------------------------------------------- - def to_hdf5(self, filename: str | Path) -> None: + def to_hdf5(self, filename: str | Path, store_raw_covariance: bool = True) -> None: """Write covariances to a standalone HDF5 file. Layout matches the "covariance/mf33" group that @@ -282,16 +486,26 @@ def to_hdf5(self, filename: str | Path) -> None: filename.parent.mkdir(parents=True, exist_ok=True) with h5py.File(filename, "w") as f: - self.write_mf33_group(f) + self.write_mf33_group(f, store_raw_covariance=store_raw_covariance) - def write_mf33_group(self, h5_group) -> None: + def write_mf33_group(self, h5_group, store_raw_covariance: bool = True, + eig_tol: float = 1e-10, regularization: float = 0.0,) -> None: """Write the "mf33" sub-tree into an already-open HDF5 group. Parameters ---------- h5_group : h5py.Group - Parent group. A child group called "mf33" will be created + Parent group. A child group called ``mf33`` will be created (or replaced) directly under it. + store_raw_covariance + If True, store full covariance matrices under ``reactions/``. + If False, only Cholesky factors are written. + eig_tol + Eigenvalue tolerance — only used as fallback when + ``self.cholesky_results`` is None. + regularization + Diagonal-loading fraction — only used as fallback when + ``self.cholesky_results`` is None. """ if "mf33" in h5_group: del h5_group["mf33"] @@ -300,6 +514,8 @@ def write_mf33_group(self, h5_group) -> None: mf33.attrs["format"] = np.bytes_("openmc.mf33.v1") mf33.attrs["source"] = np.bytes_("njoy errorr") mf33.attrs["relative"] = 1 # int flag – portable across languages + mf33.attrs["store_raw_covariance"] = int(store_raw_covariance) + mf33.attrs["cholesky_storage"] = np.bytes_("thin_rank_r") if self.mat is not None: mf33.attrs["mat"] = int(self.mat) if self.temperature_k is not None: @@ -309,115 +525,65 @@ def write_mf33_group(self, h5_group) -> None: "energy_grid_ev", data=np.asarray(self.energy_grid_ev, dtype=np.float64), ) + mf33.create_dataset( + "mts", + data=np.asarray(sorted(self.reactions.keys()), dtype=np.int32), + ) - greact = mf33.create_group("reactions") + if store_raw_covariance: + greact = mf33.create_group("reactions") gchol = mf33.create_group("cholesky") - # Compute Cholesky factors on-the-fly if not already cached - chol = self.cholesky_factors or {} + cached = self.cholesky_results or {} for mt, sec in self.reactions.items(): - gmt = greact.create_group(str(int(mt))) - gmt.attrs["ZA"] = float(sec.get("ZA", 0.0)) - gmt.attrs["AWR"] = float(sec.get("AWR", 0.0)) + mt_str = str(int(mt)) + + if store_raw_covariance: + gmt = greact.create_group(mt_str) + for attr_name in ("ZA", "AWR"): + if attr_name in sec: + gmt.attrs[attr_name] = float(sec[attr_name]) - gmt_chol = gchol.create_group(str(int(mt))) - gmt_chol.attrs["ZA"] = float(sec.get("ZA", 0.0)) - gmt_chol.attrs["AWR"] = float(sec.get("AWR", 0.0)) + gmt_chol = gchol.create_group(mt_str) + for attr_name in ("ZA", "AWR"): + if attr_name in sec: + gmt_chol.attrs[attr_name] = float(sec[attr_name]) covs: Dict[int, np.ndarray] = sec.get("COVS", {}) for mt1, M in covs.items(): M_arr = np.asarray(M, dtype=np.float64) - gmt.create_dataset( - str(int(mt1)), - data=M_arr, + ds_name = str(int(mt1)) + + # ---- raw covariance (optional) ---- + if store_raw_covariance: + gmt.create_dataset( + ds_name, + data=M_arr, + compression="gzip", + shuffle=True, + ) + + # ---- Cholesky factor (always stored) ---- + result = cached.get(mt, {}).get(mt1, None) + if result is None: + result = compute_cholesky_factor( + M_arr, tol=eig_tol, regularization=regularization, + ) + + ds = gmt_chol.create_dataset( + ds_name, + data=result.L, compression="gzip", shuffle=True, ) - # Use cached L if available, otherwise compute - L = chol.get(mt, {}).get(mt1, None) - if L is None: - L = _compute_cholesky_factor(M_arr) - gmt_chol.create_dataset( - str(int(mt1)), - data=np.asarray(L, dtype=np.float64), - compression="gzip", - shuffle=True, - ) - - @classmethod - def from_hdf5(cls, filename: str | Path, name: Optional[str] = None) -> "NeutronXSCovariances": - """Read covariances back from an HDF5 file. - - Supports two layouts: - - **standalone** file written by :meth:`to_hdf5` (``mf33`` group - at root) - - **embedded** inside an OpenMC incident-neutron file - (``/covariance/mf33``) - """ - import h5py - - filename = Path(filename) - with h5py.File(filename, "r") as f: - mf33 = cls._find_mf33_group(f) - return cls._read_mf33_group(mf33, name=name or filename.stem) - - @staticmethod - def _find_mf33_group(f) -> "h5py.Group": - """Locate the ``mf33`` group regardless of nesting depth.""" - # Standalone file: /mf33 - if "mf33" in f: - return f["mf33"] - # Embedded: //covariance/mf33 - for key in f: - g = f[key] - if isinstance(g, type(f)) and "covariance" in g: # h5py.Group - cov = g["covariance"] - if "mf33" in cov: - return cov["mf33"] - raise KeyError("No mf33 group found in HDF5 file.") - - @classmethod - def _read_mf33_group(cls, mf33, *, name: str = "") -> "NeutronXSCovariances": - """Build an instance from an ``mf33`` h5py.Group.""" - mat_val = mf33.attrs.get("mat", None) - temp_val = mf33.attrs.get("temperature_k", None) - - ek = np.asarray(mf33["energy_grid_ev"][...], dtype=np.float64) - - reactions: Dict[int, Dict[str, Any]] = {} - greact = mf33["reactions"] - for mt_str, gmt in greact.items(): - mt = int(mt_str) - covs: Dict[int, np.ndarray] = {} - for mt1_str, ds in gmt.items(): - covs[int(mt1_str)] = np.asarray(ds[...], dtype=np.float64) - - reactions[mt] = { - "MAT": int(mat_val) if mat_val is not None else 0, - "MF": 33, - "MT": mt, - "ZA": float(gmt.attrs.get("ZA", 0.0)), - "AWR": float(gmt.attrs.get("AWR", 0.0)), - "COVS": covs, - } - - # Read pre-computed Cholesky factors if present - chol: Optional[Dict[int, Dict[int, np.ndarray]]] = None - if "cholesky" in mf33: - chol = {} - for mt_str, gmt_chol in mf33["cholesky"].items(): - mt = int(mt_str) - chol[mt] = {} - for mt1_str, ds in gmt_chol.items(): - chol[mt][int(mt1_str)] = np.asarray(ds[...], dtype=np.float64) - - return cls( - name=str(name), - energy_grid_ev=ek, - reactions=reactions, - mat=int(mat_val) if mat_val is not None else None, - temperature_k=float(temp_val) if temp_val is not None else None, - cholesky_factors=chol, - ) \ No newline at end of file + # Diagnostic metadata + ds.attrs["effective_rank"] = result.effective_rank + ds.attrs["full_size"] = result.full_size + ds.attrs["method"] = np.bytes_(result.method) + ds.attrs["condition_number"] = result.condition_number + ds.attrs["negative_eig_mass"] = result.negative_eig_mass + ds.attrs["negative_eig_mass_pct"] = result.negative_eig_mass_pct + ds.attrs["regularization"] = result.regularization_applied + ds.attrs["nz_indices"] = result.nonzero_variance_indices \ No newline at end of file diff --git a/src/nuclide.cpp b/src/nuclide.cpp index d379062cfd5..08243d7b16a 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -1162,6 +1162,13 @@ void load_nuclide_covariance(const std::string& name) write_message(6, "Reading covariance data for {} from {}", name, filename); + std::vector e_bounds; + if (object_exists(mf33_group, "energy_grid_ev")) { + read_dataset(mf33_group, "energy_grid_ev", e_bounds); + } else { + warning(fmt::format("No mf33/energy_grid_ev found for {} in {}", name, filename)); + } + // For each reaction, load diagonal block: cholesky/{MT}/{MT} for (auto& rx : nuc->reactions_) { @@ -1173,17 +1180,41 @@ void load_nuclide_covariance(const std::string& name) hid_t mt_group = open_group(chol_group, mt_str.c_str()); if (object_exists(mt_group, mt_str.c_str())) { - // The dataset is 2D (n x n) + // The dataset is 2D (full_size x effective_rank) tensor::Tensor L_matrix; read_dataset(mt_group, mt_str.c_str(), L_matrix); - int n = static_cast(L_matrix.shape(0)); + // Read metadata attributes written by the Python side + hid_t dset_id = open_dataset(mt_group, mt_str.c_str()); + + int full_size = static_cast(L_matrix.shape(0)); + int effective_rank = static_cast(L_matrix.shape(1)); + + // Override with stored attributes when available + if (attribute_exists(dset_id, "full_size")) + read_attribute(dset_id, "full_size", full_size); + if (attribute_exists(dset_id, "effective_rank")) + read_attribute(dset_id, "effective_rank", effective_rank); - // Flatten the 2D matrix into the storage vector + close_dataset(dset_id); + + // Populate the CovData structure Reaction::CovData& cov = rx->cholesky_[rx->mt_]; cov.L.assign(L_matrix.data(), L_matrix.data() + L_matrix.size()); - cov.num_groups = n; - } + cov.full_size = full_size; + cov.effective_rank = effective_rank; + + // Attach boundaries so that perturb_xs can map energies + if (!e_bounds.empty()) { + if (static_cast(e_bounds.size()) == full_size + 1) { + cov.e_bounds = e_bounds; + } else { + warning(fmt::format( + "energy_grid_ev size mismatch for {} MT={}: expected {}, got {}", + name, rx->mt_, full_size + 1, e_bounds.size())); + } + } + } close_group(mt_group); } diff --git a/src/reaction.cpp b/src/reaction.cpp index d387ddeba05..535dafb2378 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -188,27 +188,32 @@ void Reaction::perturb_xs( if (it == cholesky_.end()) return; const auto& cov = it->second; - int ng = cov.num_groups; + int ng = cov.full_size; // number of energy groups (rows of L) + int r = cov.effective_rank; // number of columns of L // 1. Sample z from N(0, I) - vector z(ng); - for (int g = 0; g < ng; ++g) { - double u1 = prn(seed); + vector z(r); + for (int j = 0; j < r; ++j) { + double u1 = std::max(prn(seed), 1e-16); double u2 = prn(seed); - z[g] = std::sqrt(-2.0 * std::log(u1)) * std::cos(2.0 * M_PI * u2); + z[j] = std::sqrt(-2.0 * std::log(u1)) * std::cos(2.0 * M_PI * u2); } // 2. Correlated perturbations: delta = L * z - // L is stored column-major (lower triangular, ng x ng) + // L is stored row-major with shape (ng * r) vector delta(ng, 0.0); for (int i = 0; i < ng; ++i) - for (int j = 0; j <= i; ++j) - delta[i] += cov.L[i * ng + j] * z[j]; + for (int j = 0; j < r; ++j) + delta[i] += cov.L[i * r + j] * z[j]; // 3. Multiplicative factors per group vector factor(ng); - for (int g = 0; g < ng; ++g) - factor[g] = 1.0 + delta[g]; + for (int g = 0; g < ng; ++g){ + double f = 1.0 + delta[g]; + if (f < 0.0) f = 0.0; + if (f > 2.0) f = 2.0; + factor[g] = f; + } // 4. Apply to every temperature's CE cross section for (auto& txs : xs_) { @@ -220,9 +225,15 @@ void Reaction::perturb_xs( // Find which covariance group this energy falls in // e_bounds is sorted ascending, size ng+1 - int g = lower_bound_index( - cov.e_bounds.cbegin(), cov.e_bounds.cend(), E); - g = std::min(g, ng - 1); + if (cov.e_bounds.empty()) continue; + + if (E < cov.e_bounds.front() || E >= cov.e_bounds.back()) { + continue; // outside MG grid => factor 1 + } + + auto itb = std::upper_bound(cov.e_bounds.begin(), cov.e_bounds.end(), E); + int g = static_cast(itb - cov.e_bounds.begin()) - 1; + if (g < 0 || g >= ng) continue; txs.value[k] *= factor[g]; From 479f3cbd9fc5a41b7b12c191e98c18ed5e39cde2 Mon Sep 17 00:00:00 2001 From: Grego01 Date: Mon, 9 Mar 2026 10:45:21 -0400 Subject: [PATCH 06/22] New perturb_xs function --- include/openmc/nuclide.h | 3 ++ include/openmc/reaction.h | 14 ++++--- src/nuclide.cpp | 9 +++++ src/reaction.cpp | 84 ++++++++++++++++++++++++--------------- src/simulation.cpp | 29 ++++++++++++++ 5 files changed, 101 insertions(+), 38 deletions(-) diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index 7a8b2acadd9..1b7bc4fa271 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -129,6 +129,9 @@ class Nuclide { vector> reactions_; //!< Reactions array reaction_index_; //!< Index of each reaction vector index_inelastic_scatter_; + + // Create new set of pertubations cross sections + void reset_derived(); private: void create_derived( diff --git a/include/openmc/reaction.h b/include/openmc/reaction.h index e9b8cb76de7..39bbe8b620e 100644 --- a/include/openmc/reaction.h +++ b/include/openmc/reaction.h @@ -59,12 +59,13 @@ class Reaction { vector value; }; - int mt_; //!< ENDF MT value - double q_value_; //!< Reaction Q value in [eV] - bool scatter_in_cm_; //!< scattering system in center-of-mass? - bool redundant_; //!< redundant reaction? - vector xs_; //!< Cross section at each temperature - vector products_; //!< Reaction products + int mt_; //!< ENDF MT value + double q_value_; //!< Reaction Q value in [eV] + bool scatter_in_cm_; //!< scattering system in center-of-mass? + bool redundant_; //!< redundant reaction? + vector xs_; //!< Cross section at each temperature + vector xs_reference_; //!< Unperturbed cross section at each temperature + vector products_; //!< Reaction products struct CovData { vector L; //!< Flattened (full_size x effective_rank) factor @@ -78,6 +79,7 @@ class Reaction { //! \param[in] energy_grid Nuclide energy grid (the shared CE grid) //! \param[in,out] seed Random number seed void perturb_xs(const vector& energy_grid, uint64_t* seed); + void perturb_xs_covariance(const vector& energy_grid, uint64_t* seed); }; diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 08243d7b16a..d4aadd3f9ff 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -493,6 +493,15 @@ void Nuclide::create_derived( } } +void Nuclide::reset_derived() +{ + xs_.clear(); + fission_rx_.clear(); + n_precursor_ = 0; + create_derived(prompt_photons_.get(), delayed_photons_.get()); + // reaction_index_ and resonant scattering get overwritten, no harm +} + void Nuclide::init_grid() { int neutron = ParticleType::neutron().transport_index(); diff --git a/src/reaction.cpp b/src/reaction.cpp index 535dafb2378..4497b04ca4c 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -13,6 +13,7 @@ #include "openmc/endf.h" #include "openmc/hdf5_interface.h" #include "openmc/random_lcg.h" +#include "openmc/random_dist.h" #include "openmc/search.h" #include "openmc/secondary_uncorrelated.h" #include "openmc/settings.h" @@ -180,65 +181,84 @@ double Reaction::collapse_rate(int64_t i_temp, span energy, return xs_flux_sum; } -void Reaction::perturb_xs( - const vector& energy_grid, uint64_t* seed) +void Reaction::perturb_xs(const vector& energy_grid, uint64_t* seed) { - // Look up Cholesky factors for this reaction's MT + // Perturb cross sections using covariance data auto it = cholesky_.find(mt_); - if (it == cholesky_.end()) return; + if (it != cholesky_.end()){ + perturb_xs_covariance(energy_grid, seed); + return; + } + + // Perturb manually cross sections (no covariance data) + + double rel_std = 0.00; //TO DO NEW SETTINGS TO ADD + double factor = 1.0; + + // Normal distribution + double z = normal_variate(0.0, 1.0, seed); + factor = 1.0 + rel_std * z; - const auto& cov = it->second; - int ng = cov.full_size; // number of energy groups (rows of L) - int r = cov.effective_rank; // number of columns of L + factor = std::max(0.0, std::min(2.0, factor)); - // 1. Sample z from N(0, I) + for (auto& txs : xs_) + { + for (auto& v : txs.value) + { + v *= factor; + if (v < 0.0) v = 0.0; + } + } +} + +// Function to sample cross sections using covariance data + +void Reaction::perturb_xs_covariance(const vector& energy_grid, uint64_t* seed) +{ + const auto& cov = cholesky_[mt_]; + int ng = cov.full_size; + int r = cov.effective_rank; + + // First step is to sample from normal distribution vector z(r); - for (int j = 0; j < r; ++j) { - double u1 = std::max(prn(seed), 1e-16); - double u2 = prn(seed); - z[j] = std::sqrt(-2.0 * std::log(u1)) * std::cos(2.0 * M_PI * u2); + for (int j = 0; j < r; ++j) + { + z[j] = normal_variate(0.0, 1.0, seed); } - // 2. Correlated perturbations: delta = L * z - // L is stored row-major with shape (ng * r) + // Second step is to compute correlated pertubations vector delta(ng, 0.0); for (int i = 0; i < ng; ++i) for (int j = 0; j < r; ++j) - delta[i] += cov.L[i * r + j] * z[j]; + delta[i] += cov.L[i * r + j] * z[j]; - // 3. Multiplicative factors per group + // Third step is to obtain multiplicative factors per group vector factor(ng); - for (int g = 0; g < ng; ++g){ + for (int g = 0; g < ng; ++g) + { double f = 1.0 + delta[g]; - if (f < 0.0) f = 0.0; - if (f > 2.0) f = 2.0; - factor[g] = f; + factor[g] = std::max(0.0, std::min(2.0, f)); } - // 4. Apply to every temperature's CE cross section - for (auto& txs : xs_) { + // Apply the perturbation to every temperature CE XS + for (auto& txs : xs_) + { int n_points = txs.value.size(); - for (int k = 0; k < n_points; ++k) { - // Map CE grid index to energy + for (int k = 0; k < n_points; ++k) + { int i_grid = k + txs.threshold; double E = energy_grid[i_grid]; - // Find which covariance group this energy falls in - // e_bounds is sorted ascending, size ng+1 if (cov.e_bounds.empty()) continue; - - if (E < cov.e_bounds.front() || E >= cov.e_bounds.back()) { - continue; // outside MG grid => factor 1 - } + if (E < cov.e_bounds.front() || E >= cov.e_bounds.back()) continue; auto itb = std::upper_bound(cov.e_bounds.begin(), cov.e_bounds.end(), E); int g = static_cast(itb - cov.e_bounds.begin()) - 1; if (g < 0 || g >= ng) continue; txs.value[k] *= factor[g]; + if (txs.value[k] < 0.0) txs.value[k] = 0.0; // to avoid negative cross sections - // Keep cross sections non-negative - if (txs.value[k] < 0.0) txs.value[k] = 0.0; } } } diff --git a/src/simulation.cpp b/src/simulation.cpp index 4fad196a604..800e3f478b2 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -162,6 +162,14 @@ int openmc_simulation_init() openmc_weight_windows_import(settings::weight_windows_file.c_str()); } + if (settings::load_covariance) { + for (auto& nuc : data::nuclides) { + for (auto& rx : nuc->reactions_) { + rx->xs_reference_ = rx->xs_; + } + } + } + // Set flag indicating initialization is done simulation::initialized = true; return 0; @@ -398,6 +406,27 @@ void initialize_batch() // Add user tallies to active tallies list setup_active_tallies(); + + //=== UQ perturbation ===// + if (settings::load_covariance) { + uint64_t seed = static_cast(simulation::current_batch); + + for (auto& nuc : data::nuclides) { + for (auto& rx : nuc->reactions_) { + // 1. Restore originals + rx->xs_ = rx->xs_reference_; + + // 2. Apply fresh perturbation + rx->perturb_xs(nuc->grid_[0].energy, &seed); + } + + // 3. Create perturbed derived tables (total, absorption, fission, etc.) + // valid only for redundant reaction but not its components like for + // example MT=4. In this case, ERRORR provides covaraince for MT=4 + // but not individually for MT=51 through MT=91! + nuc->reset_derived(); + } + } } void finalize_batch() From b3a62878eeb55c45c80c75a715fe8d7335ea73b4 Mon Sep 17 00:00:00 2001 From: Grego01 Date: Mon, 16 Mar 2026 11:49:53 -0400 Subject: [PATCH 07/22] Sampling on the fly working, more tests needed --- src/simulation.cpp | 54 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/src/simulation.cpp b/src/simulation.cpp index 800e3f478b2..ab1cf7bbfd8 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -412,6 +412,27 @@ void initialize_batch() uint64_t seed = static_cast(simulation::current_batch); for (auto& nuc : data::nuclides) { + + // ------ DIAGNOSTIC: snapshot a few total-XS values BEFORE ------ + // Pick a grid point in the middle of the energy range for temperature 0 + /*int diag_idx = -1; + double total_before = -1.0; + double elastic_rx_before = -1.0; + if (!nuc->xs_.empty() && !nuc->grid_.empty()) { + int n_grid = nuc->grid_[0].energy.size(); + diag_idx = n_grid / 2; // middle of the grid + // Nuclide::xs_[0] is the derived tensor for temperature 0 + // Column 0 = XS_TOTAL + total_before = nuc->xs_[0](diag_idx, 0); + // Also grab the elastic reaction XS at that grid point + if (!nuc->reactions_.empty()) { + auto& el = nuc->reactions_[0]; // MT=2 is typically index 0 + int thr = el->xs_[0].threshold; + if (diag_idx >= thr && diag_idx - thr < (int)el->xs_[0].value.size()) + elastic_rx_before = el->xs_[0].value[diag_idx - thr]; + } + }*/ + for (auto& rx : nuc->reactions_) { // 1. Restore originals rx->xs_ = rx->xs_reference_; @@ -420,11 +441,44 @@ void initialize_batch() rx->perturb_xs(nuc->grid_[0].energy, &seed); } + // ------ DIAGNOSTIC: check elastic reaction XS after perturbation ------ + /*double elastic_rx_after = -1.0; + if (diag_idx >= 0 && !nuc->reactions_.empty()) { + auto& el = nuc->reactions_[0]; + int thr = el->xs_[0].threshold; + if (diag_idx >= thr && diag_idx - thr < (int)el->xs_[0].value.size()) + elastic_rx_after = el->xs_[0].value[diag_idx - thr]; + }*/ + // 3. Create perturbed derived tables (total, absorption, fission, etc.) // valid only for redundant reaction but not its components like for // example MT=4. In this case, ERRORR provides covaraince for MT=4 // but not individually for MT=51 through MT=91! nuc->reset_derived(); + + /*double total_after = -1.0; + if (!nuc->xs_.empty() && diag_idx >= 0) { + total_after = nuc->xs_[0](diag_idx, 0); + } + + if (diag_idx >= 0) { + double E_diag = nuc->grid_[0].energy[diag_idx]; + fmt::print("[derived-check] {} batch={} E={:.6e} eV\n", + nuc->name_, simulation::current_batch, E_diag); + fmt::print(" elastic_rx: before={:.10e} after={:.10e} ratio={:.8f}\n", + elastic_rx_before, elastic_rx_after, + (elastic_rx_before > 0 ? elastic_rx_after / elastic_rx_before : 0.0)); + fmt::print(" derived_total: before={:.10e} after={:.10e} ratio={:.8f}\n", + total_before, total_after, + (total_before > 0 ? total_after / total_before : 0.0)); + + if (total_before == total_after) { + fmt::print(" *** DERIVED TOTAL IS UNCHANGED — BUG CONFIRMED ***\n"); + } else { + fmt::print(" derived total changed by {:.6f}%\n", + 100.0 * (total_after - total_before) / total_before); + } + }*/ } } } From ad02a9b5ad00c91ac3c86072a1625bf396d991f8 Mon Sep 17 00:00:00 2001 From: Gregoire Biot Date: Fri, 27 Mar 2026 14:04:32 -0400 Subject: [PATCH 08/22] Delete openmc/data/mf33_enjoy - Copy.py:Zone.Identifier --- openmc/data/mf33_enjoy - Copy.py:Zone.Identifier | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 openmc/data/mf33_enjoy - Copy.py:Zone.Identifier diff --git a/openmc/data/mf33_enjoy - Copy.py:Zone.Identifier b/openmc/data/mf33_enjoy - Copy.py:Zone.Identifier deleted file mode 100644 index e69de29bb2d..00000000000 From 27eff53700b17e640f116558aba148884e95ceae Mon Sep 17 00:00:00 2001 From: Gregoire Biot Date: Fri, 27 Mar 2026 14:04:56 -0400 Subject: [PATCH 09/22] Delete openmc/data/xs_covariance_njoy - Copy.py:Zone.Identifier --- openmc/data/xs_covariance_njoy - Copy.py:Zone.Identifier | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 openmc/data/xs_covariance_njoy - Copy.py:Zone.Identifier diff --git a/openmc/data/xs_covariance_njoy - Copy.py:Zone.Identifier b/openmc/data/xs_covariance_njoy - Copy.py:Zone.Identifier deleted file mode 100644 index e69de29bb2d..00000000000 From c5a9e606e257750aa26190bab6f0ee5f17c24f80 Mon Sep 17 00:00:00 2001 From: Grego01 Date: Mon, 30 Mar 2026 09:41:01 -0400 Subject: [PATCH 10/22] cleaning up --- src/cross_sections.cpp | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/cross_sections.cpp b/src/cross_sections.cpp index d8c89e35871..ec9da5b8abc 100644 --- a/src/cross_sections.cpp +++ b/src/cross_sections.cpp @@ -214,13 +214,6 @@ void read_ce_cross_sections(const vector>& nuc_temps, throw std::runtime_error {openmc_err_msg}; already_read.insert(name); - - /*uint64_t xs_perturbation_seed = 42; - - auto& nuc = *data::nuclides.back(); - for (auto& rx : nuc.reactions_) { - rx->perturb_xs(nuc.grid_[0].energy, &xs_perturbation_seed); - }*/ } } From 446b08c5d23b55e9e2ceee7ce761cef0345e342c Mon Sep 17 00:00:00 2001 From: Grego01 Date: Mon, 30 Mar 2026 10:18:10 -0400 Subject: [PATCH 11/22] cleaning up --- src/simulation.cpp | 56 ---------------------------------------------- 1 file changed, 56 deletions(-) diff --git a/src/simulation.cpp b/src/simulation.cpp index ab1cf7bbfd8..9ecfe6b2ce5 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -413,26 +413,6 @@ void initialize_batch() for (auto& nuc : data::nuclides) { - // ------ DIAGNOSTIC: snapshot a few total-XS values BEFORE ------ - // Pick a grid point in the middle of the energy range for temperature 0 - /*int diag_idx = -1; - double total_before = -1.0; - double elastic_rx_before = -1.0; - if (!nuc->xs_.empty() && !nuc->grid_.empty()) { - int n_grid = nuc->grid_[0].energy.size(); - diag_idx = n_grid / 2; // middle of the grid - // Nuclide::xs_[0] is the derived tensor for temperature 0 - // Column 0 = XS_TOTAL - total_before = nuc->xs_[0](diag_idx, 0); - // Also grab the elastic reaction XS at that grid point - if (!nuc->reactions_.empty()) { - auto& el = nuc->reactions_[0]; // MT=2 is typically index 0 - int thr = el->xs_[0].threshold; - if (diag_idx >= thr && diag_idx - thr < (int)el->xs_[0].value.size()) - elastic_rx_before = el->xs_[0].value[diag_idx - thr]; - } - }*/ - for (auto& rx : nuc->reactions_) { // 1. Restore originals rx->xs_ = rx->xs_reference_; @@ -441,44 +421,8 @@ void initialize_batch() rx->perturb_xs(nuc->grid_[0].energy, &seed); } - // ------ DIAGNOSTIC: check elastic reaction XS after perturbation ------ - /*double elastic_rx_after = -1.0; - if (diag_idx >= 0 && !nuc->reactions_.empty()) { - auto& el = nuc->reactions_[0]; - int thr = el->xs_[0].threshold; - if (diag_idx >= thr && diag_idx - thr < (int)el->xs_[0].value.size()) - elastic_rx_after = el->xs_[0].value[diag_idx - thr]; - }*/ - - // 3. Create perturbed derived tables (total, absorption, fission, etc.) - // valid only for redundant reaction but not its components like for - // example MT=4. In this case, ERRORR provides covaraince for MT=4 - // but not individually for MT=51 through MT=91! nuc->reset_derived(); - /*double total_after = -1.0; - if (!nuc->xs_.empty() && diag_idx >= 0) { - total_after = nuc->xs_[0](diag_idx, 0); - } - - if (diag_idx >= 0) { - double E_diag = nuc->grid_[0].energy[diag_idx]; - fmt::print("[derived-check] {} batch={} E={:.6e} eV\n", - nuc->name_, simulation::current_batch, E_diag); - fmt::print(" elastic_rx: before={:.10e} after={:.10e} ratio={:.8f}\n", - elastic_rx_before, elastic_rx_after, - (elastic_rx_before > 0 ? elastic_rx_after / elastic_rx_before : 0.0)); - fmt::print(" derived_total: before={:.10e} after={:.10e} ratio={:.8f}\n", - total_before, total_after, - (total_before > 0 ? total_after / total_before : 0.0)); - - if (total_before == total_after) { - fmt::print(" *** DERIVED TOTAL IS UNCHANGED — BUG CONFIRMED ***\n"); - } else { - fmt::print(" derived total changed by {:.6f}%\n", - 100.0 * (total_after - total_before) / total_before); - } - }*/ } } } From ec0cf79e54a5d159379630e68f22555507e25844 Mon Sep 17 00:00:00 2001 From: Grego01 Date: Mon, 30 Mar 2026 11:32:57 -0400 Subject: [PATCH 12/22] clean up --- openmc/data/neutron.py | 127 +++-------------------------------------- 1 file changed, 9 insertions(+), 118 deletions(-) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index e32f965ac21..f0f4100e5cc 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -231,14 +231,6 @@ def urr(self, urr): @property def mg_covariance(self): - """Multigroup cross-section covariance data (your MF=33-derived product). - - Expected schema (dict-like): - - energy_bounds: (G+1,) float - - mts: (N,) int - - matrix: (...) float - - representation: optional str, e.g. "relative" - """ return self._mg_covariance @mg_covariance.setter @@ -246,10 +238,10 @@ def mg_covariance(self, cov): if cov is None: self._mg_covariance = None return - cv.check_type('mg_covariance', cov, Mapping) - for k in ('energy_bounds', 'mts', 'matrix'): - if k not in cov: - raise KeyError(f"mg_covariance missing required key '{k}'") + if not hasattr(cov, 'write_mf33_group'): + raise TypeError( + "mg_covariance must provide write_mf33_group " + "(e.g. NeutronXSCovariances)") self._mg_covariance = cov # Optional alias if you already use `data.covariance = ...` somewhere @@ -470,71 +462,7 @@ def export_to_hdf5(self, path, mode='a', libver='earliest'): # Prefer the shared writer if available (NeutronXSCovariances) if hasattr(cov, "write_mf33_group"): cov.write_mf33_group(cov_root) - else: - # Fallback: write directly from duck-typed attributes - if "mf33" in cov_root: - del cov_root["mf33"] - mf33 = cov_root.create_group("mf33") - - mf33.attrs["format"] = np.bytes_("openmc.mf33.v1") - mf33.attrs["source"] = np.bytes_("njoy errorr module") - mf33.attrs["relative"] = 1 - if getattr(cov, "mat", None) is not None: - mf33.attrs["mat"] = int(cov.mat) - if getattr(cov, "temperature_k", None) is not None: - mf33.attrs["temperature_k"] = float(cov.temperature_k) - - mf33.create_dataset( - "energy_grid_ev", - data=np.asarray(cov.energy_grid_ev, dtype=np.float64), - ) - - greact = mf33.create_group("reactions") - gchol = mf33.create_group("cholesky") - - # Grab cached Cholesky factors if available - chol = getattr(cov, "cholesky_factors", None) or {} - - for mt, sec in cov.reactions.items(): - gmt = greact.create_group(str(int(mt))) - gmt.attrs["ZA"] = float(sec.get("ZA", 0.0)) - gmt.attrs["AWR"] = float(sec.get("AWR", 0.0)) - - gmt_chol = gchol.create_group(str(int(mt))) - gmt_chol.attrs["ZA"] = float(sec.get("ZA", 0.0)) - gmt_chol.attrs["AWR"] = float(sec.get("AWR", 0.0)) - - covs = sec.get("COVS", {}) - for mt1, M in covs.items(): - M_arr = np.asarray(M, dtype=np.float64) - gmt.create_dataset( - str(int(mt1)), - data=M_arr, - compression="gzip", - shuffle=True, - ) - - # Use cached L or compute on the fly - L = chol.get(mt, {}).get(mt1, None) - if L is None: - import scipy.linalg as _la - try: - L = _la.cholesky(M_arr, lower=True) - except _la.LinAlgError: - evals, V = _la.eigh(M_arr) - max_e = evals.max() if evals.max() > 0 else 1.0 - evals[evals / max_e < 1e-10] = 0.0 - evals[evals < 0.0] = 0.0 - A_tmp = V * np.sqrt(evals)[np.newaxis, :] - _, R_tmp = _la.qr(A_tmp.T, mode='economic') - L = R_tmp.T - gmt_chol.create_dataset( - str(int(mt1)), - data=np.asarray(L, dtype=np.float64), - compression="gzip", - shuffle=True, - ) - + @classmethod def from_hdf5(cls, group_or_filename): """Generate continuous-energy neutron interaction data from HDF5 group @@ -612,47 +540,10 @@ def from_hdf5(cls, group_or_filename): # Read covariance data (per-reaction MF=33 schema) if 'covariance' in group and 'mf33' in group['covariance']: - try: - from .xs_covariance import NeutronXSCovariances - data._mg_covariance = NeutronXSCovariances._read_mf33_group( - group['covariance']['mf33'], name=name, - ) - except ImportError: - # Fallback: read into a lightweight SimpleNamespace - import types - mf33 = group['covariance']['mf33'] - ek = np.asarray(mf33['energy_grid_ev'][...], dtype=np.float64) - reactions = {} - for mt_str, gmt in mf33['reactions'].items(): - mt = int(mt_str) - covs = {} - for mt1_str, ds in gmt.items(): - covs[int(mt1_str)] = np.asarray(ds[...], dtype=np.float64) - reactions[mt] = { - 'ZA': float(gmt.attrs.get('ZA', 0.0)), - 'AWR': float(gmt.attrs.get('AWR', 0.0)), - 'COVS': covs, - } - - # Read pre-computed Cholesky factors if present - chol = None - if 'cholesky' in mf33: - chol = {} - for mt_str, gmt_chol in mf33['cholesky'].items(): - mt = int(mt_str) - chol[mt] = {} - for mt1_str, ds in gmt_chol.items(): - chol[mt][int(mt1_str)] = np.asarray(ds[...], dtype=np.float64) - - obj = types.SimpleNamespace( - energy_grid_ev=ek, - reactions=reactions, - mat=int(mf33.attrs['mat']) if 'mat' in mf33.attrs else None, - temperature_k=float(mf33.attrs['temperature_k']) if 'temperature_k' in mf33.attrs else None, - cholesky_factors=chol, - ) - data._mg_covariance = obj - + from .xs_covariance_njoy import NeutronXSCovariances + data._mg_covariance = NeutronXSCovariances._read_mf33_group( + group['covariance']['mf33'], name=name,) + return data @classmethod From 92a3177130f30af95248677eaf650bda82657b7b Mon Sep 17 00:00:00 2001 From: Grego01 Date: Tue, 31 Mar 2026 10:03:36 -0400 Subject: [PATCH 13/22] delete modifications on headers --- include/openmc/nuclide.h | 3 --- include/openmc/reaction.h | 14 -------------- include/openmc/settings.h | 1 - 3 files changed, 18 deletions(-) diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index 1b7bc4fa271..7a8b2acadd9 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -129,9 +129,6 @@ class Nuclide { vector> reactions_; //!< Reactions array reaction_index_; //!< Index of each reaction vector index_inelastic_scatter_; - - // Create new set of pertubations cross sections - void reset_derived(); private: void create_derived( diff --git a/include/openmc/reaction.h b/include/openmc/reaction.h index 39bbe8b620e..73127935f49 100644 --- a/include/openmc/reaction.h +++ b/include/openmc/reaction.h @@ -67,20 +67,6 @@ class Reaction { vector xs_reference_; //!< Unperturbed cross section at each temperature vector products_; //!< Reaction products - struct CovData { - vector L; //!< Flattened (full_size x effective_rank) factor - vector e_bounds; //!< Energy group boundaries [eV], length = full_size+1 - int full_size {0}; //!< Original matrix dimension G (number of groups) - int effective_rank {0}; //!< Number of columns retained (r <= G) - }; - std::unordered_map cholesky_; // MT → Cholesky factor - - //! Apply multigroup perturbation factors to CE cross sections - //! \param[in] energy_grid Nuclide energy grid (the shared CE grid) - //! \param[in,out] seed Random number seed - void perturb_xs(const vector& energy_grid, uint64_t* seed); - void perturb_xs_covariance(const vector& energy_grid, uint64_t* seed); - }; //============================================================================== diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 182e6a38e2d..0914a0958b8 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -202,7 +202,6 @@ extern int trigger_batch_interval; //!< Batch interval for triggers extern "C" int verbosity; //!< How verbose to make output extern double weight_cutoff; //!< Weight cutoff for Russian roulette extern double weight_survive; //!< Survival weight after Russian roulette -extern bool load_covariance; //!< Load covariance data for sampling } // namespace settings From c9b7af89690c83c9e290dc6bfb72049058a6b482 Mon Sep 17 00:00:00 2001 From: Grego01 Date: Tue, 31 Mar 2026 10:06:32 -0400 Subject: [PATCH 14/22] delete modifications on source files --- src/nuclide.cpp | 135 --------------------------------------------- src/reaction.cpp | 82 --------------------------- src/settings.cpp | 6 -- src/simulation.cpp | 27 --------- 4 files changed, 250 deletions(-) diff --git a/src/nuclide.cpp b/src/nuclide.cpp index d4aadd3f9ff..b7736c83bbc 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -493,15 +493,6 @@ void Nuclide::create_derived( } } -void Nuclide::reset_derived() -{ - xs_.clear(); - fission_rx_.clear(); - n_precursor_ = 0; - create_derived(prompt_photons_.get(), delayed_photons_.get()); - // reaction_index_ and resonant scattering get overwritten, no harm -} - void Nuclide::init_grid() { int neutron = ParticleType::neutron().transport_index(); @@ -1110,132 +1101,6 @@ extern "C" size_t nuclides_size() return data::nuclides.size(); } -//============================================================================== -// Covariance loading -//============================================================================== -void load_nuclide_covariance(const std::string& name) -{ - - // Make sure the nuclide is loaded - auto it_nuc = data::nuclide_map.find(name); - if (it_nuc == data::nuclide_map.end()) - return; - - int i_nuclide = it_nuc->second; - auto& nuc = data::nuclides[i_nuclide]; - - // Find the same neutron library file used for cross sections - LibraryKey key {Library::Type::neutron, name}; - const auto& it_lib = data::library_map.find(key); - if (it_lib == data::library_map.end()) - return; - - int idx = it_lib->second; - const auto& filename = data::libraries[idx].path_; - - // Open the standard neutron XS library file - hid_t file_id = file_open(filename, 'r'); - check_data_version(file_id); - - // Open nuclide group: "/Fe56" - if (!object_exists(file_id, name.c_str())) { - file_close(file_id); - return; - } - hid_t nuc_group = open_group(file_id, name.c_str()); - - // Walk: covariance/mf33/cholesky - if (!object_exists(nuc_group, "covariance")) { - close_group(nuc_group); - file_close(file_id); - return; - } - hid_t cov_group = open_group(nuc_group, "covariance"); - - if (!object_exists(cov_group, "mf33")) { - close_group(cov_group); - close_group(nuc_group); - file_close(file_id); - return; - } - hid_t mf33_group = open_group(cov_group, "mf33"); - - if (!object_exists(mf33_group, "cholesky")) { - close_group(mf33_group); - close_group(cov_group); - close_group(nuc_group); - file_close(file_id); - return; - } - hid_t chol_group = open_group(mf33_group, "cholesky"); - - write_message(6, "Reading covariance data for {} from {}", name, filename); - - std::vector e_bounds; - if (object_exists(mf33_group, "energy_grid_ev")) { - read_dataset(mf33_group, "energy_grid_ev", e_bounds); - } else { - warning(fmt::format("No mf33/energy_grid_ev found for {} in {}", name, filename)); - } - - // For each reaction, load diagonal block: cholesky/{MT}/{MT} - for (auto& rx : nuc->reactions_) { - - std::string mt_str = std::to_string(rx->mt_); - - if (!object_exists(chol_group, mt_str.c_str())) - continue; - - hid_t mt_group = open_group(chol_group, mt_str.c_str()); - - if (object_exists(mt_group, mt_str.c_str())) { - // The dataset is 2D (full_size x effective_rank) - tensor::Tensor L_matrix; - read_dataset(mt_group, mt_str.c_str(), L_matrix); - - // Read metadata attributes written by the Python side - hid_t dset_id = open_dataset(mt_group, mt_str.c_str()); - - int full_size = static_cast(L_matrix.shape(0)); - int effective_rank = static_cast(L_matrix.shape(1)); - - // Override with stored attributes when available - if (attribute_exists(dset_id, "full_size")) - read_attribute(dset_id, "full_size", full_size); - if (attribute_exists(dset_id, "effective_rank")) - read_attribute(dset_id, "effective_rank", effective_rank); - - close_dataset(dset_id); - - // Populate the CovData structure - Reaction::CovData& cov = rx->cholesky_[rx->mt_]; - cov.L.assign(L_matrix.data(), L_matrix.data() + L_matrix.size()); - cov.full_size = full_size; - cov.effective_rank = effective_rank; - - // Attach boundaries so that perturb_xs can map energies - if (!e_bounds.empty()) { - if (static_cast(e_bounds.size()) == full_size + 1) { - cov.e_bounds = e_bounds; - } else { - warning(fmt::format( - "energy_grid_ev size mismatch for {} MT={}: expected {}, got {}", - name, rx->mt_, full_size + 1, e_bounds.size())); - } - } - } - - close_group(mt_group); - } - - // Close groups/files - close_group(chol_group); - close_group(mf33_group); - close_group(cov_group); - close_group(nuc_group); - file_close(file_id); -} - //============================================================================== // C API //============================================================================== diff --git a/src/reaction.cpp b/src/reaction.cpp index 4497b04ca4c..c3f2daf67f2 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -181,88 +181,6 @@ double Reaction::collapse_rate(int64_t i_temp, span energy, return xs_flux_sum; } -void Reaction::perturb_xs(const vector& energy_grid, uint64_t* seed) -{ - // Perturb cross sections using covariance data - auto it = cholesky_.find(mt_); - if (it != cholesky_.end()){ - perturb_xs_covariance(energy_grid, seed); - return; - } - - // Perturb manually cross sections (no covariance data) - - double rel_std = 0.00; //TO DO NEW SETTINGS TO ADD - double factor = 1.0; - - // Normal distribution - double z = normal_variate(0.0, 1.0, seed); - factor = 1.0 + rel_std * z; - - factor = std::max(0.0, std::min(2.0, factor)); - - for (auto& txs : xs_) - { - for (auto& v : txs.value) - { - v *= factor; - if (v < 0.0) v = 0.0; - } - } -} - -// Function to sample cross sections using covariance data - -void Reaction::perturb_xs_covariance(const vector& energy_grid, uint64_t* seed) -{ - const auto& cov = cholesky_[mt_]; - int ng = cov.full_size; - int r = cov.effective_rank; - - // First step is to sample from normal distribution - vector z(r); - for (int j = 0; j < r; ++j) - { - z[j] = normal_variate(0.0, 1.0, seed); - } - - // Second step is to compute correlated pertubations - vector delta(ng, 0.0); - for (int i = 0; i < ng; ++i) - for (int j = 0; j < r; ++j) - delta[i] += cov.L[i * r + j] * z[j]; - - // Third step is to obtain multiplicative factors per group - vector factor(ng); - for (int g = 0; g < ng; ++g) - { - double f = 1.0 + delta[g]; - factor[g] = std::max(0.0, std::min(2.0, f)); - } - - // Apply the perturbation to every temperature CE XS - for (auto& txs : xs_) - { - int n_points = txs.value.size(); - for (int k = 0; k < n_points; ++k) - { - int i_grid = k + txs.threshold; - double E = energy_grid[i_grid]; - - if (cov.e_bounds.empty()) continue; - if (E < cov.e_bounds.front() || E >= cov.e_bounds.back()) continue; - - auto itb = std::upper_bound(cov.e_bounds.begin(), cov.e_bounds.end(), E); - int g = static_cast(itb - cov.e_bounds.begin()) - 1; - if (g < 0 || g >= ng) continue; - - txs.value[k] *= factor[g]; - if (txs.value[k] < 0.0) txs.value[k] = 0.0; // to avoid negative cross sections - - } - } -} - //============================================================================== // Non-member functions //============================================================================== diff --git a/src/settings.cpp b/src/settings.cpp index bb50f70f509..a1942c93cae 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -152,7 +152,6 @@ int trigger_batch_interval {1}; int verbosity {-1}; double weight_cutoff {0.25}; double weight_survive {1.0}; -bool load_covariance {false}; } // namespace settings @@ -1304,11 +1303,6 @@ void read_settings_xml(pugi::xml_node root) get_node_value_bool(root, "use_decay_photons"); } - if (check_for_node(root, "load_covariance")) { - settings::load_covariance = - get_node_value_bool(root, "load_covariance"); - } - } diff --git a/src/simulation.cpp b/src/simulation.cpp index 9ecfe6b2ce5..90b96b5c418 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -160,14 +160,6 @@ int openmc_simulation_init() // load weight windows from file if (!settings::weight_windows_file.empty()) { openmc_weight_windows_import(settings::weight_windows_file.c_str()); - } - - if (settings::load_covariance) { - for (auto& nuc : data::nuclides) { - for (auto& rx : nuc->reactions_) { - rx->xs_reference_ = rx->xs_; - } - } } // Set flag indicating initialization is done @@ -406,25 +398,6 @@ void initialize_batch() // Add user tallies to active tallies list setup_active_tallies(); - - //=== UQ perturbation ===// - if (settings::load_covariance) { - uint64_t seed = static_cast(simulation::current_batch); - - for (auto& nuc : data::nuclides) { - - for (auto& rx : nuc->reactions_) { - // 1. Restore originals - rx->xs_ = rx->xs_reference_; - - // 2. Apply fresh perturbation - rx->perturb_xs(nuc->grid_[0].energy, &seed); - } - - nuc->reset_derived(); - - } - } } void finalize_batch() From ed90799600c9261943bbcc4a3017b707b32c1356 Mon Sep 17 00:00:00 2001 From: Grego01 Date: Tue, 31 Mar 2026 10:09:45 -0400 Subject: [PATCH 15/22] delete all modifications on source files --- src/nuclide.cpp | 4 ---- src/reaction.cpp | 2 -- src/settings.cpp | 2 -- 3 files changed, 8 deletions(-) diff --git a/src/nuclide.cpp b/src/nuclide.cpp index b7736c83bbc..17d6e952c39 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -1139,10 +1139,6 @@ extern "C" int openmc_load_nuclide(const char* name, const double* temps, int n) if (settings::temperature_multipole) read_multipole_data(i_nuclide); - // Read covariance data (Cholesky factors) if requested - if (settings::load_covariance) - load_nuclide_covariance(name); - // Read elemental data, if necessary if (settings::photon_transport) { auto element = to_element(name); diff --git a/src/reaction.cpp b/src/reaction.cpp index c3f2daf67f2..c02e0cc407e 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -4,7 +4,6 @@ #include #include #include // for move -#include #include @@ -13,7 +12,6 @@ #include "openmc/endf.h" #include "openmc/hdf5_interface.h" #include "openmc/random_lcg.h" -#include "openmc/random_dist.h" #include "openmc/search.h" #include "openmc/secondary_uncorrelated.h" #include "openmc/settings.h" diff --git a/src/settings.cpp b/src/settings.cpp index a1942c93cae..ab9f9a5aafa 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -1302,8 +1302,6 @@ void read_settings_xml(pugi::xml_node root) settings::use_decay_photons = get_node_value_bool(root, "use_decay_photons"); } - - } void free_memory_settings() From d56f436703bd2cfe3b83821b5d9016478353cd86 Mon Sep 17 00:00:00 2001 From: Grego01 Date: Tue, 31 Mar 2026 10:12:27 -0400 Subject: [PATCH 16/22] delete all modifications on C++ side --- include/openmc/reaction.h | 2 -- src/simulation.cpp | 1 - 2 files changed, 3 deletions(-) diff --git a/include/openmc/reaction.h b/include/openmc/reaction.h index 73127935f49..d63b8e60d7d 100644 --- a/include/openmc/reaction.h +++ b/include/openmc/reaction.h @@ -5,7 +5,6 @@ #define OPENMC_REACTION_H #include -#include #include "hdf5.h" @@ -64,7 +63,6 @@ class Reaction { bool scatter_in_cm_; //!< scattering system in center-of-mass? bool redundant_; //!< redundant reaction? vector xs_; //!< Cross section at each temperature - vector xs_reference_; //!< Unperturbed cross section at each temperature vector products_; //!< Reaction products }; diff --git a/src/simulation.cpp b/src/simulation.cpp index 90b96b5c418..d44cb03b6a1 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -161,7 +161,6 @@ int openmc_simulation_init() if (!settings::weight_windows_file.empty()) { openmc_weight_windows_import(settings::weight_windows_file.c_str()); } - // Set flag indicating initialization is done simulation::initialized = true; return 0; From 7cd1d2b9b88ddd31d0399071bdda26cc7bfa75fd Mon Sep 17 00:00:00 2001 From: Grego01 Date: Tue, 31 Mar 2026 10:14:42 -0400 Subject: [PATCH 17/22] no boolean flag needed --- openmc/settings.py | 16 ---------------- 1 file changed, 16 deletions(-) diff --git a/openmc/settings.py b/openmc/settings.py index 224126419a5..6919afca4cc 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -485,7 +485,6 @@ def __init__(self, **kwargs): self._max_tracks = None self._max_secondaries = None self._use_decay_photons = None - self._load_covariance = None self._random_ray = {} @@ -1462,15 +1461,6 @@ def free_gas_threshold(self, free_gas_threshold: float | None): cv.check_greater_than('free gas threshold', free_gas_threshold, 0.0) self._free_gas_threshold = free_gas_threshold - @property - def load_covariance(self): - return self._load_covariance - - @load_covariance.setter - def load_covariance(self, value): - cv.check_type('load_covariance', value, bool) - self._load_covariance = value - def _create_run_mode_subelement(self, root): elem = ET.SubElement(root, "run_mode") elem.text = self._run_mode.value @@ -1835,11 +1825,6 @@ def _create_use_decay_photons_subelement(self, root): if self._use_decay_photons is not None: element = ET.SubElement(root, "use_decay_photons") element.text = str(self._use_decay_photons).lower() - - def _create_load_covariance_subelement(self, root): - if self._load_covariance is not None: - element = ET.SubElement(root, "load_covariance") - element.text = str(self._load_covariance).lower() def _create_resonance_scattering_subelement(self, root): res = self.resonance_scattering @@ -2591,7 +2576,6 @@ def to_xml_element(self, mesh_memo=None): self._create_use_decay_photons_subelement(element) self._create_source_rejection_fraction_subelement(element) self._create_free_gas_threshold_subelement(element) - self._create_load_covariance_subelement(element) # Clean the indentation in the file to be user-readable clean_indentation(element) From ad6f17d8d48571a67b4f37ba64a9025862469964 Mon Sep 17 00:00:00 2001 From: Gregoire Biot Date: Tue, 31 Mar 2026 10:17:53 -0400 Subject: [PATCH 18/22] Missing space --- src/simulation.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/simulation.cpp b/src/simulation.cpp index d44cb03b6a1..1023b3d321e 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -161,6 +161,7 @@ int openmc_simulation_init() if (!settings::weight_windows_file.empty()) { openmc_weight_windows_import(settings::weight_windows_file.c_str()); } + // Set flag indicating initialization is done simulation::initialized = true; return 0; From e90dc5ee115edcc5d9d515c147e6abe9854d246a Mon Sep 17 00:00:00 2001 From: Gregoire Biot Date: Tue, 31 Mar 2026 10:19:16 -0400 Subject: [PATCH 19/22] Fix formatting of reaction.h file --- include/openmc/reaction.h | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/include/openmc/reaction.h b/include/openmc/reaction.h index d63b8e60d7d..7389c69c1a6 100644 --- a/include/openmc/reaction.h +++ b/include/openmc/reaction.h @@ -58,12 +58,12 @@ class Reaction { vector value; }; - int mt_; //!< ENDF MT value - double q_value_; //!< Reaction Q value in [eV] - bool scatter_in_cm_; //!< scattering system in center-of-mass? - bool redundant_; //!< redundant reaction? - vector xs_; //!< Cross section at each temperature - vector products_; //!< Reaction products + int mt_; //!< ENDF MT value + double q_value_; //!< Reaction Q value in [eV] + bool scatter_in_cm_; //!< scattering system in center-of-mass? + bool redundant_; //!< redundant reaction? + vector xs_; //!< Cross section at each temperature + vector products_; //!< Reaction products }; From bcc356f3a90a9599a2ee2a0ce3830c2177b03381 Mon Sep 17 00:00:00 2001 From: Gregoire Biot Date: Tue, 31 Mar 2026 10:20:17 -0400 Subject: [PATCH 20/22] Fix formatting issues in reaction.h --- include/openmc/reaction.h | 1 - 1 file changed, 1 deletion(-) diff --git a/include/openmc/reaction.h b/include/openmc/reaction.h index 7389c69c1a6..e95db1665d7 100644 --- a/include/openmc/reaction.h +++ b/include/openmc/reaction.h @@ -64,7 +64,6 @@ class Reaction { bool redundant_; //!< redundant reaction? vector xs_; //!< Cross section at each temperature vector products_; //!< Reaction products - }; //============================================================================== From e45ef34a8cdc52c07800a4ecac005ce9004e037e Mon Sep 17 00:00:00 2001 From: Grego01 Date: Tue, 31 Mar 2026 12:45:46 -0400 Subject: [PATCH 21/22] cleaning file for reviewing --- openmc/data/xs_covariance_njoy.py | 217 ++++++------------------------ 1 file changed, 44 insertions(+), 173 deletions(-) diff --git a/openmc/data/xs_covariance_njoy.py b/openmc/data/xs_covariance_njoy.py index 9c8db342ae6..40b18ec4c06 100644 --- a/openmc/data/xs_covariance_njoy.py +++ b/openmc/data/xs_covariance_njoy.py @@ -21,7 +21,7 @@ import logging from dataclasses import dataclass, field from pathlib import Path -from typing import Any, Dict, List, Optional, Sequence, Tuple, Union +from typing import Any, Dict, List, Optional, Sequence, Tuple import numpy as np import scipy.linalg as la @@ -170,13 +170,14 @@ def flush(): return reactions + # ----------------------------------------------------------------------------- -# Cholesky factor computation (with eigendecomposition fallback) +# Covariance factor computation (eigendecomposition + QR) # ----------------------------------------------------------------------------- @dataclass -class CholeskyResult: - """Container for a Cholesky-like factorization and its diagnostics. +class CovFactorResult: + """Container for a covariance factorization and its diagnostics. Attributes ---------- @@ -189,29 +190,12 @@ class CholeskyResult: Original matrix dimension G. method : str Which path was taken: "cholesky", "eigen_qr", or "zero_matrix". - condition_number : float - Ratio of largest to smallest *positive* singular value. - np.inf when the matrix is singular. - negative_eig_mass : float - Sum of absolute values of negative eigenvalues (before clamping). - Zero means the raw matrix was already PSD. - negative_eig_mass_pct : float - negative_eig_mass as a percentage of the Frobenius norm. - regularization_applied : float - The diagonal-loading factor actually used (0.0 if none). - nonzero_variance_indices : np.ndarray - Integer indices of rows/columns that had non-zero variance. """ L: np.ndarray effective_rank: int full_size: int method: str = "eigen_qr" - condition_number: float = np.inf - negative_eig_mass: float = 0.0 - negative_eig_mass_pct: float = 0.0 - regularization_applied: float = 0.0 - nonzero_variance_indices: np.ndarray = field(default_factory=lambda: np.array([], dtype=int)) - + def _enforce_symmetry(A: np.ndarray) -> np.ndarray: """Force exact symmetry: A_sym = (A + A^T) / 2.""" @@ -219,14 +203,8 @@ def _enforce_symmetry(A: np.ndarray) -> np.ndarray: def _strip_zero_variance(A: np.ndarray): - """Remove rows/columns with zero diagonal (zero variance). - - Returns - ------- - A_reduced : np.ndarray - Square sub-matrix with non-zero-variance entries only. - nz_indices : np.ndarray - Integer indices into the original matrix for the kept rows/cols. + """ + Remove rows/columns with zero diagonal (zero variance). """ diag = np.diag(A) nz = np.flatnonzero(diag > 0.0) @@ -234,54 +212,13 @@ def _strip_zero_variance(A: np.ndarray): return np.empty((0, 0), dtype=A.dtype), nz return A[np.ix_(nz, nz)], nz - -def _regularize(A: np.ndarray, correction: float) -> np.ndarray: - """Add scaled diagonal loading: A_reg = A + correction * diag(A).""" - D = np.diag(A).copy() - return A + np.diag(D * correction) - - -def _compute_diagnostics(eigenvalues: np.ndarray, fro_norm: float) -> dict: - """Compute negative-eigenvalue mass and condition number.""" - eig_neg = eigenvalues[eigenvalues < 0.0] - neg_mass = float(-np.sum(eig_neg)) - - eig_pos = eigenvalues[eigenvalues > 0.0] - if len(eig_pos) >= 2: - cond = float(eig_pos.max() / eig_pos.min()) - elif len(eig_pos) == 1: - cond = 1.0 - else: - cond = np.inf - - neg_pct = 100.0 * neg_mass / fro_norm if fro_norm > 0 else 0.0 - return dict( - negative_eig_mass=neg_mass, - negative_eig_mass_pct=neg_pct, - condition_number=cond, - ) - - -def compute_cholesky_factor( +def compute_covariance_factor( cov_matrix: np.ndarray, tol: float = 1e-10, - regularization: float = 0.0, -) -> CholeskyResult: - """Compute a lower-triangular factor L such that L @ L.T ≈ Sigma. - Parameters - ---------- - cov_matrix : np.ndarray - Square covariance matrix (G * G). - tol : float - Eigenvalue tolerance: eigenvalues with λ / λ_max < tol are - zeroed. - regularization : float - Diagonal-loading fraction (0.0 = none, 0.005 = 0.5 %). - - Returns - ------- - CholeskyResult +) -> CovFactorResult: + """ + Compute a lower-triangular factor L such that L @ L.T ≈ Sigma. """ G = cov_matrix.shape[0] @@ -292,7 +229,7 @@ def compute_cholesky_factor( Ar, nz = _strip_zero_variance(A) if Ar.size == 0: - return CholeskyResult( + return CovFactorResult( L=np.zeros((G, 0), dtype=np.float64), effective_rank=0, full_size=G, @@ -301,16 +238,9 @@ def compute_cholesky_factor( nonzero_variance_indices=nz, ) - # --- 2. Optional regularization --- - reg_applied = 0.0 - if regularization > 0.0: - Ar = _regularize(Ar, regularization) - reg_applied = regularization - Gr = Ar.shape[0] - fro_norm = float(np.linalg.norm(Ar, "fro")) - # --- 3a. Fast path: standard Cholesky --- + # --- 2a. Fast path: standard Cholesky --- try: Lr = la.cholesky(Ar, lower=True) eig_pos = np.diag(Lr) ** 2 @@ -319,34 +249,18 @@ def compute_cholesky_factor( L_full = np.zeros((G, Gr), dtype=np.float64) L_full[nz, :] = Lr - return CholeskyResult( + return CovFactorResult( L=L_full, effective_rank=Gr, full_size=G, - method="cholesky", - condition_number=cond, - negative_eig_mass=0.0, - negative_eig_mass_pct=0.0, - regularization_applied=reg_applied, - nonzero_variance_indices=nz, + method="cholesky" ) except la.LinAlgError: pass - # --- 3b. Fallback: eigendecomposition + QR --- + # --- 2b. Fallback: eigendecomposition + QR --- eigenvalues, V = la.eigh(Ar) - diag_info = _compute_diagnostics(eigenvalues, fro_norm) - - if diag_info["negative_eig_mass"] > 0: - log.warning( - "Non-PSD block (size %d): negative eigenvalue mass = %.3e " - "(%.2f%% of Frobenius norm). Clamping to zero.", - Gr, - diag_info["negative_eig_mass"], - diag_info["negative_eig_mass_pct"], - ) - max_eig = eigenvalues.max() if eigenvalues.max() > 0 else 1.0 eigenvalues[eigenvalues / max_eig < tol] = 0.0 eigenvalues[eigenvalues < 0.0] = 0.0 @@ -355,16 +269,11 @@ def compute_cholesky_factor( r = int(pos_mask.sum()) if r == 0: - return CholeskyResult( + return CovFactorResult( L=np.zeros((G, 0), dtype=np.float64), effective_rank=0, full_size=G, method="eigen_qr", - condition_number=np.inf, - negative_eig_mass=diag_info["negative_eig_mass"], - negative_eig_mass_pct=diag_info["negative_eig_mass_pct"], - regularization_applied=reg_applied, - nonzero_variance_indices=nz, ) V_pos = V[:, pos_mask] @@ -377,19 +286,11 @@ def compute_cholesky_factor( L_full = np.zeros((G, r), dtype=np.float64) L_full[nz, :] = Lr - eig_pos_vals = eigenvalues[pos_mask] - cond = float(eig_pos_vals.max() / eig_pos_vals.min()) if eig_pos_vals.min() > 0 else np.inf - - return CholeskyResult( + return CovFactorResult( L=L_full, effective_rank=r, full_size=G, method="eigen_qr", - condition_number=cond, - negative_eig_mass=diag_info["negative_eig_mass"], - negative_eig_mass_pct=diag_info["negative_eig_mass_pct"], - regularization_applied=reg_applied, - nonzero_variance_indices=nz, ) # ----------------------------------------------------------------------------- @@ -405,7 +306,7 @@ class NeutronXSCovariances: reactions: Dict[int, Dict[str, Any]] # MT -> parsed dict from ERRORR mat: Optional[int] = None temperature_k: Optional[float] = None # what we processed at (single T) - cholesky_results: Optional[Dict[int, Dict[int, CholeskyResult]]] = None + factor_results: Optional[Dict[int, Dict[int, CovFactorResult]]] = None @classmethod def from_endf( @@ -417,9 +318,8 @@ def from_endf( mat: Optional[int] = None, temperature: float = 293.6, name: Optional[str] = None, - compute_cholesky: bool = True, + compute_factors: bool = True, eig_tol: float = 1e-10, - regularization: float = 0.0, ) -> "NeutronXSCovariances": res = generate_errorr_mf33( endf_path, @@ -433,19 +333,18 @@ def from_endf( reactions = parse_errorr_mf33_text(res["tape33"], mat=mat_used) - # Pre-compute Cholesky factors for all covariance sub-blocks - chol: Optional[Dict[int, Dict[int, CholeskyResult]]] = None - if compute_cholesky: - chol = {} + # Pre-compute covariance factors for all sub-blocks + fcache: Optional[Dict[int, Dict[int, CovFactorResult]]] = None + if compute_factors: + fcache = {} for mt_key, sec in reactions.items(): - chol[mt_key] = {} + fcache[mt_key] = {} for mt1, M in sec.get("COVS", {}).items(): - result = compute_cholesky_factor( + result = compute_covariance_factor( np.asarray(M, dtype=np.float64), tol=eig_tol, - regularization=regularization, ) - chol[mt_key][mt1] = result + fcache[mt_key][mt1] = result log.info( " MT %s -> MT1 %s: method=%-9s rank=%d/%d " "cond=%.2e neg_mass=%.2e%%", @@ -466,7 +365,7 @@ def from_endf( reactions=reactions, mat=mat_used, temperature_k=float(temperature), - cholesky_results=chol, + factor_results=fcache, ) # ----------------------------------------------------------------- @@ -474,11 +373,8 @@ def from_endf( # ----------------------------------------------------------------- def to_hdf5(self, filename: str | Path, store_raw_covariance: bool = True) -> None: - """Write covariances to a standalone HDF5 file. - - Layout matches the "covariance/mf33" group that - :meth:'IncidentNeutron.export_to_hdf5' produces, so that both - paths share one schema and one reader. + """ + Write covariances to a standalone HDF5 file. """ import h5py @@ -489,23 +385,9 @@ def to_hdf5(self, filename: str | Path, store_raw_covariance: bool = True) -> No self.write_mf33_group(f, store_raw_covariance=store_raw_covariance) def write_mf33_group(self, h5_group, store_raw_covariance: bool = True, - eig_tol: float = 1e-10, regularization: float = 0.0,) -> None: - """Write the "mf33" sub-tree into an already-open HDF5 group. - - Parameters - ---------- - h5_group : h5py.Group - Parent group. A child group called ``mf33`` will be created - (or replaced) directly under it. - store_raw_covariance - If True, store full covariance matrices under ``reactions/``. - If False, only Cholesky factors are written. - eig_tol - Eigenvalue tolerance — only used as fallback when - ``self.cholesky_results`` is None. - regularization - Diagonal-loading fraction — only used as fallback when - ``self.cholesky_results`` is None. + eig_tol: float = 1e-10) -> None: + """ + Write the "mf33" sub-tree into an already-open HDF5 group. """ if "mf33" in h5_group: del h5_group["mf33"] @@ -515,7 +397,7 @@ def write_mf33_group(self, h5_group, store_raw_covariance: bool = True, mf33.attrs["source"] = np.bytes_("njoy errorr") mf33.attrs["relative"] = 1 # int flag – portable across languages mf33.attrs["store_raw_covariance"] = int(store_raw_covariance) - mf33.attrs["cholesky_storage"] = np.bytes_("thin_rank_r") + mf33.attrs["factorization"] = np.bytes_("eigen_qr_thin") if self.mat is not None: mf33.attrs["mat"] = int(self.mat) if self.temperature_k is not None: @@ -532,9 +414,9 @@ def write_mf33_group(self, h5_group, store_raw_covariance: bool = True, if store_raw_covariance: greact = mf33.create_group("reactions") - gchol = mf33.create_group("cholesky") + gfact = mf33.create_group("factors") - cached = self.cholesky_results or {} + cached = self.factor_results or {} for mt, sec in self.reactions.items(): mt_str = str(int(mt)) @@ -545,10 +427,10 @@ def write_mf33_group(self, h5_group, store_raw_covariance: bool = True, if attr_name in sec: gmt.attrs[attr_name] = float(sec[attr_name]) - gmt_chol = gchol.create_group(mt_str) + gmt_fact = gfact.create_group(mt_str) for attr_name in ("ZA", "AWR"): if attr_name in sec: - gmt_chol.attrs[attr_name] = float(sec[attr_name]) + gmt_fact.attrs[attr_name] = float(sec[attr_name]) covs: Dict[int, np.ndarray] = sec.get("COVS", {}) for mt1, M in covs.items(): @@ -564,26 +446,15 @@ def write_mf33_group(self, h5_group, store_raw_covariance: bool = True, shuffle=True, ) - # ---- Cholesky factor (always stored) ---- + # ---- Triangular factor (always stored) ---- result = cached.get(mt, {}).get(mt1, None) if result is None: - result = compute_cholesky_factor( - M_arr, tol=eig_tol, regularization=regularization, - ) + result = compute_covariance_factor( + M_arr, tol=eig_tol,) - ds = gmt_chol.create_dataset( + gmt_fact.create_dataset( ds_name, data=result.L, compression="gzip", shuffle=True, - ) - - # Diagnostic metadata - ds.attrs["effective_rank"] = result.effective_rank - ds.attrs["full_size"] = result.full_size - ds.attrs["method"] = np.bytes_(result.method) - ds.attrs["condition_number"] = result.condition_number - ds.attrs["negative_eig_mass"] = result.negative_eig_mass - ds.attrs["negative_eig_mass_pct"] = result.negative_eig_mass_pct - ds.attrs["regularization"] = result.regularization_applied - ds.attrs["nz_indices"] = result.nonzero_variance_indices \ No newline at end of file + ) \ No newline at end of file From 8aa2e76270452a0ae646ac1650ce042e5c71bbc0 Mon Sep 17 00:00:00 2001 From: Grego01 Date: Wed, 1 Apr 2026 10:02:36 -0400 Subject: [PATCH 22/22] simplified python file --- openmc/data/xs_covariance_njoy.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/data/xs_covariance_njoy.py b/openmc/data/xs_covariance_njoy.py index 40b18ec4c06..79ba67d176e 100644 --- a/openmc/data/xs_covariance_njoy.py +++ b/openmc/data/xs_covariance_njoy.py @@ -397,7 +397,7 @@ def write_mf33_group(self, h5_group, store_raw_covariance: bool = True, mf33.attrs["source"] = np.bytes_("njoy errorr") mf33.attrs["relative"] = 1 # int flag – portable across languages mf33.attrs["store_raw_covariance"] = int(store_raw_covariance) - mf33.attrs["factorization"] = np.bytes_("eigen_qr_thin") + #mf33.attrs["factorization"] = np.bytes_("eigen_qr_thin") if self.mat is not None: mf33.attrs["mat"] = int(self.mat) if self.temperature_k is not None: