Production deterministic-transport MGXS (switch-free, real geometry + thermal in-kernel)

FINDING_production.md

@@ -0,0 +1,92 @@
+# Production deterministic-transport (#117): first validation
+
+One full-stack 1D transport solve (real layer order) -> per-region flux -> collapse each
+material with the CORRECT kernel scatter XS (incl. multiplicity) and weighting = transport
+above the free-gas threshold (400 kT, universal constant) + NR below.
+
+## 10-material result (CCFE-709, vs converged stochastic_slab)
+| material | TOTAL det/slab | ROWSUM det/slab | SHAPE det/slab |
+|---|---|---|---|
+| **Fe-56** | **0.32 / 1.41** | **0.40 / 1.47** | 0.010/0.009 |
+| **Zircaloy** | **0.56 / 0.91** | **0.92 / 1.00** | 0.014/0.013 |
+| Helium | 0.04 / 0.04 | **0.05 / 9.66** | 0.116/0.523 |
+| concrete/H2O/Li4SiO4 | win/tie | slightly worse (<0.4%) | win/tie |
+| tungsten | 3.79 / 2.57 | 5.26 / 2.34 | 0.004/0.004 |
+| steel / CuCrZr | worse | worse | ~tie |
+
+**Wins: TOTAL 7/10, scatter ROWSUM 3/10, SHAPE 3/10.**
+
+## Headline + the catch
+- **Fe-56 cracked on BOTH total (0.32) and scatter (0.40)** -- the resonant metal NR failed.
+  The deep / cross-talk materials win; this is the whole point.
+- **Near-source metals (tungsten/steel/CuCrZr) regress vs NR.** Checked vs BOTH the sphere
+  reference AND the geometry-consistent slab1d reference: tungsten bad vs both (3.79 / 4.20)
+  -> this is NOT a sphere-vs-slab artifact. The transport's *innermost-shell* (source-adjacent)
+  flux is genuinely wrong -- almost certainly the **incident-boundary source treatment** vs the
+  references' volumetric source.
+
+## Implication for "switch-free"
+"Transport everywhere" does NOT beat NR everywhere -- it trades the near-source metals for Fe-56.
+Two ways forward:
+1. **Fix the innermost-shell source treatment** (volumetric isotropic source in the source cell
+   instead of an incident boundary) -- if that fixes tungsten, switch-free is viable.
+2. **Physically-motivated division (auto-detected, not arbitrary):** NR for the source-adjacent
+   region (where the flux IS the source spectrum and NR is exact) + transport for the deeper
+   materials (cross-talk). The transport flux itself signals the regime (source-like vs degraded).
+
+For deep-penetration shielding (the target use case) the deep materials dominate, so the method
+already delivers where it matters; the near-source first wall is where NR is fine anyway.
+
+## Update: volumetric source fix tried -> no effect (diagnostic)
+Replaced the incident-boundary source with a volumetric isotropic source in the source
+cell. Result is **byte-identical** (tungsten still 3.79/5.26). This is informative: the
+collapse uses the per-region SPECTRAL SHAPE, which is set by the slowing-down (energy)
+physics; the source's spatial/angular form doesn't change it (both inject the same 14 MeV
+muir spectrum). So the near-source error is **not** the source treatment, and (verified
+earlier vs slab1d) **not** geometry.
+
+The remaining cause is the transport's **slowing-down physics: the down-scatter source is
+ELASTIC-ONLY**. Tungsten has strong inelastic scattering, so its source-region flux comes
+out too hard -> wrong collapse. (Fe-56, deep, is unaffected: its flux is set by the
+incoming spectrum from upstream, not its own inelastic.)
+
+## Update: inelastic + (n,2n)/(n,3n) down-scatter ADDED to the transport (`inel_source.py`)
+The transport source is now elastic (fine gather) **plus** inelastic + (n,xn). The inelastic
+outgoing is smooth, so it's built as a COARSE transfer matrix `M[out,in]` per material
+(incl. multiplicity) and applied to the flux in an **outer iteration** (elastic solve ->
+add inelastic source from that flux -> re-solve, x2). +71s, ~no accuracy cost from the
+coarse grid (inelastic spectra are broad).
+
+### Result: TOTAL 7/10 -> **9/10** (steel fixed)
+| material | TOTAL det/slab | ROWSUM det/slab | was (elastic-only) |
+|---|---|---|---|
+| **steel** | **0.25 / 0.48** | **0.42 / 0.51** | 0.70 / 0.89 -> now beats slab on both |
+| **Fe-56** | **0.44 / 1.41** | **0.51 / 1.47** | 0.32 / 0.40 (slightly up, still crushes slab) |
+| **Zircaloy** | **0.57 / 0.91** | **0.94 / 1.00** | win |
+| CuCrZr | 1.06 / 1.02 | 1.47 / 1.05 | ~tie total |
+| tungsten | **3.79 / 2.57** | 5.26 / 2.34 | **UNCHANGED** |
+
+**Inelastic fixed steel (the win) but left tungsten exactly unchanged.**
+
+### Why tungsten is immune (it is NOT the inelastic, and NOT a build bug)
+Confirmed `M_inel[tungsten]` builds correctly (sum 17.5, all 5 W isotopes' MT51-91 + (n,2n)
+present) and is applied to the tungsten cells -- yet the collapse doesn't move. So tungsten's
+error is **source-region geometry**, not slowing-down physics:
+- In `material_wise` (the truth) tungsten is mixed through the **whole domain** and sees the
+  fully built-up equilibrium (~1/E) slowing-down spectrum.
+- In the slab tungsten sits **only at the front 6 cm**, so its flux is the 14 MeV source
+  lightly self-slowed -- **too hard** vs the equilibrium. The deep materials don't have this
+  problem (they sit where a slowed spectrum is physical).
+- **NR assumes exactly that equilibrium 1/E**, so NR is *better* for the front wall (1.96 vs
+  transport 3.79). This is intrinsic to placing the front material only at the front.
+
+## Conclusion for #117
+- Inelastic down-scatter **completes the transport's slowing-down physics**: steel fixed,
+  Fe-56/Zircaloy/deep materials all win -> **TOTAL 9/10, the deep-penetration goal delivered.**
+- The lone holdout is **tungsten, the first wall** -- a geometry/placement effect (front
+  material under-sees the equilibrium spectrum), not a physics gap. **NR is the right tool
+  there** (1.96), and the first wall is exactly where NR is valid.
+- Clean switch-free rule that needs no arbitrary threshold: **NR for the source layer**
+  (flux == source spectrum, NR exact) **+ transport for everything downstream**. Equivalent
+  to "use NR where there is no upstream material to slow neutrons down," which the geometry
+  itself defines.

