Real-geometry deterministic transport (1D spherical Sn) for MGXS weighting

FINDING_realgeom.md

@@ -0,0 +1,95 @@
+# Real-geometry deterministic transport: status
+
+Goal: solve the deterministic transport in the ACTUAL reference geometry (concentric
+spheres, central point source) so the near-source shells get the correct (harder) flux
+that the 1D-slab prototype (#114) over-softens.
+
+## Built + validated: 1D spherical Sn (sphere_sn.py)
+Weighted-diamond curvilinear Sn with the angular-redistribution recursion + central
+point source. **Validated against the analytic point-source-in-uniform-absorber flux
+phi(r)=S exp(-Sigma_t r)/(4 pi r^2): 0.9% mean, 1.6% max** -- the curvilinear physics
+(streaming + angular redistribution + central source) is correct.
+
+## Open issue: diamond-difference negativity on the resonant problem
+On the real material stack (resonances -> optically-thick cells + the alpha angular-
+redistribution term), the diamond scheme produces negative/oscillating fluxes:
+  steel TOTAL 7.37 (NR 0.57!), broken scatter SHAPE 0.76.
+Naive clamping (set psi>=0) breaks conservation and doubles the smooth-case flux (ratio
+2.0) -- it is NOT a valid fixup. A positivity-preserving CURVILINEAR scheme is required
+(step characteristic, or weighted-diamond with a proper negative-flux fixup that also
+treats the angular edges). This is the next engineering step.
+
+## Where it stands
+- Deep/cross-talk materials are already won by the 1D-SLAB transport (#114): Fe-56 0.32,
+  H2O/Li4SiO4/concrete/He all beat slab, deterministically.
+- The spherical solver is needed only to also fix the NEAR-SOURCE shells (tungsten,
+  steel) -- which is exactly where NR is already excellent (<=2%).
+- Pragmatic alternatives to the full positivity-preserving spherical Sn:
+  (a) hybrid: NR for near-source shells (hard, source-dominated) + slab transport for the
+      deep shells (cross-talk) -- the local flux itself signals the regime;
+  (b) ship #113 (NR) as the robust geometry-independent library and #114 (slab transport)
+      as the optional deep-material enhancement.
+
+## Update: positivity-preserving step-characteristic scheme added
+The curvature coefficients alpha_{m+1/2} are provably >= 0 (a tent: 0 -> peak -> 0),
+which makes STEP-CHARACTERISTIC differencing (upwind in space AND angle)
+UNCONDITIONALLY POSITIVE for spherical Sn. Implemented as scheme='step' (default) in
+sphere_sn.py. Validation (point source in absorber): step 1.2% mean (diamond 0.9%),
+min flux > 0 (vs diamond which can go negative).
+
+Re-run on the resonant stack (CCFE-709, total %err vs material_wise):
+| material | NR | diamond(broken) | STEP | slab |
+|---|---|---|---|---|
+| tungsten | 1.96 | 5.51 | **2.13** | 2.56 |
+| steel | 0.57 | 7.37 | 2.44 | 0.48 |
+| Fe-56 | 3.26 | (neg) | **1.08** | 1.39 |
+| CuCrZr | 0.85 | (neg) | 1.20 | 1.02 |
+
+Step fixes the negativity (no more catastrophic steel 7.37) and tungsten + Fe-56 now
+BEAT slab. BUT step is 1st-order diffusive, so steel/CuCrZr are worse than NR/slab.
+The diamond scheme is 2nd-order accurate but unstable on resonances. Neither is a clean
+win: the next step for a uniform spherical solver is **diamond with a proper negative-
+flux fixup** (set-to-zero + re-solve the cell, conserving), which keeps 2nd-order
+accuracy AND positivity.
+
+## Practical comparison of paths (CCFE-709)
+- The 1D-SLAB transport (#114) + thermal-NR fallback is the most accurate for the DEEP
+  materials: Fe-56 total 0.60 / scatter 0.69 (both beat slab) -- but slab geometry
+  over-softens the NEAR-SOURCE shells.
+- This spherical-STEP solver fixes near-source tungsten (2.13 < slab 2.56) but its
+  diffusivity hurts steel/CuCrZr.
+- A diamond+fixup spherical solver should get the best of both (accurate + positive +
+  correct geometry) -- the recommended future direction for a single uniform method.
+
+## Update 2: diamond + negative-flux fixup added (scheme='fixup', now default)
+Set-to-zero & re-solve: where a diamond cell's outgoing spatial/angular edge would go
+negative, clamp it to 0 and re-solve the cell (the curvature coeffs alpha>=0 guarantee
+termination positive). Validation (point source in absorber): **fixup 0.74% mean (best;
+diamond 0.89, step 1.17), positive flux.**
+
+BUT on the RESONANT stack (CCFE-709, total %err vs material_wise):
+| material | step | fixup | slab-transport(#114) | slab(MC) |
+|---|---|---|---|---|
+| tungsten | 2.13 | 2.47 | (5.51 slab-geom) | 2.56 |
+| steel    | 2.44 | 2.53 | 0.57(NR) | 0.48 |
+| Fe-56    | 1.08 | 2.06 | **0.60** | 1.39 |
+| CuCrZr   | 1.20 | 1.63 | 0.85(NR) | 1.02 |
+
+**fixup is WORSE than step on resonances** -- at the many resonance negativities the
+set-to-zero clamp accumulates more error than step's consistent upwinding. So for the
+resonant deep-penetration problem, neither spherical scheme beats the simpler, robust
+1D-SLAB transport (#114, diamond + edge clamp): Fe-56 0.60 (slab-transport) vs 1.08
+(spherical step). The spherical geometry's only clear gain is the near-source shell
+(tungsten 2.13 vs slab-method 2.56), where NR is already excellent anyway.
+
+## Bottom line / recommendation
+A correct-geometry spherical Sn does NOT outperform the robust slab transport on the
+deep resonant materials -- the spherical central-source + resonance negativity make it
+numerically harder, and the schemes that are positive (step/fixup) are too diffusive or
+clamp-lossy. The practical best method remains:
+  * #113 (NR) as the geometry-free default, +
+  * #114 (1D-slab transport + thermal-NR fallback) for the deep cross-talk materials
+    (Fe-56 total 0.60 / scatter 0.69, both beat slab),
+  * NR for the near-source shells (already <=2%).
+A uniform high-accuracy positive spherical solver would need more advanced numerics
+(characteristic/CN, or much finer mesh) -- diminishing returns vs the hybrid above.

det_realgeom.py

@@ -0,0 +1,97 @@
+"""Real-geometry deterministic transport: 1D SPHERICAL Sn through the actual reference
+geometry (concentric shells, central point source) -> per-shell local flux -> collapse
+each material. Matches the material_wise/slab reference geometry exactly. No Monte Carlo."""
+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 sphere_sn import solve_sphere
+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)
+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 not np.any(ss > 0): continue
+        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}
+# ---- spherical mesh: shell radii (W: 0-6, then 4cm shells), subdivided ----
+rad = [0.0, 6.0] + [6.0+4.0*i for i in range(1, len(ORDER))]   # outer radii: 0,6,10,...,42
+redge = [0.0]; cmat = []
+for j in range(len(ORDER)):
+    r0, r1 = rad[j], rad[j+1]; nc = 6 if j == 0 else 4
+    for c in range(nc): redge.append(r0+(r1-r0)*(c+1)/nc); cmat.append(ORDER[j])
+redge = np.array(redge); ncell = len(cmat)
+V = 4*np.pi/3*(redge[1:]**3 - redge[:-1]**3)
+ST = np.array([MD[nm][0] for nm in cmat])
+lay_sl = []; s = 0
+for j in range(len(ORDER)):
+    nc = 6 if j == 0 else 4; lay_sl.append((ORDER[j], slice(s, s+nc))); s += nc
+mu, wmu = np.polynomial.legendre.leggauss(8); alpha = np.zeros(9)
+for m in range(8): alpha[m+1] = alpha[m] - mu[m]*wmu[m]
+print(f"[{GS}] grid {N}, {ncell} spherical cells (R={redge[-1]:.0f}cm)", flush=True)
+# ---- energy sweep ----
+phi = np.zeros((ncell, N)); t0 = time.time()
+for i in range(N-1, -1, -1):
+    Svol = np.zeros(ncell)
+    for nm, sl in lay_sl:
+        for cf, jh in MD[nm][1]:
+            j = jh[i]
+            if j > i+1: Svol[sl] += phi[sl, i+1:j] @ cf[i+1:j]
+    Svol = Svol/np.clip(V, 1e-30, None)               # downscatter source per unit volume
+    Svol[0] += inc_src[i]/V[0]                          # central point source in innermost cell
+    phi[:, i] = solve_sphere(redge, ST[:, i], Svol, mu, wmu, alpha)
+print(f"  spherical transport solve {time.time()-t0:.0f}s", flush=True)
+# ---- collapse each material vs its shell flux (spherical volume average) ----
+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])))
+print(f"{'material':9}| TOTAL: NR / TR / slab | ROWSUM: NR / TR / slab | SHAPE: TR / slab", flush=True)
+for nm, sl in lay_sl:
+    Vt = V[sl]; phi_loc = (phi[sl]*Vt[:, None]).sum(0)/Vt.sum(); st_fe = MD[nm][0]
+    phi_nr = (1.0/np.clip(grid, 1e-11, None)+inc_src)/np.clip(st_fe, 1e-30, None)
+    Msh = scatter_matrix(mats[nm], GE, source=SRC)
+    d_ = mats[nm].get_nuclide_atom_densities()
+    sab = _macroscopic({k: INC[k] for k in d_}, d_, TS, grid, 101)
+    if sab is None: sab = _macroscopic({k: INC[k] for k in d_}, d_, TS, grid, 102)
+    sscat = np.clip(st_fe-(sab if sab is not None else 0.0), 0, None)
+    Mtr = Msh*(coll(phi_loc, sscat)/np.clip(Msh.sum(1), 1e-30, None))[:, None]
+    mwt = lib("mw", nm, 't'); slt = lib("slab", nm, 't'); mws = lib("mw", nm, 's'); sls = lib("slab", nm, 's')
+    print(f"{nm:9}|  {et(coll(phi_nr,st_fe),mwt):4.2f} {et(coll(phi_loc,st_fe),mwt):4.2f} {et(slt,mwt):4.2f} | {rs(Msh,mws):4.2f} {rs(Mtr,mws):4.2f} {rs(sls,mws):4.2f} | {shp(Mtr,mws):.3f} {shp(sls,mws):.3f}", flush=True)

