Add collapsed_self_shielding MGXS generation method - #118
Add collapsed_self_shielding MGXS generation method#118jon-proximafusion wants to merge 1 commit into
Conversation
Add method="collapsed_self_shielding" to Model.convert_to_multigroup: a deterministic, NJOY/FISPACT-style multigroup generator that collapses (group-averages) each material continuous-energy data against an assumed narrow-resonance weighting flux with always-on resonance self-shielding (resolved range plus unresolved-range probability tables), and builds a deterministic P0 group-to-group scattering matrix (with a free-gas thermal kernel for light nuclides). Optional transport correction via correction="P0". Unlike the existing material_wise / stochastic_slab / infinite_medium methods it uses no Monte Carlo and no transport solve, so it introduces no statistical noise and yields a positive total cross section in every group -- which lets it feed solvers such as random ray with no zero/negative-cross-section fixups. Implementation in openmc/mgxs/collapsed_self_shielding.py (collapse_material for the vector cross sections, scatter_matrix for the P0 matrix); wired into convert_to_multigroup via _generate_collapsed_self_shielding_mgxs.
Verification vs
|
| Material | Ref noise (total) | Ref noise (scatter) | Total XS (slab) | Total XS (css) | Scatter (slab) | Scatter (css) |
|---|---|---|---|---|---|---|
| concrete | 0.04 | 0.06 | 0.21 | 0.85 | 0.22 | 1.26 |
| graphite | 0.05 | 0.05 | 0.07 | 0.70 | 0.07 | 1.12 |
| lipb_he | 0.01 | 0.23 | 0.62 | 1.06 | 0.62 | 0.25 |
| plasma | 0.00 | 100.00 | 0.03 | 0.80 | 100.00 | 0.07 |
| steel | 0.10 | 0.17 | 0.23 | 3.30 | 0.27 | 3.10 |
| tungsten | 0.13 | 0.23 | 5.91 | 8.90 | 5.30 | 7.60 |
| water_steel | 0.02 | 0.03 | 0.08 | 0.96 | 0.08 | 1.07 |
(slab is stochastic_slab, css is collapsed_self_shielding.)
Reading the table
- Both multigroup methods agree with the well-converged reference to within a few %.
stochastic_slabis generally closest: it is itself Monte Carlo (so it shares the
reference's transport physics) and its mixed-material spectrum resembles the true degraded
in-place spectra deep in the shield.collapsed_self_shieldingtrails by low single-digit %, largest on the deep steel layer
(3.3 %) and the thin near-source tungsten (8.9 %, though tungsten is hard for slab too at 5.9 %).
This is the spatial-spectrum limitation: a single assumed (1/E + source) weighting flux cannot
match every material's true in-place spectrum. Finer groups shrink it for steel (5.5 % down to
3.3 % from VITAMIN-J-42 to VITAMIN-J-175) but do not close it, because it is a weighting
difference, not a group-resolution one.collapsed_self_shielding's advantages here are structural rather than on this metric: zero
Monte-Carlo noise (e.g. plasma scatter, slab 100 % vs css 0.07 %), a positive total cross
section in every group, and deterministic reproducibility, at competitive (low single-digit %)
accuracy.
Note that material_wise-as-truth inherently favors the other Monte-Carlo method; the decisive
test of whether these MGXS differences matter is downstream transport (e.g. a random-ray vs
continuous-energy benchmark).
verification script (VITAMIN-J-175)
"""Verify collapsed_self_shielding MGXS against material_wise (reference) and
stochastic_slab on a fusion shield, VITAMIN-J groups. material_wise is run with two
seeds so its seed-to-seed noise (the reference's own uncertainty) is reported alongside
the per-material total-XS and scatter-matrix-rowsum % differences."""
import numpy as np, openmc
GROUPS, NP = "VITAMIN-J-175", 30000
edges = np.array(openmc.mgxs.GROUP_STRUCTURES[GROUPS]); edges = edges[edges <= 2.0e7]
GE = openmc.mgxs.EnergyGroups(edges)
SRC = openmc.stats.muir(e0=14.06e6, m_rat=5.0, kt=20000.0) # 14 MeV DT, Gaussian (normal) broadened
def build(seed=1):
openmc.reset_auto_ids()
dt = openmc.Material(name="plasma"); dt.add_nuclide("H2", .5); dt.add_nuclide("H3", .5); dt.set_density("g/cm3", 1e-8)
w = openmc.Material(name="tungsten"); w.add_element("W", 1.); w.set_density("g/cm3", 19.3)
water = openmc.Material(); water.add_element("H", 2.); water.add_element("O", 1.); water.set_density("g/cm3", 1.)
steel = openmc.Material(name="steel")
for el, f in [("Fe", .70), ("Cr", .18), ("Ni", .12)]: steel.add_element(el, f, "wo")
steel.set_density("g/cm3", 7.93)
lipb = openmc.Material(); lipb.add_element("Li", .17); lipb.add_element("Pb", .83); lipb.set_density("g/cm3", 9.4)
he = openmc.Material(); he.add_element("He", 1.); he.set_density("g/cm3", 0.00311) # ~5 MPa, 500 C
gr = openmc.Material(name="graphite"); gr.add_element("C", 1.); gr.set_density("g/cm3", 1.7)
conc = openmc.Material(name="concrete")
for el, f in [("O", .52), ("Si", .325), ("Ca", .06), ("Al", .033), ("Fe", .014),
("H", .01), ("Na", .017), ("Mg", .002), ("K", .019)]: conc.add_element(el, f, "wo")
conc.set_density("g/cm3", 2.3)
ws = openmc.Material.mix_materials([water, steel], [.5, .5], "vo"); ws.name = "water_steel"
lph = openmc.Material.mix_materials([lipb, he], [.8, .2], "vo"); lph.name = "lipb_he"
layers = [(dt, 50), (w, 0.1), (ws, 3), (lph, 60), (gr, 5), (steel, 2), (None, 100), (conc, 200)]
surfs, cells, r = [], [], 0.0
for fill, t in layers:
r += t; s = openmc.Sphere(r=r)
cells.append(openmc.Cell(fill=fill, region=(-s & +surfs[-1]) if surfs else -s)); surfs.append(s)
surfs[-1].boundary_type = "vacuum"
st = openmc.Settings(); st.run_mode = "fixed source"; st.batches = 10; st.particles = NP; st.seed = seed
st.source = openmc.IndependentSource(space=openmc.stats.Point(), energy=SRC)
return openmc.Model(openmc.Geometry(cells), openmc.Materials([dt, w, ws, lph, gr, steel, conc]), st)
for method, seed, path in [("material_wise", 1, "mw.h5"), ("material_wise", 2, "mw2.h5"),
("stochastic_slab", 1, "slab.h5"), ("collapsed_self_shielding", 1, "css.h5")]:
build(seed).convert_to_multigroup(method=method, groups=GE, nparticles=NP,
mgxs_path=path, overwrite_mgxs_library=True)
def read(path):
L = openmc.MGXSLibrary.from_hdf5(path)
return {x.name: (np.array(x._total[0]), np.array(x._scatter_matrix[0])[..., 0].sum(1)) for x in L.xsdatas}
ref, ref2, slab, css = (read(p) for p in ["mw.h5", "mw2.h5", "slab.h5", "css.h5"])
def pd(a, b): # mean abs % difference where the reference is significant
b = np.asarray(b); k = b > b.max() * 1e-6
return 100 * np.mean(np.abs(np.asarray(a)[k] - b[k]) / b[k]) if k.any() else float("nan")
print(f"\n{GROUPS}: per-material % difference vs material_wise (mw 2-seed noise = reference uncertainty)")
print(f"{'material':12} | mw-noise tot/scat | total: slab css | scatter: slab css")
for n in ref:
print(f"{n:12} | {pd(ref2[n][0], ref[n][0]):6.2f} {pd(ref2[n][1], ref[n][1]):6.2f} |"
f" {pd(slab[n][0], ref[n][0]):7.2f} {pd(css[n][0], ref[n][0]):7.2f} |"
f" {pd(slab[n][1], ref[n][1]):7.2f} {pd(css[n][1], ref[n][1]):7.2f}")
Add a
collapsed_self_shieldingMGXS generation methodAdds a new option to
Model.convert_to_multigroup:It generates the multigroup library deterministically, in the style of NJOY (GROUPR) and
FISPACT: each material's continuous-energy data is collapsed (group-averaged) against an
assumed narrow-resonance weighting flux with always-on resonance self-shielding (resolved
range via
1/(Sigma_t + sigma_0)weighting, unresolved range via the NJOY probability tablesalready carried in
openmc.data). It also builds a deterministic P0 group-to-group scatteringmatrix (elastic with the real CM angular distribution, discrete inelastic levels,
unit-base-interpolated continuum and (n,xn) with multiplicity, plus a free-gas thermal kernel for
light nuclides). An optional transport correction is available via
correction="P0".Unlike the existing
material_wise,stochastic_slab, andinfinite_mediummethods, it usesno Monte Carlo and no transport solve.
Why a user might prefer it
This method is not pitched as more accurate than the Monte Carlo methods (see Accuracy below).
Its value is robustness, coverage, and usability:
seed dependence.
nparticles, and no convergence to check. The Monte Carlo methods make the user choose aparticle count and verify it converged; this method has no statistical knob. Its only
convergence is the energy grid, which is deterministic and monotonic.
transport corrections, negative) cross sections in deep or unpopulated groups. This method is
always positive everywhere, which is what lets it feed deterministic solvers such as random
ray or Sn with no zero or negative cross section fixups.
material_wisenotes in its owndocstring that materials far from the source may score no tallies and come out with zero cross
sections (pointing users to
stochastic_slab). This method produces a valid library for everymaterial from its composition, wherever it sits, and without
stochastic_slab's shuffled latticegeometry distortion.
energy tail groups, which are exactly the ones that drive shielding and dose. This method is
unaffected there.
VITAMIN-J-175 or CCFE-709 for fusion) come out clean, whereas Monte Carlo needs proportionally
more particles per group to suppress noise.
transportable geometry or converged source is required to produce a material's cross sections).
statistical tolerances.
Accuracy
On a deep, resonant metal heavy fusion shield, verified against a well converged
material_wisereference (2 seed noise at most 0.23 % everywhere), this method agrees to within a few percent and
is competitive with
stochastic_slab, but it does not beat it on this metric: it trails by lowsingle digit percent on the deep resonant metals (for example steel ~3 %), because a single
assumed weighting flux cannot reproduce each material's degraded in place spectrum (a spatial
spectrum limitation that finer groups do not remove). See the verification comment for the full per
material table and the script. Capturing those deep spectra accurately requires a deterministic
transport pre pass, which is a separate and more complex approach not included here.
Implementation
openmc/mgxs/collapsed_self_shielding.py:collapse_material(vector total / absorption /capture / fission) and
scatter_matrix(P0 matrix). Pure NumPy plusopenmc.data; no newdependencies.
Model._generate_collapsed_self_shielding_mgxswires it intoconvert_to_multigroup.Notes for review
(per nuclide multigroup data is never summed, since group cross sections are flux weighted
averages that do not add cleanly).
a comment rather than committed.