prototype/det_prod.py

@@ -0,0 +1,117 @@
+"""Production deterministic-transport MGXS (switch-free) prototype.
+One full-stack 1D transport solve on the real layered geometry -> per-region local flux.
+Collapse each material against its region flux: TOTAL and SCATTER MATRIX (correct kernel
+scatter XS incl. multiplicity), with the weighting flux = transport above the fixed
+free-gas threshold (400 kT) and NR below it (NR exact in thermal; universal constant, not
+a per-problem switch). Validate all 10 materials vs material_wise + converged slab."""
+import sys, time, numpy as np, openmc, openmc.data
+from openmc.mgxs.transport_free import _macroscopic, _apply_urr, _source_pdf, _nearest_temperature
+from scatter_det import scatter_matrix
+from inel_source import coarse_grid, build_inel_transfer, inel_source_fine
+from mats import materials
+_trapz = getattr(np, 'trapezoid', None) or np.trapz
+GS = sys.argv[1] if len(sys.argv) > 1 else "CCFE-709"
+SRC = openmc.stats.muir(e0=14.06e6, m_rat=5.0, kt=20000.0)
+datalib = openmc.data.DataLibrary.from_xml(openmc.config['cross_sections'])
+_ms = materials(); mats = {m.name: m for m in _ms}; ORDER = [m.name for m in _ms]
+E = np.asarray(openmc.mgxs.GROUP_STRUCTURES[GS]); E = E[E <= 2e7]; edges = E; G = len(E)-1
+GE = openmc.mgxs.EnergyGroups(E); mid = np.sqrt(edges[:-1]*edges[1:])[::-1]; Ulg = -np.log(mid)
+E_TH = 400.0 * 8.617e-5 * 294.0                                  # free-gas threshold (universal constant)
+allnuc = set()
+for m in _ms: allnuc |= set(m.get_nuclide_atom_densities())
+INC = {}; TS = {}
+for nuc in allnuc:
+    inc = openmc.data.IncidentNeutron.from_hdf5(datalib.get_by_material(nuc, data_type='neutron')['path'])
+    INC[nuc] = inc; TS[nuc] = _nearest_temperature(inc, 294.0)
+grids = [edges] + [np.asarray(INC[n].energy[TS[n]]) for n in allnuc]
+grid = np.unique(np.concatenate(grids)); grid = grid[(grid >= edges[0]) & (grid <= edges[-1])]
+grid = np.unique(np.concatenate([grid, np.concatenate([np.geomspace(edges[g], edges[g+1], 9)[1:-1] for g in range(G)])]))
+u = np.log(grid); keep = [0]; last = u[0]
+for j in range(1, len(u)):
+    if u[j]-last >= 4e-4 or j == len(u)-1: keep.append(j); last = u[j]
+grid = grid[keep]; N = len(grid)
+dE = np.empty(N); dE[1:-1] = 0.5*(grid[2:]-grid[:-2]); dE[0] = 0.5*(grid[1]-grid[0]); dE[-1] = 0.5*(grid[-1]-grid[-2])
+inc_src = _source_pdf(SRC, grid)
+def mdat(nm):
+    d = mats[nm].get_nuclide_atom_densities()
+    st = _macroscopic({k: INC[k] for k in d}, d, TS, grid, 1); st = st + _apply_urr({k: INC[k] for k in d}, d, TS, grid, st, 294.0)[1]
+    kn = []
+    for nuc, n in d.items():
+        A = INC[nuc].atomic_weight_ratio; a = ((A-1)/(A+1))**2
+        if a >= 1-1e-9: continue
+        try: ss = n*INC[nuc][2].xs[TS[nuc]](grid)
+        except Exception: continue
+        if np.any(ss > 0): kn.append((ss/((1-a)*np.clip(grid, 1e-30, None))*dE, np.searchsorted(grid, grid/a, side='right')))
+    return st, kn
+MD = {nm: mdat(nm) for nm in ORDER}
+LAYERS = [(ORDER[j], 6.0 if j == 0 else 4.0, 6 if j == 0 else 4) for j in range(len(ORDER))]
+cells = []
+for nm, th, nc in LAYERS:
+    for _ in range(nc): cells.append((nm, th/nc))
+ncell = len(cells); dx = np.array([c[1] for c in cells]); cmat = [c[0] for c in cells]
+ST = np.array([MD[nm][0] for nm in cmat]); lay_sl = []; s = 0
+for nm, th, nc in LAYERS: lay_sl.append((nm, slice(s, s+nc))); s += nc
+mu, wmu = np.polynomial.legendre.leggauss(8); pos = mu > 0; neg = mu < 0; mu_p = mu[pos]; mu_n = np.abs(mu[neg])
+cg, be = coarse_grid(grid[0], grid[-1])                          # coarse grid for the inelastic source
+MINEL = {}
+for nm in ORDER:
+    d = mats[nm].get_nuclide_atom_densities(); MINEL[nm] = build_inel_transfer({k: INC[k] for k in d}, d, TS, cg, be)
+def solve_transport(q_inel):
+    phi = np.zeros((ncell, N))
+    for i in range(N-1, -1, -1):
+        Q = q_inel[:, i].copy()                                  # inelastic source (precomputed from prev pass)
+        for nm, sl in lay_sl:
+            for cf, jh in MD[nm][1]:                             # elastic down-scatter (fine gather)
+                j = jh[i]
+                if j > i+1: Q[sl] += phi[sl, i+1:j] @ cf[i+1:j]
+        Q[0] += inc_src[i]/dx[0]                                 # volumetric isotropic source in source cell
+        sigt = ST[:, i]; pa = np.zeros((ncell, 8)); psn = np.zeros(mu_p.size)
+        for c in range(ncell):
+            tM = 2*mu_p/dx[c]; po = np.clip((Q[c]+psn*(tM-sigt[c]))/(tM+sigt[c]), 0, None); pa[c, pos] = 0.5*(psn+po); psn = po
+        psn = np.zeros(mu_n.size)
+        for c in range(ncell-1, -1, -1):
+            tM = 2*mu_n/dx[c]; po = np.clip((Q[c]+psn*(tM-sigt[c]))/(tM+sigt[c]), 0, None); pa[c, neg] = 0.5*(psn+po); psn = po
+        phi[:, i] = pa @ wmu
+    return phi
+t0 = time.time()
+phi = solve_transport(np.zeros((ncell, N)))                      # pass 0: elastic only
+for _ in range(2):                                              # outer iterations: add inelastic source
+    q_inel = np.zeros((ncell, N))
+    for nm, sl in lay_sl:
+        for c in range(sl.start, sl.stop):
+            q_inel[c] = inel_source_fine(MINEL[nm], phi[c], grid, dE, cg, be)
+    phi = solve_transport(q_inel)
+print(f"[{GS}] grid {N}, transport+inelastic {time.time()-t0:.0f}s", flush=True)
+def coll(p, sig):
+    out = np.zeros(G)
+    for g in range(G):
+        k = (grid >= edges[g]) & (grid <= edges[g+1])
+        if k.sum() < 2: continue
+        x, pp, sg = grid[k], p[k], sig[k]; dd = _trapz(pp, x); out[g] = _trapz(sg*pp, x)/dd if dd > 0 else 0
+    return out[::-1]
+def lib(tag, nm, kind):
+    L = openmc.MGXSLibrary.from_hdf5(f"scatref_{GS}_{tag}.h5"); x = [a for a in L.xsdatas if a.name.startswith(nm)][0]
+    return np.array(x._total[0]) if kind == 't' else np.array(x._scatter_matrix[0])[..., 0]
+def et(a, b): k = np.abs(b) > 1e-9; return 100*np.mean(np.abs((a-b)[k]/b[k]))
+def rs(M, Nr): r = Nr.sum(1) > 1e-3; return 100*np.mean(np.abs((M.sum(1)-Nr.sum(1))[r]/Nr.sum(1)[r]))
+def mlg(M): r = M.sum(1); return np.where(r > 0, (M@Ulg)/np.clip(r, 1e-30, None), 0)
+def shp(M, Nr): a, b = mlg(M), mlg(Nr); r = Nr.sum(1) > 1e-3; return float(np.mean(np.abs((a-b)[r])))
+ithr = int(np.argmin(np.abs(grid-E_TH)))
+wins = {'tot': 0, 'rs': 0, 'shp': 0}
+print(f"{'material':9}| TOTAL det/slab | ROWSUM det/slab | SHAPE det/slab   (slab=converged)", flush=True)
+for nm, sl in lay_sl:
+    Vt = dx[sl]; phir = (phi[sl]*Vt[:, None]).sum(0)/Vt.sum(); st_fe = MD[nm][0]
+    phinr = (1.0/np.clip(grid, 1e-11, None)+inc_src)/np.clip(st_fe, 1e-30, None)
+    sc = phir[ithr]/max(phinr[ithr], 1e-30)
+    hyb = np.where(grid > E_TH, phir, phinr*sc)                  # transport>thr, NR below (free-gas const)
+    Mtr = scatter_matrix(mats[nm], GE, source=SRC, phi_ext=(grid, hyb), grid_ext=grid)
+    tot = coll(hyb, st_fe)
+    mwt = lib("mw", nm, 't'); mws = lib("mw", nm, 's')
+    try: slt = lib("slab10x", nm, 't'); sls = lib("slab10x", nm, 's')
+    except Exception: slt = lib("slab", nm, 't'); sls = lib("slab", nm, 's')
+    et_d, et_s = et(tot, mwt), et(slt, mwt); rs_d, rs_s = rs(Mtr, mws), rs(sls, mws); sh_d, sh_s = shp(Mtr, mws), shp(sls, mws)
+    if et_d <= et_s*1.05: wins['tot'] += 1
+    if rs_d <= rs_s*1.05: wins['rs'] += 1
+    if sh_d <= sh_s*1.05: wins['shp'] += 1
+    print(f"{nm:9}|  {et_d:5.2f} {et_s:5.2f}  |  {rs_d:5.2f} {rs_s:5.2f}  |  {sh_d:.3f} {sh_s:.3f}", flush=True)
+print(f"\nwins (det <= 1.05*slab):  TOTAL {wins['tot']}/10   ROWSUM {wins['rs']}/10   SHAPE {wins['shp']}/10", flush=True)