sphere_sn.py

@@ -0,0 +1,78 @@
+"""1D spherical Sn fixed-source solver. Schemes:
+  'diamond' : weighted-diamond (2nd order; can go negative on resonant problems)
+  'step'    : step characteristic (upwind space+angle); unconditionally positive but diffusive
+  'fixup'   : diamond + negative-flux fixup (set negative outgoing edges to 0 and re-solve the
+              cell). 2nd-order where positive, local upwind only where needed -> accurate AND
+              positive. DEFAULT.
+Curvature coeffs alpha_{m+1/2} >= 0 (tent 0->peak->0), which guarantees the fixup terminates
+positive. Validated vs analytic point-source-in-absorber S*exp(-st r)/(4 pi r^2)."""
+import numpy as np
+
+def _cell(srci, sigtV, amu, Aout, Ain, pin, c, a_hi, a_lo, pang, fixup):
+    """Solve one cell for psi_cell and the outgoing spatial (pout) + angular (pmid) edges.
+    Balance: amu*(Aout*pout - Ain*pin) + c*(a_hi*pmid - a_lo*pang) + sigtV*psc = srci,
+    with diamond closures pout=2psc-pin, pmid=2psc-pang unless an edge is clamped to 0."""
+    co = cm = False
+    for _ in range(3):
+        coef = sigtV + (0.0 if co else 2*amu*Aout) + (0.0 if cm else 2*c*a_hi)
+        const = srci + amu*Ain*pin + c*a_lo*pang + (0.0 if co else amu*Aout*pin) + (0.0 if cm else c*a_hi*pang)
+        psc = const/(coef + 1e-30)
+        pout = 0.0 if co else 2*psc - pin
+        pmid = 0.0 if cm else 2*psc - pang
+        if not fixup or (pout >= 0 and pmid >= 0):
+            break
+        if pout < 0: co = True
+        if pmid < 0: cm = True
+    return psc, pout, pmid
+
+def solve_sphere(redge, sigt, Svol, mu, w, alpha, scheme='fixup'):
+    I = len(sigt); M = len(mu)
+    A = 4*np.pi*redge**2
+    V = 4*np.pi/3*(redge[1:]**3 - redge[:-1]**3)
+    src = 0.5*Svol*V
+    fixup = (scheme == 'fixup'); step = (scheme == 'step')
+    psi = np.zeros((I, M))
+    # starting direction mu=-1 (no angular redistribution), inward
+    pin = 0.0; psi_edge = np.zeros(I)
+    for i in range(I-1, -1, -1):
+        if step:
+            psc = (src[i] + A[i+1]*pin)/(A[i] + sigt[i]*V[i] + 1e-30); pout = psc
+        else:
+            psc, pout, _ = _cell(src[i], sigt[i]*V[i], 1.0, A[i], A[i+1], pin, 0.0, 0.0, 0.0, 0.0, fixup)
+        psi_edge[i] = psc; pin = pout
+    for m in range(M):
+        a_lo = alpha[m]; a_hi = alpha[m+1]; cur = np.zeros(I); amu = abs(mu[m])
+        if mu[m] < 0:                                              # inward: Aout=A[i], Ain=A[i+1]
+            pin = 0.0
+            for i in range(I-1, -1, -1):
+                c = (A[i+1]-A[i])/w[m]
+                if step:
+                    psc = (src[i] + amu*A[i+1]*pin + c*a_lo*psi_edge[i])/(amu*A[i] + c*a_hi + sigt[i]*V[i] + 1e-30); pout = psc; cur[i] = psc; pin = pout; psi_edge[i] = psc
+                else:
+                    psc, pout, pmid = _cell(src[i], sigt[i]*V[i], amu, A[i], A[i+1], pin, c, a_hi, a_lo, psi_edge[i], fixup)
+                    cur[i] = psc; pin = pout; psi_edge[i] = pmid
+        else:                                                     # outward: Aout=A[i+1], Ain=A[i]
+            pin = 0.0
+            for i in range(I):
+                c = (A[i+1]-A[i])/w[m]
+                if step:
+                    psc = (src[i] + amu*A[i]*pin + c*a_lo*psi_edge[i])/(amu*A[i+1] + c*a_hi + sigt[i]*V[i] + 1e-30); pout = psc; cur[i] = psc; pin = pout; psi_edge[i] = psc
+                else:
+                    psc, pout, pmid = _cell(src[i], sigt[i]*V[i], amu, A[i+1], A[i], pin, c, a_hi, a_lo, psi_edge[i], fixup)
+                    cur[i] = psc; pin = pout; psi_edge[i] = pmid
+        psi[:, m] = cur
+    return psi @ w
+
+if __name__ == "__main__":
+    R = 20.0; I = 400; st = 0.1; S = 1.0
+    redge = np.linspace(0, R, I+1); rc = 0.5*(redge[1:]+redge[:-1])
+    sigt = np.full(I, st); Svol = np.zeros(I)
+    V = 4*np.pi/3*(redge[1:]**3-redge[:-1]**3); Svol[0] = S/V[0]
+    mu, w = np.polynomial.legendre.leggauss(16); alpha = np.zeros(17)
+    for m in range(16): alpha[m+1] = alpha[m] - mu[m]*w[m]
+    ana = S*np.exp(-st*rc)/(4*np.pi*rc**2); msk = (rc > 2) & (rc < 16)
+    print("spherical Sn validation (point source in uniform absorber):")
+    for sch in ('diamond', 'step', 'fixup'):
+        phi = solve_sphere(redge, sigt, Svol, mu, w, alpha, scheme=sch)
+        err = np.abs(phi[msk]-ana[msk])/ana[msk]
+        print(f"  {sch:8}: mean {100*err.mean():.2f}%  max {100*err.max():.2f}%  min phi {phi.min():.2e}")