From 9a29211959eb3525b61fda74f04b9edce3aa4a7b Mon Sep 17 00:00:00 2001 From: shimwell Date: Fri, 26 Jun 2026 14:57:23 +0200 Subject: [PATCH 01/14] WIP: transport-free multigroup vector-XS collapse core First slice of the convert_to_multigroup(method="transport_free") feature (see shimwell/openmc#112). Deterministic per-material collapse of CE data against a narrow-resonance self-shielded flux phi = (1/E)/Sigma_t, no Monte Carlo / nparticles. Vector XS only (total/absorption/capture/fission); scatter matrices, the source-shaped + slowing-down weighting, and the convert_to_multigroup wiring come next. --- openmc/mgxs/transport_free.py | 128 ++++++++++++++++++++++++++++++++++ 1 file changed, 128 insertions(+) create mode 100644 openmc/mgxs/transport_free.py diff --git a/openmc/mgxs/transport_free.py b/openmc/mgxs/transport_free.py new file mode 100644 index 00000000000..8cb657ae1d8 --- /dev/null +++ b/openmc/mgxs/transport_free.py @@ -0,0 +1,128 @@ +"""Deterministic, transport-free multigroup cross-section collapse. + +Collapses continuous-energy nuclear data against an *assumed* weighting flux +with narrow-resonance (Bondarenko) self-shielding -- no Monte Carlo, no +transport solve, no ``nparticles``. Each material is collapsed *directly* +(its own macroscopic flux); per-nuclide multigroup data are never combined, +because multigroup cross sections are flux-weighted averages and do not add +cleanly (see shimwell/openmc#112, design principle 7). + +This is the Phase-1 vector-cross-section core (total / absorption / capture / +fission / nu-fission). Scatter matrices and the 0-D slowing-down weighting +(option 3) come later. +""" +from __future__ import annotations + +import numpy as np + +import openmc +import openmc.data +from .groups import EnergyGroups + +# numpy >= 2.0 renamed trapz -> trapezoid +_trapz = getattr(np, "trapezoid", None) or np.trapz + + +def _as_energy_groups(groups) -> EnergyGroups: + from openmc.mgxs import GROUP_STRUCTURES # defined in openmc/mgxs/__init__.py + if isinstance(groups, EnergyGroups): + return groups + if isinstance(groups, str): + return EnergyGroups(GROUP_STRUCTURES[groups]) + return EnergyGroups(np.asarray(groups, dtype=float)) + + +def _nearest_temperature(inc: "openmc.data.IncidentNeutron", temperature: float) -> str: + temps = np.array([float(t[:-1]) for t in inc.temperatures]) + return inc.temperatures[int(np.argmin(np.abs(temps - temperature)))] + + +def _macroscopic(incs, dens, temp_str, grid, mt): + """Macroscopic pointwise xs (1/cm) for reaction ``mt`` on ``grid``; None if absent.""" + total = np.zeros_like(grid) + present = False + for nuc, n in dens.items(): + inc = incs[nuc] + try: + xs = inc[mt].xs[temp_str[nuc]] + except KeyError: + continue + total += n * xs(grid) + present = True + return total if present else None + + +def collapse_material(material, groups, temperature=294.0, cross_sections=None, + self_shield=True): + """Transport-free macroscopic multigroup cross sections for one material. + + Parameters + ---------- + material : openmc.Material + groups : EnergyGroups | str | sequence of float + temperature : float + Target temperature [K] (nearest available data temperature is used). + self_shield : bool + If True, use the narrow-resonance self-shielded flux phi = w(E)/Sigma_t(E); + if False, use the unshielded smooth flux phi = w(E) (infinite dilution). + + Returns + ------- + dict with 'group_edges' (eV, ascending) and macroscopic group XS arrays + (1/cm), ordered group 1 = highest energy (OpenMC convention). + """ + groups = _as_energy_groups(groups) + edges = np.asarray(groups.group_edges, dtype=float) + emin, emax = edges[0], edges[-1] + + if cross_sections is None: + cross_sections = openmc.config['cross_sections'] + datalib = openmc.data.DataLibrary.from_xml(cross_sections) + + dens = material.get_nuclide_atom_densities() # {nuclide: atom/b-cm} + + incs, temp_str, grids = {}, {}, [edges] + for nuc in dens: + entry = datalib.get_by_material(nuc, data_type='neutron') + inc = openmc.data.IncidentNeutron.from_hdf5(entry['path']) + incs[nuc] = inc + ts = _nearest_temperature(inc, temperature) + temp_str[nuc] = ts + grids.append(np.asarray(inc.energy[ts])) + + grid = np.unique(np.concatenate(grids)) + grid = grid[(grid >= emin) & (grid <= emax)] + + sigma_t = _macroscopic(incs, dens, temp_str, grid, 1) + if sigma_t is None: + raise ValueError("no total cross section (MT=1) found for material") + + # Weighting flux: smooth part w(E) = 1/E (asymptotic slowing-down), + # narrow-resonance self-shielded by the material's own total. + w = 1.0 / np.clip(grid, 1e-11, None) + phi = w / np.clip(sigma_t, 1e-30, None) if self_shield else w + + reactions = { + 'total': sigma_t, + 'absorption': _macroscopic(incs, dens, temp_str, grid, 101), + 'capture': _macroscopic(incs, dens, temp_str, grid, 102), + 'fission': _macroscopic(incs, dens, temp_str, grid, 18), + } + + G = groups.num_groups + out = {'group_edges': edges} + for name, sig in reactions.items(): + if sig is None: + continue + sig_g = np.zeros(G) + for g in range(G): # g: ascending in energy + lo, hi = edges[g], edges[g + 1] + m = (grid >= lo) & (grid <= hi) + if m.sum() < 2: + continue + x, p, s = grid[m], phi[m], sig[m] + den = _trapz(p, x) + sig_g[g] = _trapz(s * p, x) / den if den > 0 else 0.0 + # OpenMC orders groups with group 1 = highest energy → reverse + out[name] = sig_g[::-1] + return out From b4991205a94e29f9d33a2e056a2046aa51e22739 Mon Sep 17 00:00:00 2001 From: shimwell Date: Fri, 26 Jun 2026 15:20:18 +0200 Subject: [PATCH 02/14] Add source-aware weighting (optional source term) to transport-free collapse w(E) = 1/E + normalized source PDF (Normal/Discrete/Tabular/Mixture), so the fast groups are weighted by the actual DD/DT-Muir/TT source when provided; source remains optional (generic 1/E default). --- openmc/mgxs/transport_free.py | 28 ++++++++++++++++++++++++++-- 1 file changed, 26 insertions(+), 2 deletions(-) diff --git a/openmc/mgxs/transport_free.py b/openmc/mgxs/transport_free.py index 8cb657ae1d8..e3d8bd054cb 100644 --- a/openmc/mgxs/transport_free.py +++ b/openmc/mgxs/transport_free.py @@ -37,6 +37,26 @@ def _nearest_temperature(inc: "openmc.data.IncidentNeutron", temperature: float) return inc.temperatures[int(np.argmin(np.abs(temps - temperature)))] +def _source_pdf(dist, grid): + """Evaluate a source energy distribution as a normalized PDF on ``grid``.""" + import openmc.stats as st + pdf = np.zeros_like(grid) + if isinstance(dist, st.Normal): # e.g. openmc.stats.muir() + mu, sig = float(dist.mean_value), float(dist.std_dev) + pdf = np.exp(-0.5 * ((grid - mu) / sig) ** 2) / (sig * np.sqrt(2 * np.pi)) + elif isinstance(dist, st.Discrete): # mono lines (e.g. DD 2.45 MeV) + for xi, pi in zip(np.atleast_1d(dist.x), np.atleast_1d(dist.p)): + j = int(np.clip(np.searchsorted(grid, xi), 1, len(grid) - 1)) + pdf[j - 1] += pi / max(grid[j] - grid[j - 1], 1e-30) + elif isinstance(dist, st.Tabular): # TT continuum / arbitrary + pdf = np.interp(grid, dist.x, dist.p, left=0.0, right=0.0) + elif isinstance(dist, st.Mixture): # mixtures + for p, d in zip(dist.probability, dist.distribution): + pdf = pdf + p * _source_pdf(d, grid) + integral = _trapz(pdf, grid) + return pdf / integral if integral > 0 else pdf + + def _macroscopic(incs, dens, temp_str, grid, mt): """Macroscopic pointwise xs (1/cm) for reaction ``mt`` on ``grid``; None if absent.""" total = np.zeros_like(grid) @@ -53,7 +73,7 @@ def _macroscopic(incs, dens, temp_str, grid, mt): def collapse_material(material, groups, temperature=294.0, cross_sections=None, - self_shield=True): + self_shield=True, source=None): """Transport-free macroscopic multigroup cross sections for one material. Parameters @@ -97,9 +117,13 @@ def collapse_material(material, groups, temperature=294.0, cross_sections=None, if sigma_t is None: raise ValueError("no total cross section (MT=1) found for material") - # Weighting flux: smooth part w(E) = 1/E (asymptotic slowing-down), + # Weighting flux: smooth part w(E) = 1/E (asymptotic slowing-down), optionally + # sharpened in the fast groups by the source spectrum (added as a normalized + # PDF — at high E the source dominates 1/E, below it 1/E dominates), then # narrow-resonance self-shielded by the material's own total. w = 1.0 / np.clip(grid, 1e-11, None) + if source is not None: + w = w + _source_pdf(source, grid) phi = w / np.clip(sigma_t, 1e-30, None) if self_shield else w reactions = { From 8f079036003a3794b37d396859f2dab9224e0166 Mon Sep 17 00:00:00 2001 From: shimwell Date: Fri, 26 Jun 2026 17:21:42 +0200 Subject: [PATCH 03/14] Add 0-D slowing-down weighting (option 3) to transport-free collapse Solve the infinite-medium slowing-down balance with real energy-transfer kernels from openmc.data (analytic elastic, discrete inelastic levels, tabulated continuum / (n,xn)) on a coarse lethargy grid via an exact high->low sweep, then self-shield on the fine total. Selectable via weighting='slowing_down' (NR remains the default). Empirically this does not beat NR on total XS, but the transfer kernels are the basis for the scatter matrix. --- openmc/mgxs/transport_free.py | 166 ++++++++++++++++++++++++++++++++-- 1 file changed, 157 insertions(+), 9 deletions(-) diff --git a/openmc/mgxs/transport_free.py b/openmc/mgxs/transport_free.py index e3d8bd054cb..edb0a5c6e03 100644 --- a/openmc/mgxs/transport_free.py +++ b/openmc/mgxs/transport_free.py @@ -72,8 +72,148 @@ def _macroscopic(incs, dens, temp_str, grid, mt): return total if present else None +def _outgoing(dist): + """Classify a secondary-neutron energy distribution for the transfer kernel. + + Returns ('delta', threshold, mass_ratio) for a discrete inelastic level + (E_out = mass_ratio*(E_in - threshold)), or ('tab', incident_energies, + [Tabular,...]) for a tabulated continuum / (n,xn) distribution, or None. + """ + try: + from openmc.data import (UncorrelatedAngleEnergy, CorrelatedAngleEnergy, + LevelInelastic) + except Exception: + return None + if isinstance(dist, CorrelatedAngleEnergy): + return ('tab', np.asarray(dist.energy, dtype=float), dist.energy_out) + ed = dist.energy if isinstance(dist, UncorrelatedAngleEnergy) else getattr(dist, 'energy', None) + if isinstance(ed, LevelInelastic): + return ('delta', float(ed.threshold), float(ed.mass_ratio)) + if ed is not None and hasattr(ed, 'energy_out') and hasattr(ed, 'energy'): + return ('tab', np.asarray(ed.energy, dtype=float), ed.energy_out) + return None + + +def _slowing_down_weight(incs, dens, temp_str, fine_grid, sigma_t_fine, + source, per_decade=40): + """0-D infinite-medium slowing-down weighting flux on ``fine_grid``. + + Solves the energy-domain neutron balance + Sigma_t(E) phi(E) = S(E) + sum_r integral Sigma_s,r(E') f_r(E'->E) phi(E') dE' + on a coarse lethargy grid (the slowing-down *source* is smooth), using real + energy-transfer kernels from ``openmc.data``: analytic elastic (MT=2), + discrete inelastic levels (MT=51-90), and tabulated continuum / (n,xn) + (MT=91/16/17). Strictly down-scatter -> a single high->low energy sweep is + exact. The smooth slowing-down source is then divided by the *fine* total to + reintroduce resonance self-shielding. No transport, no Monte Carlo. + + Thermal up-scatter (S(alpha,beta)) is not modelled -> valid in the + fast/epithermal range (the populated range for fast fusion shields). + """ + emin, emax = fine_grid[0], fine_grid[-1] + nb = int(max(per_decade * np.log10(emax / emin), 40)) + be = np.logspace(np.log10(emin), np.log10(emax), nb + 1) # bin edges + cg = np.sqrt(be[:-1] * be[1:]) # bin centres + nb = len(cg) + + # dilute (1/E-weighted) bin-averaged total, smooth -> no erratic resonance sampling + w = 1.0 / np.clip(fine_grid, 1e-11, None) + def _cum(y): + c = np.zeros_like(fine_grid) + c[1:] = np.cumsum(0.5 * (y[1:] + y[:-1]) * np.diff(fine_grid)) + return c + Ni = np.interp(be, fine_grid, _cum(sigma_t_fine * w)) + Di = np.interp(be, fine_grid, _cum(w)) + sigt_c = np.diff(Ni) / np.clip(np.diff(Di), 1e-300, None) + + # transfer matrix M[l, k] = scatter rate into bin l per unit flux in bin k + M = np.zeros((nb, nb)) + for nuc, dn in dens.items(): + inc = incs[nuc] + ts = temp_str[nuc] + alpha = ((inc.atomic_weight_ratio - 1.0) / (inc.atomic_weight_ratio + 1.0)) ** 2 + for mt, r in inc.reactions.items(): + is_el = (mt == 2) + if not (is_el or (51 <= mt <= 91) or mt in (16, 17)): + continue + try: + sig = dn * r.xs[ts](cg) + except Exception: + continue + if not np.any(sig > 0): + continue + if is_el: + lo = alpha * cg + for k in range(nb): + if sig[k] <= 0: + continue + width = cg[k] - lo[k] + if width <= 0: # alpha ~ 1 (heavy): no loss + M[k, k] += sig[k] + continue + ov = np.clip(np.minimum(be[1:], cg[k]) - np.maximum(be[:-1], lo[k]), 0, None) + M[:, k] += sig[k] * ov / width + continue + prods = [p for p in r.products if p.particle == 'neutron'] + if not prods: + continue + try: + og = _outgoing(prods[0].distribution[0]) + except Exception: + og = None + if og is None: + continue + try: + mult = np.atleast_1d(np.asarray(prods[0].yield_(cg), dtype=float)) + if mult.size == 1: + mult = np.full(nb, float(mult[0])) + except Exception: + mult = np.ones(nb) + if og[0] == 'delta': + thr, mr = og[1], og[2] + eout = mr * (cg - thr) + for k in range(nb): + if sig[k] <= 0 or cg[k] <= thr or eout[k] < be[0]: + continue + l = min(max(np.searchsorted(be, eout[k]) - 1, 0), nb - 1) + M[l, k] += sig[k] * mult[k] + else: # 'tab' + ein, eos = og[1], og[2] + for k in range(nb): + if sig[k] <= 0 or cg[k] < ein[0]: + continue + t = eos[min(np.searchsorted(ein, cg[k]), len(eos) - 1)] + cx, cp = np.asarray(t.x), np.asarray(t.p) + cc = np.zeros_like(cx) + cc[1:] = np.cumsum(0.5 * (cp[1:] + cp[:-1]) * np.diff(cx)) + if cc[-1] <= 0: + continue + Wb = np.diff(np.interp(be, cx, cc / cc[-1], left=0.0, right=1.0)) + M[:, k] += sig[k] * mult[k] * Wb + + # external (or generic top-energy) source on the coarse grid + S = _source_pdf(source, cg) if source is not None else np.zeros(nb) + if S.sum() <= 0: + S = np.zeros(nb) + S[-1] = 1.0 + + # exact downward sweep (strictly down-scatter) + phi = np.zeros(nb) + for k in range(nb - 1, -1, -1): + inscat = M[k, k + 1:].dot(phi[k + 1:]) if k + 1 < nb else 0.0 + denom = sigt_c[k] - M[k, k] # remove within-bin self-scatter + phi[k] = (S[k] + inscat) / (denom if denom > 1e-30 else max(sigt_c[k], 1e-30)) + + qtot = phi * sigt_c # smooth slowing-down source + good = qtot > 0 + if good.sum() < 2: # degenerate -> fall back to 1/E + return (1.0 / np.clip(fine_grid, 1e-11, None)) / np.clip(sigma_t_fine, 1e-30, None) + qf = np.exp(np.interp(np.log(fine_grid), np.log(cg[good]), np.log(qtot[good]))) + return qf / np.clip(sigma_t_fine, 1e-30, None) # self-shield on the fine total + + def collapse_material(material, groups, temperature=294.0, cross_sections=None, - self_shield=True, source=None): + self_shield=True, source=None, weighting='nr', sd_per_decade=40): """Transport-free macroscopic multigroup cross sections for one material. Parameters @@ -117,14 +257,22 @@ def collapse_material(material, groups, temperature=294.0, cross_sections=None, if sigma_t is None: raise ValueError("no total cross section (MT=1) found for material") - # Weighting flux: smooth part w(E) = 1/E (asymptotic slowing-down), optionally - # sharpened in the fast groups by the source spectrum (added as a normalized - # PDF — at high E the source dominates 1/E, below it 1/E dominates), then - # narrow-resonance self-shielded by the material's own total. - w = 1.0 / np.clip(grid, 1e-11, None) - if source is not None: - w = w + _source_pdf(source, grid) - phi = w / np.clip(sigma_t, 1e-30, None) if self_shield else w + # Weighting flux. Two options: + # 'nr' narrow-resonance: phi = w(E)/Sigma_t(E) with smooth part + # w = 1/E (+ source PDF in the fast groups). Cheap; assumes the + # between-collision flux recovers 1/E. + # 'slowing_down' option 3: solve the 0-D infinite-medium slowing-down balance + # with real elastic+inelastic+(n,xn) transfer kernels, then + # self-shield on the fine total. Most accurate (the deterministic + # equivalent of the MC infinite_medium spectrum). + if weighting == 'slowing_down': + phi = _slowing_down_weight(incs, dens, temp_str, grid, sigma_t, source, + per_decade=sd_per_decade) + else: + w = 1.0 / np.clip(grid, 1e-11, None) + if source is not None: + w = w + _source_pdf(source, grid) + phi = w / np.clip(sigma_t, 1e-30, None) if self_shield else w reactions = { 'total': sigma_t, From 2d293e5fed7e1a67e38f786cfb01c28bc2dbabf3 Mon Sep 17 00:00:00 2001 From: shimwell Date: Fri, 26 Jun 2026 17:42:43 +0200 Subject: [PATCH 04/14] Add URR probability-table self-shielding to transport-free collapse In the unresolved resonance range the pointwise data is the infinitely-dilute average, so phi=1/Sigma_t applies no self-shielding there. _apply_urr() restores the band structure from the URR probability tables (LSSF=1 factors) and computes the Bondarenko self-shielded effective total/absorption/capture, with sigma_0 the per-resonant-nuclide background from the rest of the material. On by default (use_urr=True). Empirically (VITAMIN-J-42, vs converged material_wise): tungsten total 2.54%->0.36%, absorption 2.63%->1.36%; non-URR materials unchanged. --- openmc/mgxs/transport_free.py | 83 +++++++++++++++++++++++++++++++++-- 1 file changed, 79 insertions(+), 4 deletions(-) diff --git a/openmc/mgxs/transport_free.py b/openmc/mgxs/transport_free.py index edb0a5c6e03..8922c8511f5 100644 --- a/openmc/mgxs/transport_free.py +++ b/openmc/mgxs/transport_free.py @@ -212,8 +212,69 @@ def _cum(y): return qf / np.clip(sigma_t_fine, 1e-30, None) # self-shield on the fine total +def _apply_urr(incs, dens, temp_str, grid, sigma_t_smooth, temperature): + """Unresolved-resonance self-shielding via probability tables (Bondarenko). + + In the unresolved range the pointwise data is the infinitely-dilute (smooth) + average, so phi = 1/Sigma_t applies *no* self-shielding there. The probability + tables restore the band structure: for each band b (probability p_b, micro + total sigma_t,b), the flux ~ 1/(sigma_t,b + sigma_0), giving the self-shielded + effective micro xs = sum_b p_b sigma_x,b/(sigma_t,b+sigma_0) + / sum_b p_b/(sigma_t,b+sigma_0), + with sigma_0 the per-resonant-nuclide background from the rest of the material. + Tables here are LSSF=1 factors on the smooth xs (``multiply_smooth``). + + Returns {mt: macroscopic delta array (1/cm)} to ADD to the dilute macroscopic + cross sections, for mt in (1, 101, 102). + """ + delta = {1: np.zeros_like(grid), 101: np.zeros_like(grid), 102: np.zeros_like(grid)} + for nuc, n in dens.items(): + inc = incs[nuc] + urr = getattr(inc, 'urr', None) + if not urr: + continue + cand = [t for t in urr if urr.get(t) is not None and t in inc.temperatures] + if not cand: + continue + ts = min(cand, key=lambda t: abs(float(t[:-1]) - temperature)) + pt = urr[ts] + if not getattr(pt, 'multiply_smooth', False): + continue # only LSSF=1 factor tables + e = np.asarray(pt.energy, dtype=float) + tab = np.asarray(pt.table, dtype=float) # [nE, 6, nbands] + idx = np.where((grid >= e[0]) & (grid <= e[-1]))[0] + if idx.size == 0: + continue + st = inc[1].xs[ts](grid) # smooth micro total + sc = inc[102].xs[ts](grid) # smooth micro capture + try: + sa = inc[101].xs[ts](grid) + except Exception: + sa = sc + for gi in idx: + j = int(np.clip(np.searchsorted(e, grid[gi]), 1, len(e) - 1)) + row = tab[j - 1] if abs(e[j - 1] - grid[gi]) <= abs(e[j] - grid[gi]) else tab[j] + p = np.diff(np.concatenate(([0.0], row[0]))) # band probabilities + sig0 = (sigma_t_smooth[gi] - n * st[gi]) / n # micro background (others) + if sig0 < 0: + sig0 = 0.0 + stb = st[gi] * row[1] # band micro total + w = p / (stb + sig0) + denom = w.sum() + if denom <= 0: + continue + st_eff = float((w * stb).sum() / denom) + sc_eff = float((w * (sc[gi] * row[4])).sum() / denom) + fc = sc_eff / sc[gi] if sc[gi] > 0 else 1.0 # shield absorption like capture + delta[1][gi] += n * (st_eff - st[gi]) + delta[102][gi] += n * (sc_eff - sc[gi]) + delta[101][gi] += n * sa[gi] * (fc - 1.0) + return delta + + def collapse_material(material, groups, temperature=294.0, cross_sections=None, - self_shield=True, source=None, weighting='nr', sd_per_decade=40): + self_shield=True, source=None, weighting='nr', sd_per_decade=40, + use_urr=True): """Transport-free macroscopic multigroup cross sections for one material. Parameters @@ -257,6 +318,20 @@ def collapse_material(material, groups, temperature=294.0, cross_sections=None, if sigma_t is None: raise ValueError("no total cross section (MT=1) found for material") + sigma_a = _macroscopic(incs, dens, temp_str, grid, 101) + sigma_c = _macroscopic(incs, dens, temp_str, grid, 102) + sigma_f = _macroscopic(incs, dens, temp_str, grid, 18) + + # Unresolved-resonance self-shielding (probability tables): correct the dilute + # total / absorption / capture band-by-band in the URR before weighting. + if self_shield and use_urr: + d = _apply_urr(incs, dens, temp_str, grid, sigma_t, temperature) + sigma_t = sigma_t + d[1] + if sigma_a is not None: + sigma_a = sigma_a + d[101] + if sigma_c is not None: + sigma_c = sigma_c + d[102] + # Weighting flux. Two options: # 'nr' narrow-resonance: phi = w(E)/Sigma_t(E) with smooth part # w = 1/E (+ source PDF in the fast groups). Cheap; assumes the @@ -276,9 +351,9 @@ def collapse_material(material, groups, temperature=294.0, cross_sections=None, reactions = { 'total': sigma_t, - 'absorption': _macroscopic(incs, dens, temp_str, grid, 101), - 'capture': _macroscopic(incs, dens, temp_str, grid, 102), - 'fission': _macroscopic(incs, dens, temp_str, grid, 18), + 'absorption': sigma_a, + 'capture': sigma_c, + 'fission': sigma_f, } G = groups.num_groups From ea0a95238e95c225861a96203ade1acf2c92476b Mon Sep 17 00:00:00 2001 From: shimwell Date: Fri, 26 Jun 2026 18:47:03 +0200 Subject: [PATCH 05/14] Add intermediate-resonance (IR) weighting option to transport-free collapse weighting='ir' uses phi = w/[Sigma_t - sum_i (1-lambda_i) Sigma_s,i], with a per-nuclide lambda (default mass proxy 1-alpha, overridable via ir_lambda). lambda=1 reproduces the NR flux exactly (verified to machine zero). NON-DEFAULT: NR remains the default. Empirically (VITAMIN-J-42, vs converged material_wise) the parameter-free mass-proxy lambda does NOT improve on NR+URR and over-shields absorption (steel 1.9%->7.7%, Fe56 6.9%->14.3%): 1-alpha is the scatterer's slowing-down weight, whereas absorber self-shielding needs the Goldstein-Cohen lambda which is ~1 for narrow resolved resonances. A rigorous per-group lambda needs resonance-parameter fits and is sub-percent in fast spectra anyway. Kept as an option (overridable lambda) with this caveat; the recommended configuration is NR + URR probability tables. --- openmc/mgxs/transport_free.py | 44 +++++++++++++++++++++++++++-------- 1 file changed, 34 insertions(+), 10 deletions(-) diff --git a/openmc/mgxs/transport_free.py b/openmc/mgxs/transport_free.py index 8922c8511f5..6daa69819b2 100644 --- a/openmc/mgxs/transport_free.py +++ b/openmc/mgxs/transport_free.py @@ -274,7 +274,7 @@ def _apply_urr(incs, dens, temp_str, grid, sigma_t_smooth, temperature): def collapse_material(material, groups, temperature=294.0, cross_sections=None, self_shield=True, source=None, weighting='nr', sd_per_decade=40, - use_urr=True): + use_urr=True, ir_lambda=None): """Transport-free macroscopic multigroup cross sections for one material. Parameters @@ -332,14 +332,18 @@ def collapse_material(material, groups, temperature=294.0, cross_sections=None, if sigma_c is not None: sigma_c = sigma_c + d[102] - # Weighting flux. Two options: - # 'nr' narrow-resonance: phi = w(E)/Sigma_t(E) with smooth part - # w = 1/E (+ source PDF in the fast groups). Cheap; assumes the - # between-collision flux recovers 1/E. - # 'slowing_down' option 3: solve the 0-D infinite-medium slowing-down balance - # with real elastic+inelastic+(n,xn) transfer kernels, then - # self-shield on the fine total. Most accurate (the deterministic - # equivalent of the MC infinite_medium spectrum). + # Weighting flux. Options: + # 'nr' narrow-resonance: phi = w(E)/Sigma_t(E), w = 1/E (+ source PDF). + # 'ir' intermediate resonance: phi = w(E)/[Sigma_t - sum_i (1-lambda_i) + # Sigma_s,i], i.e. only a fraction lambda_i of each nuclide's + # scattering moderates. lambda_i=1 recovers NR exactly; lambda_i=0 + # is wide-resonance. The default per-nuclide lambda is the mass + # proxy 1-alpha (alpha=((A-1)/(A+1))^2) -- a documented kinematic + # proxy for the *scatterer's* slowing-down weight, NOT the rigorous + # per-group Goldstein-Cohen parameter; override via `ir_lambda` + # {nuclide: lambda}. Applied on the (URR-corrected) Sigma_t. + # 'slowing_down' option 3: solve the 0-D slowing-down balance with real transfer + # kernels, then self-shield on the fine total. if weighting == 'slowing_down': phi = _slowing_down_weight(incs, dens, temp_str, grid, sigma_t, source, per_decade=sd_per_decade) @@ -347,7 +351,27 @@ def collapse_material(material, groups, temperature=294.0, cross_sections=None, w = 1.0 / np.clip(grid, 1e-11, None) if source is not None: w = w + _source_pdf(source, grid) - phi = w / np.clip(sigma_t, 1e-30, None) if self_shield else w + if weighting == 'ir' and self_shield: + removed = np.zeros_like(grid) # sum_i (1-lambda_i) Sigma_s,i + for nuc, n in dens.items(): + inc = incs[nuc] + A = inc.atomic_weight_ratio + lam = (ir_lambda or {}).get(nuc, 4.0 * A / (A + 1.0) ** 2) # 1 - alpha + if lam >= 1.0: + continue + ts = temp_str[nuc] + sti = inc[1].xs[ts](grid) + try: + sai = inc[101].xs[ts](grid) + except KeyError: + try: + sai = inc[102].xs[ts](grid) + except KeyError: + sai = np.zeros_like(grid) + removed += (1.0 - lam) * n * np.clip(sti - sai, 0.0, None) + phi = w / np.clip(sigma_t - removed, 1e-30, None) # positive: = Sigma_a + sum lam_i Sigma_s,i + else: + phi = w / np.clip(sigma_t, 1e-30, None) if self_shield else w reactions = { 'total': sigma_t, From 7527d5eaa380ecf6a798590724cfc9c9933a6603 Mon Sep 17 00:00:00 2001 From: shimwell Date: Fri, 26 Jun 2026 21:09:50 +0200 Subject: [PATCH 06/14] Always self-shield in transport-free collapse (drop self_shield and use_urr knobs) Self-shielding is not optional for a real material: the infinitely-dilute flux is never the more-accurate choice for transport, so neither the master self_shield toggle nor the use_urr sub-flag should be user-facing. Both are removed; resolved-resonance self-shielding (phi=w/Sigma_t) and unresolved probability-table self-shielding are now always applied. Keeps the public API knob-free (most-accurate-by-default), consistent with not exposing weight/scatter_order/correction. The raw dilute macroscopic xs remain available internally for diagnostics. --- openmc/mgxs/transport_free.py | 31 ++++++++++++++++--------------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/openmc/mgxs/transport_free.py b/openmc/mgxs/transport_free.py index 6daa69819b2..49f6de77f0a 100644 --- a/openmc/mgxs/transport_free.py +++ b/openmc/mgxs/transport_free.py @@ -273,8 +273,8 @@ def _apply_urr(incs, dens, temp_str, grid, sigma_t_smooth, temperature): def collapse_material(material, groups, temperature=294.0, cross_sections=None, - self_shield=True, source=None, weighting='nr', sd_per_decade=40, - use_urr=True, ir_lambda=None): + source=None, weighting='nr', sd_per_decade=40, + ir_lambda=None): """Transport-free macroscopic multigroup cross sections for one material. Parameters @@ -283,9 +283,10 @@ def collapse_material(material, groups, temperature=294.0, cross_sections=None, groups : EnergyGroups | str | sequence of float temperature : float Target temperature [K] (nearest available data temperature is used). - self_shield : bool - If True, use the narrow-resonance self-shielded flux phi = w(E)/Sigma_t(E); - if False, use the unshielded smooth flux phi = w(E) (infinite dilution). + + Self-shielding (resolved resonances via phi = w/Sigma_t, plus the unresolved + range via probability tables) is always applied: a real material is never at + infinite dilution, so there is no knob to turn it off. Returns ------- @@ -323,14 +324,14 @@ def collapse_material(material, groups, temperature=294.0, cross_sections=None, sigma_f = _macroscopic(incs, dens, temp_str, grid, 18) # Unresolved-resonance self-shielding (probability tables): correct the dilute - # total / absorption / capture band-by-band in the URR before weighting. - if self_shield and use_urr: - d = _apply_urr(incs, dens, temp_str, grid, sigma_t, temperature) - sigma_t = sigma_t + d[1] - if sigma_a is not None: - sigma_a = sigma_a + d[101] - if sigma_c is not None: - sigma_c = sigma_c + d[102] + # total / absorption / capture band-by-band in the URR before weighting. Always + # applied -- it is part of self-shielding, which a real material always has. + d = _apply_urr(incs, dens, temp_str, grid, sigma_t, temperature) + sigma_t = sigma_t + d[1] + if sigma_a is not None: + sigma_a = sigma_a + d[101] + if sigma_c is not None: + sigma_c = sigma_c + d[102] # Weighting flux. Options: # 'nr' narrow-resonance: phi = w(E)/Sigma_t(E), w = 1/E (+ source PDF). @@ -351,7 +352,7 @@ def collapse_material(material, groups, temperature=294.0, cross_sections=None, w = 1.0 / np.clip(grid, 1e-11, None) if source is not None: w = w + _source_pdf(source, grid) - if weighting == 'ir' and self_shield: + if weighting == 'ir': removed = np.zeros_like(grid) # sum_i (1-lambda_i) Sigma_s,i for nuc, n in dens.items(): inc = incs[nuc] @@ -371,7 +372,7 @@ def collapse_material(material, groups, temperature=294.0, cross_sections=None, removed += (1.0 - lam) * n * np.clip(sti - sai, 0.0, None) phi = w / np.clip(sigma_t - removed, 1e-30, None) # positive: = Sigma_a + sum lam_i Sigma_s,i else: - phi = w / np.clip(sigma_t, 1e-30, None) if self_shield else w + phi = w / np.clip(sigma_t, 1e-30, None) reactions = { 'total': sigma_t, From 8108d526876e17b26bb1cd77bbdbc07e340bef0e Mon Sep 17 00:00:00 2001 From: shimwell Date: Fri, 26 Jun 2026 21:38:43 +0200 Subject: [PATCH 07/14] Remove IR and slowing-down weighting options -- NR + self-shielding is the method Neither improved the MGXS, and both are out of scope for a focused MGXS generator. IR with the mass-proxy lambda over-shields absorption, and a rigorous Goldstein-Cohen lambda can neither be supplied by users nor derived automatically. Slowing-down does not beat NR on total XS; its only value was the scatter-matrix transfer kernels, which we are not pursuing (random ray consumes the vector MGXS directly). Leaves a single deterministic, knob-free, noise-free generator: collapse_material(material, groups, source=None) -- 1/E (+optional source) narrow-resonance weighting with always-on self-shielding (resolved + URR). NR+URR output is byte-identical to before; module ~140 lines lighter. --- openmc/mgxs/transport_free.py | 197 +++------------------------------- 1 file changed, 12 insertions(+), 185 deletions(-) diff --git a/openmc/mgxs/transport_free.py b/openmc/mgxs/transport_free.py index 49f6de77f0a..b330562b7ab 100644 --- a/openmc/mgxs/transport_free.py +++ b/openmc/mgxs/transport_free.py @@ -7,9 +7,9 @@ because multigroup cross sections are flux-weighted averages and do not add cleanly (see shimwell/openmc#112, design principle 7). -This is the Phase-1 vector-cross-section core (total / absorption / capture / -fission / nu-fission). Scatter matrices and the 0-D slowing-down weighting -(option 3) come later. +Vector cross sections only (total / absorption / capture / fission). Weighting +is a 1/E (+ optional source) narrow-resonance flux; self-shielding (resolved +resonances plus the unresolved range via probability tables) is always applied. """ from __future__ import annotations @@ -72,146 +72,6 @@ def _macroscopic(incs, dens, temp_str, grid, mt): return total if present else None -def _outgoing(dist): - """Classify a secondary-neutron energy distribution for the transfer kernel. - - Returns ('delta', threshold, mass_ratio) for a discrete inelastic level - (E_out = mass_ratio*(E_in - threshold)), or ('tab', incident_energies, - [Tabular,...]) for a tabulated continuum / (n,xn) distribution, or None. - """ - try: - from openmc.data import (UncorrelatedAngleEnergy, CorrelatedAngleEnergy, - LevelInelastic) - except Exception: - return None - if isinstance(dist, CorrelatedAngleEnergy): - return ('tab', np.asarray(dist.energy, dtype=float), dist.energy_out) - ed = dist.energy if isinstance(dist, UncorrelatedAngleEnergy) else getattr(dist, 'energy', None) - if isinstance(ed, LevelInelastic): - return ('delta', float(ed.threshold), float(ed.mass_ratio)) - if ed is not None and hasattr(ed, 'energy_out') and hasattr(ed, 'energy'): - return ('tab', np.asarray(ed.energy, dtype=float), ed.energy_out) - return None - - -def _slowing_down_weight(incs, dens, temp_str, fine_grid, sigma_t_fine, - source, per_decade=40): - """0-D infinite-medium slowing-down weighting flux on ``fine_grid``. - - Solves the energy-domain neutron balance - Sigma_t(E) phi(E) = S(E) + sum_r integral Sigma_s,r(E') f_r(E'->E) phi(E') dE' - on a coarse lethargy grid (the slowing-down *source* is smooth), using real - energy-transfer kernels from ``openmc.data``: analytic elastic (MT=2), - discrete inelastic levels (MT=51-90), and tabulated continuum / (n,xn) - (MT=91/16/17). Strictly down-scatter -> a single high->low energy sweep is - exact. The smooth slowing-down source is then divided by the *fine* total to - reintroduce resonance self-shielding. No transport, no Monte Carlo. - - Thermal up-scatter (S(alpha,beta)) is not modelled -> valid in the - fast/epithermal range (the populated range for fast fusion shields). - """ - emin, emax = fine_grid[0], fine_grid[-1] - nb = int(max(per_decade * np.log10(emax / emin), 40)) - be = np.logspace(np.log10(emin), np.log10(emax), nb + 1) # bin edges - cg = np.sqrt(be[:-1] * be[1:]) # bin centres - nb = len(cg) - - # dilute (1/E-weighted) bin-averaged total, smooth -> no erratic resonance sampling - w = 1.0 / np.clip(fine_grid, 1e-11, None) - def _cum(y): - c = np.zeros_like(fine_grid) - c[1:] = np.cumsum(0.5 * (y[1:] + y[:-1]) * np.diff(fine_grid)) - return c - Ni = np.interp(be, fine_grid, _cum(sigma_t_fine * w)) - Di = np.interp(be, fine_grid, _cum(w)) - sigt_c = np.diff(Ni) / np.clip(np.diff(Di), 1e-300, None) - - # transfer matrix M[l, k] = scatter rate into bin l per unit flux in bin k - M = np.zeros((nb, nb)) - for nuc, dn in dens.items(): - inc = incs[nuc] - ts = temp_str[nuc] - alpha = ((inc.atomic_weight_ratio - 1.0) / (inc.atomic_weight_ratio + 1.0)) ** 2 - for mt, r in inc.reactions.items(): - is_el = (mt == 2) - if not (is_el or (51 <= mt <= 91) or mt in (16, 17)): - continue - try: - sig = dn * r.xs[ts](cg) - except Exception: - continue - if not np.any(sig > 0): - continue - if is_el: - lo = alpha * cg - for k in range(nb): - if sig[k] <= 0: - continue - width = cg[k] - lo[k] - if width <= 0: # alpha ~ 1 (heavy): no loss - M[k, k] += sig[k] - continue - ov = np.clip(np.minimum(be[1:], cg[k]) - np.maximum(be[:-1], lo[k]), 0, None) - M[:, k] += sig[k] * ov / width - continue - prods = [p for p in r.products if p.particle == 'neutron'] - if not prods: - continue - try: - og = _outgoing(prods[0].distribution[0]) - except Exception: - og = None - if og is None: - continue - try: - mult = np.atleast_1d(np.asarray(prods[0].yield_(cg), dtype=float)) - if mult.size == 1: - mult = np.full(nb, float(mult[0])) - except Exception: - mult = np.ones(nb) - if og[0] == 'delta': - thr, mr = og[1], og[2] - eout = mr * (cg - thr) - for k in range(nb): - if sig[k] <= 0 or cg[k] <= thr or eout[k] < be[0]: - continue - l = min(max(np.searchsorted(be, eout[k]) - 1, 0), nb - 1) - M[l, k] += sig[k] * mult[k] - else: # 'tab' - ein, eos = og[1], og[2] - for k in range(nb): - if sig[k] <= 0 or cg[k] < ein[0]: - continue - t = eos[min(np.searchsorted(ein, cg[k]), len(eos) - 1)] - cx, cp = np.asarray(t.x), np.asarray(t.p) - cc = np.zeros_like(cx) - cc[1:] = np.cumsum(0.5 * (cp[1:] + cp[:-1]) * np.diff(cx)) - if cc[-1] <= 0: - continue - Wb = np.diff(np.interp(be, cx, cc / cc[-1], left=0.0, right=1.0)) - M[:, k] += sig[k] * mult[k] * Wb - - # external (or generic top-energy) source on the coarse grid - S = _source_pdf(source, cg) if source is not None else np.zeros(nb) - if S.sum() <= 0: - S = np.zeros(nb) - S[-1] = 1.0 - - # exact downward sweep (strictly down-scatter) - phi = np.zeros(nb) - for k in range(nb - 1, -1, -1): - inscat = M[k, k + 1:].dot(phi[k + 1:]) if k + 1 < nb else 0.0 - denom = sigt_c[k] - M[k, k] # remove within-bin self-scatter - phi[k] = (S[k] + inscat) / (denom if denom > 1e-30 else max(sigt_c[k], 1e-30)) - - qtot = phi * sigt_c # smooth slowing-down source - good = qtot > 0 - if good.sum() < 2: # degenerate -> fall back to 1/E - return (1.0 / np.clip(fine_grid, 1e-11, None)) / np.clip(sigma_t_fine, 1e-30, None) - qf = np.exp(np.interp(np.log(fine_grid), np.log(cg[good]), np.log(qtot[good]))) - return qf / np.clip(sigma_t_fine, 1e-30, None) # self-shield on the fine total - - def _apply_urr(incs, dens, temp_str, grid, sigma_t_smooth, temperature): """Unresolved-resonance self-shielding via probability tables (Bondarenko). @@ -273,8 +133,7 @@ def _apply_urr(incs, dens, temp_str, grid, sigma_t_smooth, temperature): def collapse_material(material, groups, temperature=294.0, cross_sections=None, - source=None, weighting='nr', sd_per_decade=40, - ir_lambda=None): + source=None): """Transport-free macroscopic multigroup cross sections for one material. Parameters @@ -333,46 +192,14 @@ def collapse_material(material, groups, temperature=294.0, cross_sections=None, if sigma_c is not None: sigma_c = sigma_c + d[102] - # Weighting flux. Options: - # 'nr' narrow-resonance: phi = w(E)/Sigma_t(E), w = 1/E (+ source PDF). - # 'ir' intermediate resonance: phi = w(E)/[Sigma_t - sum_i (1-lambda_i) - # Sigma_s,i], i.e. only a fraction lambda_i of each nuclide's - # scattering moderates. lambda_i=1 recovers NR exactly; lambda_i=0 - # is wide-resonance. The default per-nuclide lambda is the mass - # proxy 1-alpha (alpha=((A-1)/(A+1))^2) -- a documented kinematic - # proxy for the *scatterer's* slowing-down weight, NOT the rigorous - # per-group Goldstein-Cohen parameter; override via `ir_lambda` - # {nuclide: lambda}. Applied on the (URR-corrected) Sigma_t. - # 'slowing_down' option 3: solve the 0-D slowing-down balance with real transfer - # kernels, then self-shield on the fine total. - if weighting == 'slowing_down': - phi = _slowing_down_weight(incs, dens, temp_str, grid, sigma_t, source, - per_decade=sd_per_decade) - else: - w = 1.0 / np.clip(grid, 1e-11, None) - if source is not None: - w = w + _source_pdf(source, grid) - if weighting == 'ir': - removed = np.zeros_like(grid) # sum_i (1-lambda_i) Sigma_s,i - for nuc, n in dens.items(): - inc = incs[nuc] - A = inc.atomic_weight_ratio - lam = (ir_lambda or {}).get(nuc, 4.0 * A / (A + 1.0) ** 2) # 1 - alpha - if lam >= 1.0: - continue - ts = temp_str[nuc] - sti = inc[1].xs[ts](grid) - try: - sai = inc[101].xs[ts](grid) - except KeyError: - try: - sai = inc[102].xs[ts](grid) - except KeyError: - sai = np.zeros_like(grid) - removed += (1.0 - lam) * n * np.clip(sti - sai, 0.0, None) - phi = w / np.clip(sigma_t - removed, 1e-30, None) # positive: = Sigma_a + sum lam_i Sigma_s,i - else: - phi = w / np.clip(sigma_t, 1e-30, None) + # Narrow-resonance weighting flux: phi = w(E)/Sigma_t(E). The smooth part + # w = 1/E (asymptotic slowing-down) is optionally sharpened in the fast groups + # by the source spectrum, then self-shielded by the material's own + # (URR-corrected) total. + w = 1.0 / np.clip(grid, 1e-11, None) + if source is not None: + w = w + _source_pdf(source, grid) + phi = w / np.clip(sigma_t, 1e-30, None) reactions = { 'total': sigma_t, From f96b34c7e0819ebcd2fa8f12a5d1dbc6bf2c6992 Mon Sep 17 00:00:00 2001 From: shimwell Date: Fri, 26 Jun 2026 22:56:57 +0200 Subject: [PATCH 08/14] Add deterministic P0 group-to-group scattering matrix scatter_matrix(material, groups, source=None) builds the P0 Sigma_s,g->g' [G_in,G_out] (OpenMC ordering) needed to actually feed a P0 solver such as random ray -- completing the deterministic library alongside the vector XS. Kernels (no Monte Carlo): elastic with the real CM angular distribution (E_out=E*(A^2+2A mu+1)/(A+1)^2 folded over f(mu)); discrete inelastic levels spread over their lab energy range [base +/- amp]; continuum / (n,xn) via unit-base interpolation of the secondary-energy distributions, summed over all neutron products with multiplicity (yield_). A per-group panel sub-grid resolves the elastic in-group/down-scatter split where pointwise data is sparse, and upscatter artifacts are folded into the diagonal. Validated vs cached MC scatter matrices (P0): tungsten ties stochastic_slab at fine fusion groups (CCFE-709: 10.8% vs 10.6% element-wise); steel/Fe56 ~1.4-1.5x slab; row-sums competitive. Same NR+URR weighting flux as collapse_material. --- openmc/mgxs/transport_free.py | 243 +++++++++++++++++++++++++++++++++- 1 file changed, 240 insertions(+), 3 deletions(-) diff --git a/openmc/mgxs/transport_free.py b/openmc/mgxs/transport_free.py index b330562b7ab..eb579919f13 100644 --- a/openmc/mgxs/transport_free.py +++ b/openmc/mgxs/transport_free.py @@ -7,9 +7,13 @@ because multigroup cross sections are flux-weighted averages and do not add cleanly (see shimwell/openmc#112, design principle 7). -Vector cross sections only (total / absorption / capture / fission). Weighting -is a 1/E (+ optional source) narrow-resonance flux; self-shielding (resolved -resonances plus the unresolved range via probability tables) is always applied. +``collapse_material`` gives the vector cross sections (total / absorption / +capture / fission); ``scatter_matrix`` gives the P0 group-to-group scattering +matrix (elastic with the real CM angular distribution, discrete inelastic levels, +and unit-base-interpolated continuum / (n,xn) with multiplicity) -- together a +complete MG library for a P0 solver such as random ray. Weighting is a 1/E +(+ optional source) narrow-resonance flux; self-shielding (resolved resonances +plus the unresolved range via probability tables) is always applied. """ from __future__ import annotations @@ -225,3 +229,236 @@ def collapse_material(material, groups, temperature=294.0, cross_sections=None, # OpenMC orders groups with group 1 = highest energy → reverse out[name] = sig_g[::-1] return out + + +# ---------------------------------------------------------------------------- +# P0 group-to-group scattering matrix (deterministic) +# ---------------------------------------------------------------------------- + +def _curve_frac(x, p, edges): + """Group-integrate a lin-lin density (x, p) over ascending edges -> len-G, normalized.""" + x = np.asarray(x, float); p = np.asarray(p, float) + if x.size < 2: + return None + cc = np.zeros_like(x); cc[1:] = np.cumsum(0.5 * (p[1:] + p[:-1]) * np.diff(x)) + tot = cc[-1] + if tot <= 0: + return None + return np.diff(np.interp(edges, x, cc, left=0.0, right=tot)) / tot + + +def _tab_frac(t, edges): + """Group fractions for an ``openmc.stats.Tabular`` outgoing distribution.""" + try: + cc = np.asarray(t.cdf(), float); x = np.asarray(t.x, float) + if cc[-1] <= 0: + return None + return np.diff(np.interp(edges, x, cc, left=0.0, right=cc[-1])) / cc[-1] + except Exception: + return _curve_frac(getattr(t, 'x', []), getattr(t, 'p', []), edges) + + +def _unitbase_frac(d_lo, d_hi, f, edges): + """Unit-base interpolation of two Tabulars at fraction ``f``, group-integrated. + + Both outgoing distributions are mapped to the unit interval, blended there + (valid because the two span different outgoing-energy ranges), then mapped back + so every outgoing energy stays within the interpolated kinematic range. + """ + xl, pl = np.asarray(d_lo.x, float), np.asarray(d_lo.p, float) + xh, ph = np.asarray(d_hi.x, float), np.asarray(d_hi.p, float) + if xl.size < 2 or xh.size < 2: + return _tab_frac(d_lo if f < 0.5 else d_hi, edges) + Ll, Lh = xl[-1] - xl[0], xh[-1] - xh[0] + if Ll <= 0 or Lh <= 0: + return _tab_frac(d_lo if f < 0.5 else d_hi, edges) + E0 = xl[0] + f * (xh[0] - xl[0]); E1 = xl[-1] + f * (xh[-1] - xl[-1]) + if E1 <= E0: + return None + ul = (xl - xl[0]) / Ll; uh = (xh - xh[0]) / Lh + u = np.union1d(ul, uh) + pu = (1 - f) * np.interp(u, ul, pl * Ll, left=0, right=0) + f * np.interp(u, uh, ph * Lh, left=0, right=0) + return _curve_frac(E0 + u * (E1 - E0), pu, edges) + + +def _incident_code(dist, k): + """Interpolation code on the incident interval starting at index ``k`` (2=lin-lin, 1=histogram).""" + bps = getattr(dist, 'breakpoints', None); itp = getattr(dist, 'interpolation', None) + if bps is None or itp is None: + return 2 + for r, bp in enumerate(np.atleast_1d(bps)): + if k + 1 <= bp: + return int(np.atleast_1d(itp)[r]) + return 2 + + +def _tabdist(dist): + """Return the object exposing .energy (incident) + .energy_out (Tabulars), or None.""" + if hasattr(dist, 'energy_out') and hasattr(dist, 'energy') and \ + not isinstance(getattr(dist, 'energy'), openmc.data.LevelInelastic): + return dist + ed = getattr(dist, 'energy', None) + if ed is not None and hasattr(ed, 'energy_out') and hasattr(ed, 'energy'): + return ed + return None + + +def _unitbase_anchors(td, edges, per=4): + """Precompute unit-base group fractions at sub-sampled incident energies.""" + ein = np.asarray(td.energy, float); eos = td.energy_out; K = len(eos) + Ea, GF = [], [] + G = len(edges) - 1 + for k in range(K - 1): + hist = _incident_code(td, k) == 1 + for s in range(per): + f = s / per + gf = _tab_frac(eos[k], edges) if hist else _unitbase_frac(eos[k], eos[k + 1], f, edges) + Ea.append(ein[k] + f * (ein[k + 1] - ein[k])); GF.append(gf if gf is not None else np.zeros(G)) + gf = _tab_frac(eos[-1], edges) + Ea.append(ein[-1]); GF.append(gf if gf is not None else np.zeros(G)) + return np.array(Ea), np.array(GF) + + +def scatter_matrix(material, groups, temperature=294.0, cross_sections=None, source=None): + """Deterministic P0 group-to-group scattering matrix for one material. + + Returns the macroscopic Sigma_s,g->g' (1/cm) as an ``[G_in, G_out]`` array in + OpenMC ordering (group 1 = highest energy), matching ``openmc.XSdata`` P0 + (Legendre order 0). Kernels: elastic with the real CM angular distribution, + discrete inelastic levels spread over their lab energy range, and unit-base + interpolated continuum / (n,xn) distributions summed over all neutron products + with their multiplicity. No Monte Carlo. Same weighting flux as + :func:`collapse_material`. + """ + groups = _as_energy_groups(groups) + edges = np.asarray(groups.group_edges, float) + emin, emax = edges[0], edges[-1]; G = groups.num_groups + if cross_sections is None: + cross_sections = openmc.config['cross_sections'] + datalib = openmc.data.DataLibrary.from_xml(cross_sections) + dens = material.get_nuclide_atom_densities() + incs, temp_str, grids = {}, {}, [edges] + for nuc in dens: + entry = datalib.get_by_material(nuc, data_type='neutron') + inc = openmc.data.IncidentNeutron.from_hdf5(entry['path']); incs[nuc] = inc + ts = _nearest_temperature(inc, temperature); temp_str[nuc] = ts + grids.append(np.asarray(inc.energy[ts])) + grid = np.unique(np.concatenate(grids)); grid = grid[(grid >= emin) & (grid <= emax)] + # panel sub-grid: interior points in every group so the incoming integration + # (esp. the elastic in-group vs down-scatter split) is resolved where data is sparse. + sub = np.concatenate([np.geomspace(edges[g], edges[g + 1], 17)[1:-1] for g in range(G)]) + grid = np.unique(np.concatenate([grid, sub])) + + sigma_t = _macroscopic(incs, dens, temp_str, grid, 1) + sigma_t = sigma_t + _apply_urr(incs, dens, temp_str, grid, sigma_t, temperature)[1] + w = 1.0 / np.clip(grid, 1e-11, None) + if source is not None: + w = w + _source_pdf(source, grid) + phi = w / np.clip(sigma_t, 1e-30, None) + dwid = np.empty_like(grid) + dwid[1:-1] = 0.5 * (grid[2:] - grid[:-2]); dwid[0] = 0.5*(grid[1]-grid[0]); dwid[-1] = 0.5*(grid[-1]-grid[-2]) + wt = phi * dwid + gi = np.clip(np.searchsorted(edges, grid, side='right') - 1, 0, G - 1) + M = np.zeros((G, G)); denom = np.zeros(G); np.add.at(denom, gi, wt) + + def overlap(lo, hi): + if hi <= lo: + return None + return np.clip(np.minimum(edges[1:], hi) - np.maximum(edges[:-1], lo), 0, None) / (hi - lo) + + for nuc, n in dens.items(): + inc = incs[nuc]; ts = temp_str[nuc]; A = inc.atomic_weight_ratio; a1 = (A + 1.0) ** 2 + for mt, r in inc.reactions.items(): + is_el = (mt == 2) + if not (is_el or (51 <= mt <= 91) or mt in (16, 17)): + continue + try: + xs = r.xs[ts](grid) + except Exception: + continue + nz = np.nonzero(xs > 0)[0] + if nz.size == 0: + continue + src = n * xs * wt + if is_el: # elastic: CM angular distribution + try: + ang = r.products[0].distribution[0].angle; aE = np.asarray(ang.energy) + except Exception: + ang = None + for i in nz: + E = grid[i]; mu = fp = None + if ang is not None: + t = ang.mu[min(np.searchsorted(aE, E), len(aE) - 1)] + if hasattr(t, 'x') and hasattr(t, 'p'): + mu = np.asarray(t.x, float); fp = np.asarray(t.p, float) + if mu is None or mu.size < 2: + mu = np.linspace(-1, 1, 33); fp = np.full(33, 0.5) + wmu = np.empty_like(mu); wmu[1:-1] = 0.5*(mu[2:]-mu[:-2]); wmu[0]=0.5*(mu[1]-mu[0]); wmu[-1]=0.5*(mu[-1]-mu[-2]) + pw = fp * wmu; s = pw.sum() + if s <= 0: + M[gi[i], gi[i]] += src[i]; continue + go = np.clip(np.searchsorted(edges, E * (A*A + 2*A*mu + 1.0) / a1) - 1, 0, G - 1) + np.add.at(M[gi[i]], go, src[i] * pw / s) + continue + for prod in r.products: # inelastic: all neutron products + if prod.particle != 'neutron': + continue + try: + yld = np.atleast_1d(np.asarray(prod.yield_(grid), float)) + if yld.size == 1: + yld = np.full(len(grid), float(yld[0])) + except Exception: + yld = np.ones(len(grid)) + dists = prod.distribution; nd = len(dists) + appl = getattr(prod, 'applicability', None) + for k, dist in enumerate(dists): + if nd > 1 and appl and k < len(appl): + try: + app = np.atleast_1d(np.asarray(appl[k](grid), float)) + if app.size == 1: + app = np.full(len(grid), float(app[0])) + except Exception: + app = np.ones(len(grid)) + else: + app = np.ones(len(grid)) + ed = getattr(dist, 'energy', None) + if isinstance(ed, openmc.data.LevelInelastic): # discrete level + thr, mr = float(ed.threshold), float(ed.mass_ratio) + for i in nz: + E = grid[i] + if E <= thr: + continue + ecm = mr * (E - thr) + if ecm <= 0: + continue + base = ecm + E / a1; amp = 2.0 * np.sqrt(E * ecm) / (A + 1.0) + of = overlap(base - amp, base + amp); wgt = src[i] * yld[i] * app[i] + if of is None: + M[gi[i], min(max(np.searchsorted(edges, base) - 1, 0), G - 1)] += wgt + else: + M[gi[i]] += wgt * of + continue + td = _tabdist(dist) # continuum / (n,xn) + if td is None: + continue + Ea, GF = _unitbase_anchors(td, edges) + ein0 = float(np.asarray(td.energy, float)[0]) + for i in nz: + E = grid[i] + if E < ein0: + continue + j = np.searchsorted(Ea, E) + if j <= 0: + gf = GF[0] + elif j >= len(Ea): + gf = GF[-1] + else: + fr = (E - Ea[j-1]) / (Ea[j] - Ea[j-1]); gf = (1 - fr) * GF[j-1] + fr * GF[j] + M[gi[i]] += src[i] * yld[i] * app[i] * gf + + for g in range(G - 1): # fold upscatter artifacts into the diagonal + up = M[g, g+1:].sum() + if up: + M[g, g] += up; M[g, g+1:] = 0.0 + M = M / np.clip(denom[:, None], 1e-30, None) + return M[::-1, ::-1] # -> OpenMC ordering (group 1 = high E) From 65a553a732adc5d6b5c12a427fe53c1171466ee5 Mon Sep 17 00:00:00 2001 From: shimwell Date: Sat, 27 Jun 2026 00:11:39 +0200 Subject: [PATCH 09/14] Add TC-P0 transport-correction moment to scatter_matrix (return_p1) scatter_matrix(..., return_p1=True) also returns the P1 outscatter moment Sigma_s1,g = sum_g' mu_lab(g->g') * sigma_s0(g->g'), accumulated from the real CM elastic angular distribution. A transport-corrected (TC-P0) library for a P0 solver such as random ray is then sigma_tr = sigma_t - Sigma_s1 with the in-group diagonal reduced by the same Sigma_s1 (random ray expects the correction in the library, not the solver; negative diagonals are expected and handled by its diagonal stabilization). Verified Sigma_s1>=0, sigma_tr<=sigma_t; the correction is large for forward-peaked elastic (tungsten ~50% of sigma_t at 14 MeV). Default behaviour unchanged (return_p1=False). --- openmc/mgxs/transport_free.py | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/openmc/mgxs/transport_free.py b/openmc/mgxs/transport_free.py index eb579919f13..4db66a2a9f4 100644 --- a/openmc/mgxs/transport_free.py +++ b/openmc/mgxs/transport_free.py @@ -319,7 +319,8 @@ def _unitbase_anchors(td, edges, per=4): return np.array(Ea), np.array(GF) -def scatter_matrix(material, groups, temperature=294.0, cross_sections=None, source=None): +def scatter_matrix(material, groups, temperature=294.0, cross_sections=None, source=None, + return_p1=False): """Deterministic P0 group-to-group scattering matrix for one material. Returns the macroscopic Sigma_s,g->g' (1/cm) as an ``[G_in, G_out]`` array in @@ -360,6 +361,7 @@ def scatter_matrix(material, groups, temperature=294.0, cross_sections=None, sou wt = phi * dwid gi = np.clip(np.searchsorted(edges, grid, side='right') - 1, 0, G - 1) M = np.zeros((G, G)); denom = np.zeros(G); np.add.at(denom, gi, wt) + M1 = np.zeros((G, G)) if return_p1 else None # P1 (mu_lab-weighted) elastic outscatter def overlap(lo, hi): if hi <= lo: @@ -397,8 +399,12 @@ def overlap(lo, hi): pw = fp * wmu; s = pw.sum() if s <= 0: M[gi[i], gi[i]] += src[i]; continue - go = np.clip(np.searchsorted(edges, E * (A*A + 2*A*mu + 1.0) / a1) - 1, 0, G - 1) - np.add.at(M[gi[i]], go, src[i] * pw / s) + kin = A*A + 2*A*mu + 1.0 + go = np.clip(np.searchsorted(edges, E * kin / a1) - 1, 0, G - 1) + pwn = pw / s + np.add.at(M[gi[i]], go, src[i] * pwn) + if M1 is not None: # mu_lab = (1+A mu)/sqrt(A^2+2A mu+1) + np.add.at(M1[gi[i]], go, src[i] * pwn * (1.0 + A*mu) / np.sqrt(kin)) continue for prod in r.products: # inelastic: all neutron products if prod.particle != 'neutron': @@ -461,4 +467,10 @@ def overlap(lo, hi): if up: M[g, g] += up; M[g, g+1:] = 0.0 M = M / np.clip(denom[:, None], 1e-30, None) + if return_p1: + # P1 outscatter moment Sigma_s1,g per group, for a transport-corrected (TC-P0) + # library: emit sigma_tr = sigma_t - Sigma_s1 and subtract Sigma_s1 from the + # in-group diagonal. Random ray expects the correction applied here, not in the solver. + sigma_s1 = (M1.sum(1) / np.clip(denom, 1e-30, None))[::-1] + return M[::-1, ::-1], sigma_s1 return M[::-1, ::-1] # -> OpenMC ordering (group 1 = high E) From 7b7164a7eeba9fb08edd52e402e8b1ccb701c91e Mon Sep 17 00:00:00 2001 From: shimwell Date: Sat, 27 Jun 2026 03:26:17 +0200 Subject: [PATCH 10/14] Add free-gas thermal scattering kernel for light nuclides (scatter matrix) scatter_matrix(..., thermal=True) replaces the static target-at-rest elastic kernel with a free-gas (ideal-gas) energy-transfer kernel below 5 eV for light nuclides (A<=20), capturing thermal up-scatter and broadening that the static kernel misses. _freegas_gf() integrates over the target Maxwellian (isotropic CM); verified against analytic Wigner-Wilkins for A=1 to ~1-2%. Real thermal up-scatter is preserved (the upscatter-to-diagonal fold is skipped below 5 eV). Cuts the H2O/concrete scatter-matrix shape error ~2-3x (mean lethargy-gain vs material_wise: H2O 3.07->1.43, concrete 2.12->0.75); neutral for metals/heavy nuclides (A>20, unchanged). Uses free-gas (not bound S(alpha,beta)) -- matches OpenMC's default thermal treatment for materials without an S(a,b) table; a residual remains vs MC for hydrogenous moderators. --- openmc/mgxs/transport_free.py | 48 ++++++++++++++++++++++++++++++++--- 1 file changed, 45 insertions(+), 3 deletions(-) diff --git a/openmc/mgxs/transport_free.py b/openmc/mgxs/transport_free.py index 4db66a2a9f4..8b10d85d234 100644 --- a/openmc/mgxs/transport_free.py +++ b/openmc/mgxs/transport_free.py @@ -319,8 +319,34 @@ def _unitbase_anchors(td, edges, per=4): return np.array(Ea), np.array(GF) +def _freegas_gf(E, A, kT, edges, nv=28, nmu=14): + """Free-gas (ideal-gas) elastic energy-transfer group fractions at incident energy E. + + Quadrature over the target Maxwellian (speed v_t, cosine mu_t); for each target the + isotropic-CM elastic scatter gives the lab outgoing energy uniform over + [0.5(Vcm-w)^2, 0.5(Vcm+w)^2]. Captures thermal up-scatter and broadening that the + static (target-at-rest) kernel misses; reduces to Wigner-Wilkins for A=1 (verified + to ~1-2%). E and kT in eV; speeds in sqrt(eV) with m_n = 1. + """ + vn = np.sqrt(2.0 * E); vth = np.sqrt(2.0 * kT / max(A, 1e-9)) + vt = np.linspace(0.02 * vth, 6.0 * vth, nv) + wv = vt**2 * np.exp(-A * vt**2 / (2.0 * kT)) + mu = np.linspace(-1.0, 1.0, nmu) + VT = vt[:, None]; MU = mu[None, :] + vrel = np.sqrt(np.clip(vn*vn + VT*VT - 2*vn*VT*MU, 0, None)) + Vcm = np.sqrt(np.clip(vn*vn + A*A*VT*VT + 2*A*vn*VT*MU, 0, None)) / (A + 1.0) + w = (A / (A + 1.0)) * vrel + Elo = (0.5 * (Vcm - w)**2).ravel(); Ehi = (0.5 * (Vcm + w)**2).ravel() + wgt = ((wv[:, None]) * np.ones_like(mu)[None, :] * vrel).ravel() + rng = np.clip(Ehi - Elo, 1e-30, None); tot = wgt.sum() + if tot <= 0: + return None + cdf = (np.clip((edges[:, None] - Elo[None, :]) / rng[None, :], 0, 1) * wgt[None, :]).sum(1) / tot + return np.diff(cdf) + + def scatter_matrix(material, groups, temperature=294.0, cross_sections=None, source=None, - return_p1=False): + return_p1=False, thermal=True): """Deterministic P0 group-to-group scattering matrix for one material. Returns the macroscopic Sigma_s,g->g' (1/cm) as an ``[G_in, G_out]`` array in @@ -387,8 +413,20 @@ def overlap(lo, hi): ang = r.products[0].distribution[0].angle; aE = np.asarray(ang.energy) except Exception: ang = None + E_TH = 5.0; kT = 8.617e-5 * temperature # free-gas thermal for light nuclides + do_fg = thermal and A <= 20.0 + tgf = {} + if do_fg: + mid = np.sqrt(edges[:-1] * edges[1:]) + for g in np.where(mid < E_TH)[0]: + gf = _freegas_gf(mid[g], A, kT, edges) + if gf is not None: + tgf[int(g)] = gf for i in nz: - E = grid[i]; mu = fp = None + E = grid[i] + if do_fg and E < E_TH and gi[i] in tgf: # thermal: free-gas energy transfer + M[gi[i]] += src[i] * tgf[gi[i]]; continue + mu = fp = None if ang is not None: t = ang.mu[min(np.searchsorted(aE, E), len(aE) - 1)] if hasattr(t, 'x') and hasattr(t, 'p'): @@ -462,7 +500,11 @@ def overlap(lo, hi): fr = (E - Ea[j-1]) / (Ea[j] - Ea[j-1]); gf = (1 - fr) * GF[j-1] + fr * GF[j] M[gi[i]] += src[i] * yld[i] * app[i] * gf - for g in range(G - 1): # fold upscatter artifacts into the diagonal + # fold upscatter artifacts into the diagonal, but keep real thermal up-scatter (free-gas, below E_TH) + _midf = np.sqrt(edges[:-1] * edges[1:]) + for g in range(G - 1): + if thermal and _midf[g] < 5.0: + continue up = M[g, g+1:].sum() if up: M[g, g] += up; M[g, g+1:] = 0.0 From 90af72415d3bba33638c2e86a1b42a208a471714 Mon Sep 17 00:00:00 2001 From: shimwell Date: Sat, 27 Jun 2026 07:21:40 +0200 Subject: [PATCH 11/14] Use 400 kT free-gas threshold (match OpenMC default), not hardcoded 5 eV The free-gas thermal kernel now activates below 400*kT (= OpenMC's settings.free_gas_threshold default, ~10 eV at 294 K) instead of a hardcoded 5 eV, and the upscatter-fold uses the same cutoff. Temperature-scaled and consistent with the MC reference. When wired into convert_to_multigroup this should read settings.free_gas_threshold from the model rather than the 400 default. (Small numerical effect, since 5-10 eV is already well above thermal, but it is the correct, temperature-scaled value.) --- openmc/mgxs/transport_free.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/openmc/mgxs/transport_free.py b/openmc/mgxs/transport_free.py index 8b10d85d234..93660d29528 100644 --- a/openmc/mgxs/transport_free.py +++ b/openmc/mgxs/transport_free.py @@ -413,8 +413,9 @@ def overlap(lo, hi): ang = r.products[0].distribution[0].angle; aE = np.asarray(ang.energy) except Exception: ang = None - E_TH = 5.0; kT = 8.617e-5 * temperature # free-gas thermal for light nuclides - do_fg = thermal and A <= 20.0 + kT = 8.617e-5 * temperature # free-gas thermal for light nuclides + E_TH = 400.0 * kT # = OpenMC free_gas_threshold (default 400 kT); + do_fg = thermal and A <= 20.0 # when model-driven, read settings.free_gas_threshold tgf = {} if do_fg: mid = np.sqrt(edges[:-1] * edges[1:]) @@ -501,9 +502,9 @@ def overlap(lo, hi): M[gi[i]] += src[i] * yld[i] * app[i] * gf # fold upscatter artifacts into the diagonal, but keep real thermal up-scatter (free-gas, below E_TH) - _midf = np.sqrt(edges[:-1] * edges[1:]) + _midf = np.sqrt(edges[:-1] * edges[1:]); _eth = 400.0 * 8.617e-5 * temperature for g in range(G - 1): - if thermal and _midf[g] < 5.0: + if thermal and _midf[g] < _eth: continue up = M[g, g+1:].sum() if up: From 103bd8ea6c52b02676a3ea06472d4b352b1f2df6 Mon Sep 17 00:00:00 2001 From: shimwell Date: Sat, 27 Jun 2026 07:43:28 +0200 Subject: [PATCH 12/14] Fix elastic scatter binning: spread angular distribution via CDF, not deltas The elastic kernel deposited one delta per angular-grid mu point. For heavy nuclides (E_out ~ E, near-diagonal) this is fine, but for hydrogen the wide [0,E] down-scatter put every mu~-1 point (E_out~0) into the LOWEST group instead of spreading smoothly down the slowing-down range -- a ~14% spurious dump at the grid floor for fast incoming groups. Now the outgoing energy (monotonic in mu) is spread across outgoing groups via the angular CDF, integrated over group edges. Huge improvement for hydrogenous materials (CCFE-709 scatter-shape, mean lethargy-gain vs material_wise): H2O 1.37 -> 0.016, concrete 0.72 -> 0.035 -- both now BEAT stochastic_slab (0.041, 0.106). Metals unchanged (already near-diagonal). The P1/TC-P0 moment is spread consistently (row-sum unchanged). --- openmc/mgxs/transport_free.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/openmc/mgxs/transport_free.py b/openmc/mgxs/transport_free.py index 93660d29528..557e615c034 100644 --- a/openmc/mgxs/transport_free.py +++ b/openmc/mgxs/transport_free.py @@ -434,16 +434,20 @@ def overlap(lo, hi): mu = np.asarray(t.x, float); fp = np.asarray(t.p, float) if mu is None or mu.size < 2: mu = np.linspace(-1, 1, 33); fp = np.full(33, 0.5) - wmu = np.empty_like(mu); wmu[1:-1] = 0.5*(mu[2:]-mu[:-2]); wmu[0]=0.5*(mu[1]-mu[0]); wmu[-1]=0.5*(mu[-1]-mu[-2]) - pw = fp * wmu; s = pw.sum() - if s <= 0: + if mu[0] > mu[-1]: + mu = mu[::-1]; fp = fp[::-1] + # E_out is monotonic in mu, so spread f(mu) SMOOTHLY across outgoing + # groups via the angular CDF -- not one delta per mu point, which piles + # hydrogen's wide (mu~-1 -> E_out~0) down-scatter into the lowest group. + eout = E * (A*A + 2*A*mu + 1.0) / a1 + cmu = np.concatenate(([0.0], np.cumsum(0.5*(fp[1:]+fp[:-1])*np.diff(mu)))) + if cmu[-1] <= 0: M[gi[i], gi[i]] += src[i]; continue - kin = A*A + 2*A*mu + 1.0 - go = np.clip(np.searchsorted(edges, E * kin / a1) - 1, 0, G - 1) - pwn = pw / s - np.add.at(M[gi[i]], go, src[i] * pwn) + M[gi[i]] += src[i] * np.diff(np.interp(edges, eout, cmu/cmu[-1], left=0.0, right=1.0)) if M1 is not None: # mu_lab = (1+A mu)/sqrt(A^2+2A mu+1) - np.add.at(M1[gi[i]], go, src[i] * pwn * (1.0 + A*mu) / np.sqrt(kin)) + mulab = (1.0 + A*mu) / np.sqrt(A*A + 2*A*mu + 1.0) + cm1 = np.concatenate(([0.0], np.cumsum(0.5*(fp[1:]*mulab[1:]+fp[:-1]*mulab[:-1])*np.diff(mu)))) + M1[gi[i]] += src[i] * np.diff(np.interp(edges, eout, cm1, left=0.0, right=cm1[-1])) / cmu[-1] continue for prod in r.products: # inelastic: all neutron products if prod.particle != 'neutron': From ddee0a93d188b7554526d60372f9837c5236c590 Mon Sep 17 00:00:00 2001 From: shimwell Date: Sat, 27 Jun 2026 09:12:45 +0200 Subject: [PATCH 13/14] Self-shield the elastic scatter matrix in the URR (Bondarenko f-factor) The scatter matrix used the dilute (infinite-dilution) elastic cross section in the unresolved resonance range, while the vector total was already URR-corrected. That left the scatter row-sum too high across the URR for nuclides with probability tables. Apply the same probability-table band average to the elastic channel (table column 2), giving a per-nuclide micro elastic self-shielding factor f_el(E)=/sigma_el,smooth, and multiply the elastic XS by it before building the transfer matrix -- the SCALE/AMPX approach of self-shielding the 2D elastic matrix via a Bondarenko f-factor. Effect (CCFE-709 scatter row-sum vs material_wise): tungsten 3.60 -> 1.95 (now beats stochastic_slab 2.34), Zircaloy 1.02 -> 0.94 (beats slab 1.00). Materials whose nuclides have no URR table (steel/Fe-56/SiC/Li4SiO4/H2O/concrete/He) are unchanged -- they are fully resolved, already self-shielded by phi=1/Sigma_t. --- openmc/mgxs/transport_free.py | 47 +++++++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/openmc/mgxs/transport_free.py b/openmc/mgxs/transport_free.py index 557e615c034..a9aee98a345 100644 --- a/openmc/mgxs/transport_free.py +++ b/openmc/mgxs/transport_free.py @@ -136,6 +136,50 @@ def _apply_urr(incs, dens, temp_str, grid, sigma_t_smooth, temperature): return delta +def _urr_elastic_factor(incs, dens, temp_str, grid, sigma_t_smooth, temperature): + """Per-nuclide micro elastic self-shielding factor f_el(E) = /sigma_el,smooth + in the unresolved range (1.0 elsewhere), from the same probability-table band average + as :func:`_apply_urr` but for the elastic channel (table column 2). Lets the scatter + matrix self-shield its elastic consistently with the URR-corrected vector total -- + otherwise the scatter row-sum is the dilute (too-high) elastic across the URR. + """ + fac = {} + for nuc, n in dens.items(): + inc = incs[nuc] + urr = getattr(inc, 'urr', None) + if not urr: + continue + cand = [t for t in urr if urr.get(t) is not None and t in inc.temperatures] + if not cand: + continue + ts = min(cand, key=lambda t: abs(float(t[:-1]) - temperature)) + pt = urr[ts] + if not getattr(pt, 'multiply_smooth', False): + continue + e = np.asarray(pt.energy, float); tab = np.asarray(pt.table, float) + idx = np.where((grid >= e[0]) & (grid <= e[-1]))[0] + if idx.size == 0: + continue + st = inc[1].xs[ts](grid); se = inc[2].xs[ts](grid) + f = np.ones_like(grid) + for gj in idx: + if se[gj] <= 0: + continue + j = int(np.clip(np.searchsorted(e, grid[gj]), 1, len(e) - 1)) + row = tab[j - 1] if abs(e[j - 1] - grid[gj]) <= abs(e[j] - grid[gj]) else tab[j] + p = np.diff(np.concatenate(([0.0], row[0]))) + sig0 = (sigma_t_smooth[gj] - n * st[gj]) / n + if sig0 < 0: + sig0 = 0.0 + wgt = p / (st[gj] * row[1] + sig0) + den = wgt.sum() + if den <= 0: + continue + f[gj] = float((wgt * (se[gj] * row[2])).sum() / den) / se[gj] + fac[nuc] = f + return fac + + def collapse_material(material, groups, temperature=294.0, cross_sections=None, source=None): """Transport-free macroscopic multigroup cross sections for one material. @@ -382,6 +426,7 @@ def scatter_matrix(material, groups, temperature=294.0, cross_sections=None, sou if source is not None: w = w + _source_pdf(source, grid) phi = w / np.clip(sigma_t, 1e-30, None) + fel = _urr_elastic_factor(incs, dens, temp_str, grid, sigma_t, temperature) # URR elastic self-shielding dwid = np.empty_like(grid) dwid[1:-1] = 0.5 * (grid[2:] - grid[:-2]); dwid[0] = 0.5*(grid[1]-grid[0]); dwid[-1] = 0.5*(grid[-1]-grid[-2]) wt = phi * dwid @@ -404,6 +449,8 @@ def overlap(lo, hi): xs = r.xs[ts](grid) except Exception: continue + if is_el and nuc in fel: # self-shield elastic in the URR + xs = xs * fel[nuc] nz = np.nonzero(xs > 0)[0] if nz.size == 0: continue From 50776f3f74daa56aaef34111c3373385da941afb Mon Sep 17 00:00:00 2001 From: shimwell Date: Mon, 29 Jun 2026 09:53:24 +0200 Subject: [PATCH 14/14] Wire transport_free into convert_to_multigroup + random-ray fusion benchmark Add method="transport_free" to Model.convert_to_multigroup, dispatching to a new _generate_transport_free_mgxs that deterministically builds the MGXS library from the openmc.mgxs.transport_free module (NR+URR self-shielded vector XS + deterministic P0 scatter matrix; optional TC-P0 via correction="P0"). No Monte Carlo, no transport solve; positive total XS in every group, so it feeds random ray with no fixups. examples/transport_free_random_ray/ adds an end-to-end RR-vs-CE fusion-shield benchmark (neutronics-workshop tokamak materials: 14.1 MeV plasma source -> W -> steel -> Li -> concrete) comparing random ray fed by transport_free vs stochastic_slab against a continuous-energy reference. Finding: at fine fusion groups (CCFE-709) transport_free matches stochastic_slab accuracy (parity, within a few %) and beats it at coarse groups; transport_free is strictly more robust (deterministic, noise-free, valid in every group, while slab needed a negative-XS fixup to run in random ray). --- examples/transport_free_random_ray/README.md | 71 ++++++ .../transport_free_random_ray/rr_analyze.py | 106 +++++++++ .../rr_fusion_bench.py | 205 ++++++++++++++++++ openmc/model/model.py | 80 ++++++- 4 files changed, 460 insertions(+), 2 deletions(-) create mode 100644 examples/transport_free_random_ray/README.md create mode 100644 examples/transport_free_random_ray/rr_analyze.py create mode 100644 examples/transport_free_random_ray/rr_fusion_bench.py diff --git a/examples/transport_free_random_ray/README.md b/examples/transport_free_random_ray/README.md new file mode 100644 index 00000000000..9cfa6b5b5c7 --- /dev/null +++ b/examples/transport_free_random_ray/README.md @@ -0,0 +1,71 @@ +# Transport-free MGXS → random ray: fusion-shield benchmark + +A decisive end-to-end test of the `transport_free` MGXS generation method (added to +`Model.convert_to_multigroup(method="transport_free")`): generate a multigroup library +deterministically (no Monte Carlo, no transport solve), feed it to OpenMC's **random ray** +solver, and compare the resulting flux against a **continuous-energy Monte Carlo reference** +on a deep-penetration fusion shield. The same is done with `stochastic_slab` (the existing +MC-based method) so the two can be compared head-to-head against the CE truth. + +## Geometry + +Spherical shield using the [neutronics-workshop](https://fusion-energy.github.io/neutronics-workshop) +tokamak materials and thicknesses: a 14.1 MeV plasma source in the centre, then a tungsten +armor layer, a steel structural layer, a lithium breeder blanket, and a concrete bioshield. + +![benchmark geometry](rr_benchmark_geometry.png) + +## Running + +```bash +python rr_fusion_bench.py ce # continuous-energy reference (Monte Carlo) +python rr_fusion_bench.py rr stochastic_slab # random ray fed by the slab MGXS +python rr_fusion_bench.py rr transport_free # random ray fed by the transport-free MGXS +python rr_analyze.py # volume-averaged comparison vs CE +``` + +Append a group-structure name (e.g. `CCFE-709`) to any command to change the energy mesh. + +## Metrics + +CE flux tallies are volume-*integrated* (track length); random ray +`volume_normalized_flux_tallies` are volume-*averaged*. `rr_analyze.py` converts both to a +common volume-averaged region spectrum using analytic shell volumes, then reports: + +- **spectrum-shape error** vs CE per region (each spectrum normalized to unit sum — + normalization-independent), +- **deep-attenuation error**: how accurately the flux suppression deep in the shield + (φ_region / φ_first-wall) matches CE. + +## Result + +**Robustness.** `transport_free` produces a library that is **valid in every group** +(positive total XS, no Monte Carlo noise), so it feeds random ray with no fixups. +`stochastic_slab` produced zero / negative total cross sections (unpopulated deep-thermal +groups + a transport-correction artifact in lithium) that random ray rejects; the harness +applies an identical positive-floor fixup to both for a fair comparison. The noise is also +visible directly: in the deep-thermal tail of the near-source regions `stochastic_slab` +swings over tens of decades of Monte Carlo noise while `transport_free` is smooth (see +`rr_spectra_*.png`). + +**Accuracy vs CE.** The verdict is group-structure dependent: + +- At a **coarse** structure (`VITAMIN-J-42`, 40 groups) `transport_free` is markedly more + accurate deep in the shield — the deep-attenuation error vs CE is ~41 % versus ~234 % for + `stochastic_slab`, because the slab's noisy/biased coarse-group cross sections badly + mis-predict the flux suppression. +- At a **fine fusion** structure (`CCFE-709`, 650 groups) the two methods **converge to + near-parity** (flux-weighted spectral error within a few % of each other in every region): + `stochastic_slab` keeps a small edge at the near-source first wall, `transport_free` is + slightly better at the deepest point. Finer groups make the within-group weighting nearly + irrelevant, so the methods agree — consistent with the collapse-metric studies on PR #113. + +So for the random ray workflow, `transport_free` **matches `stochastic_slab`'s accuracy at +the fine groups used for fusion shielding, beats it at coarse groups, and is strictly more +robust** (deterministic, noise-free, valid in every group). + +**Caveat.** The random ray solve here is moderately resolved (a handful of flat source +regions per layer); absolute errors vs CE (~20–40 % flux-weighted at `CCFE-709`) are +dominated by that spatial resolution and are common to both methods. A higher-resolution +random ray solve (more source regions, linear sources, more rays) would sharpen the absolute +accuracy verdict; the *relative* method comparison is unaffected. diff --git a/examples/transport_free_random_ray/rr_analyze.py b/examples/transport_free_random_ray/rr_analyze.py new file mode 100644 index 00000000000..df7b089de40 --- /dev/null +++ b/examples/transport_free_random_ray/rr_analyze.py @@ -0,0 +1,106 @@ +"""Consistent volume-averaged comparison of RR(slab) and RR(transport_free) vs the CE +reference, reading the three statepoints directly. CE flux is volume-INTEGRATED +(track-length); RR volume_normalized flux is volume-AVERAGED. Convert both to a common +volume-AVERAGED region spectrum using analytic shell volumes, then score.""" +import sys, glob, numpy as np, openmc +# normalize argv so rr_fusion_bench's import-time group parsing sees the group name +_grp = sys.argv[1] if len(sys.argv) > 1 else "VITAMIN-J-42" +sys.argv = [sys.argv[0], "ce", _grp] +import rr_fusion_bench as B + +EDGES, NG, LAYOUT = B.EDGES, B.NG, B.LAYOUT +GROUPS = B.GROUPS + +# rebuild geometry to recover per-cell radii/volumes and the region->cells map +def cell_geometry(): + r_prev = 0.0; cells = [] # list of (region, r_in, r_out) + for name, r_out, nsub in LAYOUT: + rs = np.linspace(r_prev, r_out, nsub + 1) + for k in range(nsub): + cells.append((name, rs[k], rs[k + 1])) + r_prev = r_out + vols = np.array([4/3*np.pi*(ro**3 - ri**3) for _, ri, ro in cells]) + regions = [c[0] for c in cells] + return regions, vols + +def per_cell_flux(sp_path): + sp = openmc.StatePoint(sp_path) + t = sp.get_tally(name="flux_spectrum") + return t.get_values(scores=["flux"]).reshape(-1, NG) # [cell, group] low->high E + +def region_avg(flux_cell, regions, vols, integrated): + """Volume-averaged region spectrum [NG], OpenMC ordering (group1=high E). + integrated=True: flux_cell is vol-integrated (CE) -> sum/V. False: vol-averaged (RR) -> V-weighted mean.""" + names = [n for n, _, _ in LAYOUT] + out = {} + for nm in names: + idx = [i for i, r in enumerate(regions) if r == nm] + V = vols[idx].sum() + if integrated: + spec = flux_cell[idx].sum(0) / V + else: + spec = (flux_cell[idx] * vols[idx][:, None]).sum(0) / V + out[nm] = spec[::-1] + return out + +def shape_err(a, b): + a = a / a.sum(); b = b / b.sum(); k = b > b.max() * 1e-6 + return 100 * np.mean(np.abs((a - b)[k] / b[k])) + +regions, vols = cell_geometry() +ce_sp = sorted(glob.glob("run_ce/statepoint.*.h5"))[-1] +sl_sp = sorted(glob.glob("run_rr_stochastic_slab/statepoint.*.h5"))[-1] +tf_sp = sorted(glob.glob("run_rr_transport_free/statepoint.*.h5"))[-1] +ce = region_avg(per_cell_flux(ce_sp), regions, vols, integrated=True) +sl = region_avg(per_cell_flux(sl_sp), regions, vols, integrated=False) +tf = region_avg(per_cell_flux(tf_sp), regions, vols, integrated=False) + +names = [n for n, _, _ in LAYOUT] +fw = "tungsten" +print(f"=== RR-vs-CE fusion shield, {GROUPS} ({NG} groups), correction=None ===") +print(f"statepoints: CE={ce_sp} slab={sl_sp} tf={tf_sp}\n") +print(f"{'region':9}| spectrum-shape err vs CE (%) | attenuation phi_region/phi_tungsten") +print(f"{'':9}| slab tf (winner) | CE slab tf") +slab_wins = tf_wins = 0 +for nm in names: + if nm == "plasma": + continue + se_sl, se_tf = shape_err(sl[nm], ce[nm]), shape_err(tf[nm], ce[nm]) + win = "tf" if se_tf < se_sl else "slab" + if nm in ("steel", "lithium", "concrete"): # deep regions (downstream of resonant metals) + tf_wins += se_tf < se_sl; slab_wins += se_tf >= se_sl + a_ce = ce[nm].sum() / ce[fw].sum() + a_sl = sl[nm].sum() / sl[fw].sum() + a_tf = tf[nm].sum() / tf[fw].sum() + print(f"{nm:9}| {se_sl:6.2f} {se_tf:6.2f} ({win:4}) | {a_ce:.3e} {a_sl:.3e} {a_tf:.3e}") + +# headline: attenuation accuracy (deep flux suppression vs CE) + deep spectrum +def atten_err(m): + return 100 * np.mean([abs((m[nm].sum()/m[fw].sum()) / (ce[nm].sum()/ce[fw].sum()) - 1) + for nm in ("steel", "lithium", "concrete")]) +print(f"\nDeep-attenuation error vs CE (mean over steel/Li/concrete): " + f"slab {atten_err(sl):.1f}% transport_free {atten_err(tf):.1f}%") +deep = "concrete" +print(f"DEEP ({deep}) spectrum-shape error vs CE: slab {shape_err(sl[deep], ce[deep]):.2f}% " + f"tf {shape_err(tf[deep], ce[deep]):.2f}%") +print(f"\nDeep-region (steel/Li/concrete) spectrum-shape wins: transport_free {tf_wins}/3, slab {slab_wins}/3") + +# spectra plot: normalized flux/group, CE vs RR(slab) vs RR(transport_free) +try: + import matplotlib; matplotlib.use("Agg"); import matplotlib.pyplot as plt + mid = np.sqrt(EDGES[:-1] * EDGES[1:])[::-1] + fig, axs = plt.subplots(1, 2, figsize=(13, 5)) + for ax, nm in zip(axs, ("tungsten", "concrete")): + for lab, d, c, ls, lw in [("CE ref", ce, "k", "-", 1.6), + ("RR slab", sl, "tab:orange", "--", 1.2), + ("RR transport_free", tf, "tab:blue", "--", 1.2)]: + y = d[nm] / d[nm].sum() + ax.step(mid, y, where="mid", label=lab, color=c, ls=ls, lw=lw) + ax.set_xscale("log"); ax.set_yscale("log") + ax.set_title(f"{nm} region - normalized flux/group ({GROUPS})") + ax.set_xlabel("E [eV]"); ax.set_ylabel("normalized flux") + ax.legend(fontsize=8); ax.grid(alpha=0.3, which="both") + fig.tight_layout(); fig.savefig(f"rr_spectra_{GROUPS}.png", dpi=130, bbox_inches="tight") + print(f"saved rr_spectra_{GROUPS}.png") +except Exception as e: + print("plot skipped:", e) diff --git a/examples/transport_free_random_ray/rr_fusion_bench.py b/examples/transport_free_random_ray/rr_fusion_bench.py new file mode 100644 index 00000000000..2caac1ba627 --- /dev/null +++ b/examples/transport_free_random_ray/rr_fusion_bench.py @@ -0,0 +1,205 @@ +"""RR-vs-CE fusion-shield benchmark: does transport_free MGXS beat stochastic_slab when +fed to OpenMC's random ray solver, judged against a continuous-energy reference? + +Geometry: spherical fusion shield using the neutronics-workshop tokamak materials/thicknesses +(14.1 MeV plasma source -> tungsten armor -> steel structure -> lithium blanket -> concrete +bioshield). Deep-penetration: the deep-region spectrum depends on the upstream metals' MGXS. + +Metrics (normalization-independent, so CE-vs-RR is apples-to-apples): + - deep-region group flux SHAPE (each spectrum normalized to unit sum) + - radial attenuation profile (region flux / first-wall flux) + +usage: rr_fusion_bench.py {ce | rr | compare} [GROUPS] +""" +import sys, warnings, numpy as np, openmc +warnings.simplefilter("ignore") + +_gpos = 3 if (len(sys.argv) > 1 and sys.argv[1] == "rr") else 2 +GROUPS = sys.argv[_gpos] if len(sys.argv) > _gpos else "VITAMIN-J-42" +_full = np.asarray(openmc.mgxs.GROUP_STRUCTURES[GROUPS]); EDGES = _full[_full <= 2.0e7] +GE = openmc.mgxs.EnergyGroups(EDGES); NG = GE.num_groups +SRC_E = 14.1e6 + +# ---- materials (workshop compositions) ------------------------------------------------- +def materials(): + plasma = openmc.Material(name="plasma") # near-void DT-ish gas, source region + plasma.add_element("H", 1.0); plasma.set_density("g/cm3", 1e-5) + w = openmc.Material(name="tungsten"); w.add_element("W", 1.0); w.set_density("g/cm3", 19.3) + steel = openmc.Material(name="steel") # SS316-ish + for el, f in [("Fe", 0.66), ("Cr", 0.17), ("Ni", 0.12), ("Mo", 0.025), ("Mn", 0.02)]: + steel.add_element(el, f, "wo") + steel.set_density("g/cm3", 7.93) + li = openmc.Material(name="lithium"); li.add_element("Li", 1.0); li.set_density("g/cm3", 0.534) + conc = openmc.Material(name="concrete") + for el, f in [("O", 0.52), ("Si", 0.325), ("Ca", 0.06), ("Al", 0.033), + ("Fe", 0.014), ("H", 0.01), ("Na", 0.017), ("Mg", 0.002), ("K", 0.019)]: + conc.add_element(el, f, "wo") + conc.set_density("g/cm3", 2.3) + return [plasma, w, steel, li, conc] + +# region outer radii [cm]; sub = radial subdivisions (FSRs for RR spatial resolution) +LAYOUT = [("plasma", 30.0, 2), ("tungsten", 50.0, 3), ("steel", 70.0, 3), + ("lithium", 120.0, 5), ("concrete", 170.0, 5)] + +def build_model(energy_mode="continuous-energy"): + mats = {m.name: m for m in materials()} + surfs = []; r_prev = 0.0; cells = []; region_cells = {} + for name, r_out, nsub in LAYOUT: + rs = np.linspace(r_prev, r_out, nsub + 1)[1:] + cl = [] + for k, r in enumerate(rs): + s = openmc.Sphere(r=r) + inner = -surfs[-1] if surfs else None + cell = openmc.Cell(fill=mats[name], region=(-s & +surfs[-1]) if surfs else -s) + surfs.append(s); cells.append(cell); cl.append(cell) + region_cells[name] = cl + r_prev = r_out + surfs[-1].boundary_type = "vacuum" + geom = openmc.Geometry(cells) + st = openmc.Settings(); st.run_mode = "fixed source"; st.energy_mode = energy_mode + st.batches = 100 if energy_mode == "continuous-energy" else 100 + st.particles = 200000 if energy_mode == "continuous-energy" else 2000 + # physical source: uniform over plasma region, 14.1 MeV (CE) / top group (MG) + if energy_mode == "continuous-energy": + e_dist = openmc.stats.Discrete([SRC_E], [1.0]) + else: + mid = np.sqrt(EDGES[:-1] * EDGES[1:]); g = int(np.argmin(np.abs(mid - SRC_E))) + e_dist = openmc.stats.Discrete([mid[g]], [1.0]) + space = openmc.stats.Point() # placeholder; constrained below for RR + st.source = openmc.IndependentSource(energy=e_dist, angle=openmc.stats.Isotropic()) + model = openmc.Model(geom, openmc.Materials(list(mats.values())), st) + return model, region_cells + +def add_tally(model, region_cells): + """Flux per group in each region (volume-summed over its sub-cells).""" + t = openmc.Tally(name="flux_spectrum") + cf = openmc.CellFilter([c for cl in region_cells.values() for c in cl]) + ef = openmc.EnergyFilter(EDGES) + t.filters = [cf, ef]; t.scores = ["flux"] + model.tallies = openmc.Tallies([t]) + return [c.id for cl in region_cells.values() for c in cl] + +def region_flux(sp_path, cell_ids, region_cells): + """Return {region: group_flux[NG]} summing the sub-cells of each region.""" + sp = openmc.StatePoint(sp_path) + t = sp.get_tally(name="flux_spectrum") + fl = t.get_values(scores=["flux"]).reshape(len(cell_ids), NG) # [cell, group] low->high E + id_index = {cid: i for i, cid in enumerate(cell_ids)} + out = {} + for name, cl in region_cells.items(): + idx = [id_index[c.id] for c in cl] + out[name] = fl[idx].sum(0)[::-1] # -> high E = group 1 (OpenMC ordering) + return out + +# ---- stages ---------------------------------------------------------------------------- +def run_ce(): + model, rc = build_model("continuous-energy") + cell_ids = add_tally(model, rc) + sp = model.run(cwd="run_ce") + fl = region_flux(sp, cell_ids, rc) + np.savez(f"rr_ref_ce_{GROUPS}.npz", cell_ids=cell_ids, + **{k: v for k, v in fl.items()}) + print("CE reference done ->", f"rr_ref_ce_{GROUPS}.npz") + for k, v in fl.items(): + print(f" {k:9} integral flux = {v.sum():.4e}") + +def _sanitize(path, floor=1e-5): + """RR rejects zero/negative total XS. Clip total (and absorption) to a small floor + where needed; applied identically to every method so the comparison stays fair. + Returns the count of clipped (material, group) entries.""" + lib = openmc.MGXSLibrary.from_hdf5(path); n = 0 + for x in lib.xsdatas: + tot = np.array(x._total[0]) + bad = tot <= 0 + if bad.any(): + n += int(bad.sum()); tot[bad] = floor + x._total[0] = tot + ab = np.array(x._absorption[0]); ab[bad] = np.maximum(ab[bad], 0.0); x._absorption[0] = ab + if n: + lib.export_to_hdf5(path) + return n + +def run_rr(method): + model, rc = build_model("continuous-energy") # start CE, then convert + src_dist = model.settings.source[0].energy + model.convert_to_multigroup(method=method, groups=GE, + mgxs_path=f"rr_mgxs_{method}_{GROUPS}.h5", + overwrite_mgxs_library=True, correction=None, + source_energy=openmc.stats.Discrete([SRC_E], [1.0])) + nclip = _sanitize(f"rr_mgxs_{method}_{GROUPS}.h5") + print(f" [{method}] sanitized {nclip} zero/negative total-XS entries") + # constrain the physical source to the plasma material (RR fixed-source domain) + plasma_mat = [m for m in model.materials if m.name == "plasma"][0] + model.settings.source = [openmc.IndependentSource( + energy=src_dist, constraints={"domains": [plasma_mat]}, strength=1.0)] + cell_ids = add_tally(model, rc) + model.convert_to_random_ray() + model.settings.random_ray["volume_normalized_flux_tallies"] = True + model.settings.particles = 4000; model.settings.batches = 150; model.settings.inactive = 50 + sp = model.run(cwd=f"run_rr_{method}") + fl = region_flux(sp, cell_ids, rc) + np.savez(f"rr_res_{method}_{GROUPS}.npz", **{k: v for k, v in fl.items()}) + print(f"RR ({method}) done ->", f"rr_res_{method}_{GROUPS}.npz") + for k, v in fl.items(): + print(f" {k:9} integral flux = {v.sum():.4e}") + +def compare(): + ref = np.load(f"rr_ref_ce_{GROUPS}.npz") + regions = [n for n, _, _ in LAYOUT] + methods = ["stochastic_slab", "transport_free"] + res = {m: np.load(f"rr_res_{m}_{GROUPS}.npz") for m in methods} + def shape_err(a, b): # normalized-spectrum mean abs % error + a = a / a.sum(); b = b / b.sum(); k = b > b.max() * 1e-6 + return 100 * np.mean(np.abs((a - b)[k] / b[k])) + print(f"\n=== RR-vs-CE, {GROUPS} ({NG} groups) ===") + print(f"{'region':9}| spectrum shape err (%) vs CE | atten ratio (region/tungsten) CE / slab / tf") + fw = "tungsten" + for nm in regions: + if nm == "plasma": continue + ce = ref[nm] + line = f"{nm:9}| slab {shape_err(res['stochastic_slab'][nm], ce):6.2f} tf {shape_err(res['transport_free'][nm], ce):6.2f} |" + a_ce = ce.sum() / ref[fw].sum() + a_sl = res['stochastic_slab'][nm].sum() / res['stochastic_slab'][fw].sum() + a_tf = res['transport_free'][nm].sum() / res['transport_free'][fw].sum() + line += f" {a_ce:.3e} / {a_sl:.3e} / {a_tf:.3e}" + print(line) + # headline: deep concrete spectrum, who wins + deep = "concrete" + sl = shape_err(res['stochastic_slab'][deep], ref[deep]) + tf = shape_err(res['transport_free'][deep], ref[deep]) + print(f"\nDEEP ({deep}) spectrum shape error vs CE: slab {sl:.2f}% transport_free {tf:.2f}% " + f"-> {'transport_free WINS' if tf < sl else 'stochastic_slab wins'}") + +def plot_geometry(): + import matplotlib; matplotlib.use("Agg"); import matplotlib.pyplot as plt + import matplotlib.patches as mp + model, rc = build_model("continuous-energy") + mats = {m.name: m for m in model.materials} + colmap = {"plasma": (255, 235, 150), "tungsten": (90, 90, 110), "steel": (140, 150, 165), + "lithium": (190, 225, 235), "concrete": (205, 180, 150)} + colors = {mats[n]: c for n, c in colmap.items()} + R = LAYOUT[-1][1] + fig, ax = plt.subplots(figsize=(7.2, 7.2), dpi=150) + model.plot(basis="xy", width=(2*R+40, 2*R+40), pixels=(1400, 1400), + colors=colors, color_by="material", axes=ax, legend=False) + ax.set_title("RR-vs-CE fusion shield benchmark (XY slice)\n" + "14.1 MeV plasma source -> W armor -> steel -> Li blanket -> concrete", fontsize=11) + ax.set_xlabel("x [cm]"); ax.set_ylabel("y [cm]") + r_prev = 0.0; handles = [] + for name, r_out, _ in LAYOUT: + handles.append(mp.Patch(facecolor=np.array(colmap[name])/255.0, edgecolor="k", + label=f"{name} ({r_prev:.0f}-{r_out:.0f} cm)")) + r_prev = r_out + ax.legend(handles=handles, loc="upper right", fontsize=8, framealpha=0.9) + ax.add_patch(plt.Circle((0, 0), LAYOUT[0][1], fill=False, ls="--", lw=1.2, ec="red")) + ax.plot(0, 0, marker="*", color="red", ms=14) + fig.tight_layout(); fig.savefig("rr_benchmark_geometry.png", bbox_inches="tight") + print("saved rr_benchmark_geometry.png") + +if __name__ == "__main__": + stage = sys.argv[1] + if stage == "ce": run_ce() + elif stage == "rr": run_rr(sys.argv[2] if len(sys.argv) > 2 else "transport_free") + elif stage == "compare": compare() + elif stage == "plot": plot_geometry() + else: print(__doc__) diff --git a/openmc/model/model.py b/openmc/model/model.py index d927b65ae64..009d00cb157 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -2703,8 +2703,14 @@ def convert_to_multigroup( Parameters ---------- - method : {"material_wise", "stochastic_slab", "infinite_medium"}, optional - Method to generate the MGXS. + method : {"material_wise", "stochastic_slab", "infinite_medium", \ + "transport_free"}, optional + Method to generate the MGXS. "transport_free" collapses each + material's pointwise data deterministically (no Monte Carlo, no + transport solve) against a narrow-resonance weighting flux with + resonance self-shielding, and builds a deterministic P0 scattering + matrix; it introduces no Monte Carlo noise and is valid in every + group. See :func:`openmc.mgxs.transport_free.collapse_material`. groups : openmc.mgxs.EnergyGroups, str, or sequence of float, optional Energy group structure for the MGXS. Can be an :class:`openmc.mgxs.EnergyGroups` object, a string name of a @@ -2793,6 +2799,9 @@ def convert_to_multigroup( self._generate_stochastic_slab_mgxs( groups, nparticles, mgxs_path, correction, tmpdir, source_energy, temperatures, temperature_settings) + elif method == "transport_free": + self._generate_transport_free_mgxs( + groups, mgxs_path, correction, source_energy, temperatures) else: raise ValueError( f'MGXS generation method "{method}" not recognized') @@ -2809,6 +2818,73 @@ def convert_to_multigroup( self.settings.energy_mode = 'multi-group' + def _generate_transport_free_mgxs( + self, groups, mgxs_path, correction, source_energy, temperatures + ): + """Deterministically generate an MGXS library (no Monte Carlo, no transport solve). + + Each material's pointwise cross sections are collapsed against a narrow-resonance + (1/E + source) weighting flux with always-on resonance self-shielding (resolved + range by 1/(sigma_t + sigma_0); unresolved range via NJOY probability tables), + together with a deterministic P0 group-to-group scattering matrix (including a + free-gas thermal kernel for light nuclides). Unlike "material_wise" and + "stochastic_slab", this introduces no Monte Carlo noise and yields a positive + total cross section in every group. See :mod:`openmc.mgxs.transport_free`. + + Parameters + ---------- + groups : openmc.mgxs.EnergyGroups + Energy group structure for the MGXS. + mgxs_path : PathLike + Filename for the MGXS HDF5 file. + correction : str or None + "P0" applies the transport-corrected outscatter P0 correction + (sigma_tr = sigma_t - Sigma_s1, with Sigma_s1 removed from the in-group + diagonal); None emits the full P0 matrix. + source_energy : openmc.stats.Univariate or None + Weighting-flux source spectrum. If None, the model's source energy + distribution is used; if the model has no source, a pure 1/E weight is used. + temperatures : Sequence[float] or None + Temperatures to generate MGXS at. Defaults to room temperature. + """ + from openmc.mgxs.transport_free import collapse_material, scatter_matrix + + # Weighting source spectrum: explicit argument, else the model's source energy + # distribution, else None (pure 1/E slowing-down weighting). + src = source_energy + if src is None and self.settings.source: + s0 = self.settings.source + s0 = s0[0] if isinstance(s0, (list, tuple)) else s0 + src = getattr(s0, 'energy', None) + + temps = list(temperatures) if temperatures else [294.0] + transport_correct = (correction == "P0") + diag = np.arange(groups.num_groups) + + mgxs_lib = openmc.MGXSLibrary(energy_groups=groups) + for material in self.materials: + xsd = openmc.XSdata(material.name, groups, temperatures=temps) + xsd.order = 0 + for temperature in temps: + coll = collapse_material(material, groups, temperature=temperature, + source=src) + total = np.asarray(coll['total'], float) + scat = scatter_matrix(material, groups, temperature=temperature, + source=src, return_p1=transport_correct) + if transport_correct: + matrix, sigma_s1 = scat + matrix = np.array(matrix, float) + total = total - sigma_s1 + matrix[diag, diag] -= sigma_s1 + else: + matrix = np.asarray(scat, float) + xsd.set_total(total, temperature=temperature) + xsd.set_absorption(np.asarray(coll['absorption'], float), + temperature=temperature) + xsd.set_scatter_matrix(matrix[:, :, np.newaxis], temperature=temperature) + mgxs_lib.add_xsdata(xsd) + mgxs_lib.export_to_hdf5(mgxs_path) + def convert_to_random_ray(self): """Convert a multigroup model to use random ray.