prototype/inel_source.py

@@ -0,0 +1,60 @@
+"""Inelastic + (n,2n)/(n,3n) down-scatter source for the deterministic transport.
+Builds a COARSE transfer matrix M[out,in] (macroscopic, incl. multiplicity) per material
+-- inelastic outgoing is smooth so a coarse grid + interpolation is accurate -- and applies
+it to the (fine) flux to give the inelastic slowing-down source, added to the elastic source
+in an outer iteration."""
+import numpy as np
+from sd_flux import _outgoing
+
+def coarse_grid(emin, emax, per_decade=50):
+    nb = int(max(per_decade*np.log10(emax/emin), 40))
+    be = np.logspace(np.log10(emin), np.log10(emax), nb+1)
+    return np.sqrt(be[:-1]*be[1:]), be                          # centers, edges
+
+def build_inel_transfer(incs, dens, ts, cg, be):
+    """M[out,in] (1/cm) = macroscopic inelastic+(n,xn) transfer rate per unit flux, incl mult."""
+    nb = len(cg); M = np.zeros((nb, nb))
+    for nuc, dn in dens.items():
+        inc = incs[nuc]; t = ts[nuc]
+        for mt, r in inc.reactions.items():
+            if not ((51 <= mt <= 91) or mt in (16, 17)): continue
+            try: sig = dn*r.xs[t](cg)
+            except Exception: continue
+            if not np.any(sig > 0): continue
+            prods = [p for p in r.products if p.particle == 'neutron']
+            if not prods: continue
+            try:
+                mult = np.atleast_1d(np.asarray(prods[0].yield_(cg), float))
+                if mult.size == 1: mult = np.full(nb, float(mult[0]))
+            except Exception:
+                mult = np.ones(nb)
+            try: og = _outgoing(prods[0].distribution[0])
+            except Exception: og = None
+            if og is None: continue
+            if og[0] == 'delta':                                # discrete level
+                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' continuum / (n,xn)
+                ein, eos = og[1], og[2]
+                for k in range(nb):
+                    if sig[k] <= 0 or cg[k] < ein[0]: continue
+                    tb = eos[min(np.searchsorted(ein, cg[k]), len(eos)-1)]
+                    cx, cp = np.asarray(tb.x, float), np.asarray(tb.p, float)
+                    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
+    return M
+
+def inel_source_fine(M, phi_fine, grid, dE, cg, be):
+    """Given coarse transfer M and a fine flux phi_fine[N], return the inelastic source
+    density on the fine grid (same units as the elastic down-scatter source)."""
+    binidx = np.clip(np.searchsorted(be, grid)-1, 0, len(cg)-1)
+    Phi = np.zeros(len(cg)); np.add.at(Phi, binidx, phi_fine*dE)   # bin-integrated flux
+    src_bin = M @ Phi                                              # total rate into each out bin /vol
+    width = np.diff(be)
+    q_coarse = src_bin/np.clip(width, 1e-30, None)                 # rate density per bin
+    return q_coarse[binidx]                                        # piecewise-constant onto fine grid

prototype/mats.py

@@ -0,0 +1,22 @@
+"""Shared 10-material test matrix for transport-free MGXS validation.
+Spans heavy shield, structural, breeder, coolant and hydrogenous materials."""
+import openmc
+
+def materials():
+    out = []
+    w = openmc.Material(name="tungsten"); w.add_element("W", 1.0); w.set_density("g/cm3", 19.3); out.append(w)
+    s = openmc.Material(name="steel"); s.add_element("Fe",.70); s.add_element("Cr",.18); s.add_element("Ni",.12); s.set_density("g/cm3",7.9); out.append(s)
+    f = openmc.Material(name="Fe56"); f.add_nuclide("Fe56",1.0); f.set_density("g/cm3",7.87); out.append(f)
+    z = openmc.Material(name="CuCrZr"); z.add_element("Cu",0.9905,percent_type="wo"); z.add_element("Cr",0.008,percent_type="wo"); z.add_element("Zr",0.0015,percent_type="wo"); z.set_density("g/cm3",8.9); out.append(z)
+    zr = openmc.Material(name="Zircaloy"); zr.add_element("Zr",0.9823,percent_type="wo"); zr.add_element("Sn",0.0145,percent_type="wo"); zr.add_element("Fe",0.0021,percent_type="wo"); zr.add_element("Cr",0.001,percent_type="wo"); zr.set_density("g/cm3",6.56); out.append(zr)
+    sic = openmc.Material(name="SiC"); sic.add_elements_from_formula("SiC"); sic.set_density("g/cm3",3.21); out.append(sic)
+    lso = openmc.Material(name="Li4SiO4"); lso.add_elements_from_formula("Li4SiO4"); lso.set_density("g/cm3",2.40); out.append(lso)  # fusion breeder ceramic
+    cc = openmc.Material(name="concrete")
+    for el, wo in [("H",.01),("O",.529),("Si",.337),("Ca",.044),("Al",.034),("Fe",.014),("Na",.016),("K",.013),("Mg",.002),("C",.001)]:
+        cc.add_element(el, wo, percent_type="wo")
+    cc.set_density("g/cm3", 2.3); out.append(cc)
+    h2o = openmc.Material(name="H2O"); h2o.add_elements_from_formula("H2O"); h2o.set_density("g/cm3",1.0); out.append(h2o)
+    he = openmc.Material(name="Helium"); he.add_element("He",1.0); he.set_density("g/cm3",0.00311); out.append(he)  # 5 MPa, 500 C ideal-gas
+    return out
+
+NAMES = [m.name for m in materials()]