From 71a7fb09eb195c9988c715813679152f568a8ce9 Mon Sep 17 00:00:00 2001 From: shimwell Date: Sat, 27 Jun 2026 21:50:50 +0200 Subject: [PATCH 1/4] Start real-geometry deterministic transport branch Supersedes the 1D-slab prototype (#114): solve the deterministic transport in the ACTUAL reference geometry (concentric spheres, central source) -> 1D spherical Sn, so the near-source shells get the correct (harder) flux instead of the slab's over-softened backscatter. Each material collapsed against its own shell's local flux. No Monte Carlo. Goal: beat slab everywhere incl. the near-source metals. From c9b4eba805130d6eafe39eaf5f5cfd3f9c6fed5c Mon Sep 17 00:00:00 2001 From: shimwell Date: Sat, 27 Jun 2026 22:05:21 +0200 Subject: [PATCH 2/4] Add validated 1D spherical Sn solver + real-geometry transport (negativity open) sphere_sn.py: weighted-diamond curvilinear Sn, central point source, VALIDATED vs analytic point-source-in-absorber (0.9% mean). det_realgeom.py: uses it for the real concentric-sphere geometry. OPEN: diamond negativity on the resonant problem (steel 7.37; naive clamp breaks conservation) -> needs positivity-preserving curvilinear scheme (step characteristic). Deep materials already won by the slab transport (#114); spherical only needed for near-source shells where NR is already excellent. See FINDING_realgeom.md. --- FINDING_realgeom.md | 31 +++++++++++++++ det_realgeom.py | 97 +++++++++++++++++++++++++++++++++++++++++++++ sphere_sn.py | 72 +++++++++++++++++++++++++++++++++ 3 files changed, 200 insertions(+) create mode 100644 FINDING_realgeom.md create mode 100644 det_realgeom.py create mode 100644 sphere_sn.py diff --git a/FINDING_realgeom.md b/FINDING_realgeom.md new file mode 100644 index 00000000000..5936fa87f01 --- /dev/null +++ b/FINDING_realgeom.md @@ -0,0 +1,31 @@ +# 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. diff --git a/det_realgeom.py b/det_realgeom.py new file mode 100644 index 00000000000..84726408d0b --- /dev/null +++ b/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) diff --git a/sphere_sn.py b/sphere_sn.py new file mode 100644 index 00000000000..0508afaa856 --- /dev/null +++ b/sphere_sn.py @@ -0,0 +1,72 @@ +"""1D spherical Sn (weighted-diamond, curvilinear angular redistribution) fixed-source +solver, with validation against the analytic point-source-in-uniform-absorber flux +phi(r) = S*exp(-Sigma_t r)/(4 pi r^2). Used for the real-geometry deterministic transport.""" +import numpy as np + +def solve_sphere(redge, sigt, Svol, mu, w, alpha): + """One-group 1D spherical Sn fixed-source solve. + redge: cell edges r_0..r_I (r_0=0 center). sigt[I], Svol[I] = isotropic source/vol. + mu,w: Gauss-Legendre nodes/weights ([-1,1], sum w=2), ascending. alpha[M+1] curvature edges. + Returns phi[I] = sum_m w_m psi_{i,m}. Vacuum at outer; symmetry at center. + Down-scatter only within a group -> single source iteration (no within-group scatter).""" + I = len(sigt); M = len(mu) + A = 4*np.pi*redge**2 # surface area at each edge (A[0]=0) + V = 4*np.pi/3*(redge[1:]**3 - redge[:-1]**3) + src = 0.5*Svol*V # isotropic angular source per cell + psi = np.zeros((I, M)) # cell-avg angular flux + # --- starting direction mu = -1 (no angular redistribution), sweep inward --- + psi_edge_m = np.zeros(I+1) # angular-edge (m-1/2) cell-avg flux, init mu=-1 + pin = 0.0 # vacuum at outer edge (incoming for inward) + psm = np.zeros(I) + for i in range(I-1, -1, -1): # inward: in=outer edge, out=inner edge + # -(A_{i+1}psi_out_outeredge ...): mu=-1 streaming -1*(A[i+1]*pin - A[i]*pout) + # balance: -1*(A[i+1]*pin - A[i]*pout) + sigt*V*psi = src ; psi=0.5(pin+pout) + # pout = 2 psi - pin + den = A[i] + sigt[i]*V[i] + 1e-30 + psi_s = (src[i] + 0.5*(A[i+1]+A[i])*pin*0 + A[i+1]*0 + pin*(A[i] *0) ) # placeholder + # solve: -(A[i+1]pin - A[i](2psi-pin)) + sigt V psi = src + # = -A[i+1]pin + 2A[i]psi - A[i]pin + sigt V psi = src + # psi(2A[i] + sigt V) = src + (A[i+1]+A[i]) pin + psi_s = (src[i] + (A[i+1]+A[i])*pin) / (2*A[i] + sigt[i]*V[i] + 1e-30) + psm[i] = psi_s; pin = 2*psi_s - pin # pout becomes next inner cell's incoming + psi_edge_prev = psm.copy() # psi at mu-edge 1/2 (the mu=-1 start) + # --- ordinate sweep m=0..M-1 (mu ascending: negatives first=inward, then positives=outward) --- + for m in range(M): + a_lo = alpha[m]; a_hi = alpha[m+1] + cur = np.zeros(I) + if mu[m] < 0: # inward sweep + pin = 0.0 + for i in range(I-1, -1, -1): + # streaming mu(A[i+1]pin - A[i]pout), pout=2psi-pin + c_ang = (A[i+1]-A[i])/w[m] + num = src[i] - mu[m]*(A[i+1]+A[i])*pin + c_ang*(a_hi+a_lo)*psi_edge_prev[i] + den = -2*mu[m]*A[i] + 2*c_ang*a_hi + sigt[i]*V[i] + 1e-30 + ps = num/den; cur[i] = ps; pin = 2*ps - pin + else: # outward sweep + pin = 0.0 # symmetry at center: incoming = outgoing of mu=-mu; approx 0 net at r=0 + for i in range(I): + c_ang = (A[i+1]-A[i])/w[m] + num = src[i] + mu[m]*(A[i+1]+A[i])*pin + c_ang*(a_hi+a_lo)*psi_edge_prev[i] + den = 2*mu[m]*A[i+1] + 2*c_ang*a_hi + sigt[i]*V[i] + 1e-30 + ps = num/den; cur[i] = ps; pin = 2*ps - pin + psi[:, m] = cur + psi_edge_prev = 2*cur - psi_edge_prev # angular-edge recursion psi_{m+1/2}=2psi_m-psi_{m-1/2} + return psi @ w # scalar flux + +if __name__ == "__main__": + # validation: uniform absorber, point source at center -> phi = S exp(-st r)/(4pi r^2) + 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] # point source in cell 0 + M = 16; mu, w = np.polynomial.legendre.leggauss(M) + alpha = np.zeros(M+1) + for m in range(M): alpha[m+1] = alpha[m] - mu[m]*w[m] + phi = solve_sphere(redge, sigt, Svol, mu, w, alpha) + ana = S*np.exp(-st*rc)/(4*np.pi*rc**2) + msk = (rc > 2) & (rc < 16) + err = np.abs(phi[msk]-ana[msk])/ana[msk] + print(f"spherical Sn validation (uniform absorber, point source):") + print(f" mean rel err vs analytic exp(-st r)/(4pi r^2): {100*err.mean():.1f}% max {100*err.max():.1f}%") + for rr in (4, 8, 12): + j = np.argmin(np.abs(rc-rr)); print(f" r={rr}: Sn {phi[j]:.4e} analytic {ana[j]:.4e} ratio {phi[j]/ana[j]:.3f}") From 4762b3b60000ff395a1b043e69b3b2fd0552acf8 Mon Sep 17 00:00:00 2001 From: shimwell Date: Sat, 27 Jun 2026 23:09:53 +0200 Subject: [PATCH 3/4] Add positivity-preserving step-characteristic scheme to spherical Sn alpha_{m+1/2}>=0 (tent) => step (upwind space+angle) is unconditionally positive. scheme='step' (default) in sphere_sn.py; validated (1.2% vs analytic, positive flux). Fixes the diamond negativity on the resonant stack: tungsten 2.13 (beats slab 2.56), Fe-56 1.08 (beats slab 1.39); but step's 1st-order diffusivity hurts steel/CuCrZr. Next: diamond + negative-flux fixup for accuracy + positivity. See FINDING_realgeom.md. --- FINDING_realgeom.md | 31 +++++++++++++ sphere_sn.py | 106 ++++++++++++++++++++------------------------ 2 files changed, 79 insertions(+), 58 deletions(-) diff --git a/FINDING_realgeom.md b/FINDING_realgeom.md index 5936fa87f01..78810cdd9c5 100644 --- a/FINDING_realgeom.md +++ b/FINDING_realgeom.md @@ -29,3 +29,34 @@ treats the angular edges). This is the next engineering step. 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. diff --git a/sphere_sn.py b/sphere_sn.py index 0508afaa856..48c09d3c212 100644 --- a/sphere_sn.py +++ b/sphere_sn.py @@ -1,72 +1,62 @@ -"""1D spherical Sn (weighted-diamond, curvilinear angular redistribution) fixed-source -solver, with validation against the analytic point-source-in-uniform-absorber flux -phi(r) = S*exp(-Sigma_t r)/(4 pi r^2). Used for the real-geometry deterministic transport.""" +"""1D spherical Sn fixed-source solver. Two schemes: + scheme='diamond' : weighted-diamond (2nd order, can go negative on resonant problems) + scheme='step' : step characteristic (upwind in space AND angle). Since the curvature + coefficients alpha_{m+1/2} are >= 0 (tent: 0 -> peak -> 0), step is + UNCONDITIONALLY POSITIVE -- robust for resonant/optically-thick cells. +Validated against the analytic point-source-in-uniform-absorber flux S*exp(-st r)/(4 pi r^2).""" import numpy as np -def solve_sphere(redge, sigt, Svol, mu, w, alpha): - """One-group 1D spherical Sn fixed-source solve. - redge: cell edges r_0..r_I (r_0=0 center). sigt[I], Svol[I] = isotropic source/vol. - mu,w: Gauss-Legendre nodes/weights ([-1,1], sum w=2), ascending. alpha[M+1] curvature edges. - Returns phi[I] = sum_m w_m psi_{i,m}. Vacuum at outer; symmetry at center. - Down-scatter only within a group -> single source iteration (no within-group scatter).""" +def solve_sphere(redge, sigt, Svol, mu, w, alpha, scheme='step'): I = len(sigt); M = len(mu) - A = 4*np.pi*redge**2 # surface area at each edge (A[0]=0) + A = 4*np.pi*redge**2 V = 4*np.pi/3*(redge[1:]**3 - redge[:-1]**3) - src = 0.5*Svol*V # isotropic angular source per cell - psi = np.zeros((I, M)) # cell-avg angular flux - # --- starting direction mu = -1 (no angular redistribution), sweep inward --- - psi_edge_m = np.zeros(I+1) # angular-edge (m-1/2) cell-avg flux, init mu=-1 - pin = 0.0 # vacuum at outer edge (incoming for inward) - psm = np.zeros(I) - for i in range(I-1, -1, -1): # inward: in=outer edge, out=inner edge - # -(A_{i+1}psi_out_outeredge ...): mu=-1 streaming -1*(A[i+1]*pin - A[i]*pout) - # balance: -1*(A[i+1]*pin - A[i]*pout) + sigt*V*psi = src ; psi=0.5(pin+pout) - # pout = 2 psi - pin - den = A[i] + sigt[i]*V[i] + 1e-30 - psi_s = (src[i] + 0.5*(A[i+1]+A[i])*pin*0 + A[i+1]*0 + pin*(A[i] *0) ) # placeholder - # solve: -(A[i+1]pin - A[i](2psi-pin)) + sigt V psi = src - # = -A[i+1]pin + 2A[i]psi - A[i]pin + sigt V psi = src - # psi(2A[i] + sigt V) = src + (A[i+1]+A[i]) pin - psi_s = (src[i] + (A[i+1]+A[i])*pin) / (2*A[i] + sigt[i]*V[i] + 1e-30) - psm[i] = psi_s; pin = 2*psi_s - pin # pout becomes next inner cell's incoming - psi_edge_prev = psm.copy() # psi at mu-edge 1/2 (the mu=-1 start) - # --- ordinate sweep m=0..M-1 (mu ascending: negatives first=inward, then positives=outward) --- + src = 0.5*Svol*V + psi = np.zeros((I, M)) + # starting direction mu=-1 (no angular redistribution), inward sweep + pin = 0.0; psi_edge = np.zeros(I) + for i in range(I-1, -1, -1): + if scheme == 'step': # psi_cell = psi_out (inner edge) + psc = (src[i] + A[i+1]*pin) / (A[i] + sigt[i]*V[i] + 1e-30) + psi_edge[i] = psc; pin = psc + else: + psc = (src[i] + (A[i+1]+A[i])*pin) / (2*A[i] + sigt[i]*V[i] + 1e-30) + psi_edge[i] = psc; pin = 2*psc - pin for m in range(M): - a_lo = alpha[m]; a_hi = alpha[m+1] - cur = np.zeros(I) - if mu[m] < 0: # inward sweep + a_lo = alpha[m]; a_hi = alpha[m+1]; cur = np.zeros(I) + if mu[m] < 0: # inward: incoming = outer edge pin = 0.0 for i in range(I-1, -1, -1): - # streaming mu(A[i+1]pin - A[i]pout), pout=2psi-pin - c_ang = (A[i+1]-A[i])/w[m] - num = src[i] - mu[m]*(A[i+1]+A[i])*pin + c_ang*(a_hi+a_lo)*psi_edge_prev[i] - den = -2*mu[m]*A[i] + 2*c_ang*a_hi + sigt[i]*V[i] + 1e-30 - ps = num/den; cur[i] = ps; pin = 2*ps - pin - else: # outward sweep - pin = 0.0 # symmetry at center: incoming = outgoing of mu=-mu; approx 0 net at r=0 + c = (A[i+1]-A[i])/w[m] + if scheme == 'step': + psc = (src[i] - mu[m]*A[i+1]*pin + c*a_lo*psi_edge[i]) / (-mu[m]*A[i] + c*a_hi + sigt[i]*V[i] + 1e-30) + cur[i] = psc; pin = psc + else: + psc = (src[i] - mu[m]*(A[i+1]+A[i])*pin + c*(a_hi+a_lo)*psi_edge[i]) / (-2*mu[m]*A[i] + 2*c*a_hi + sigt[i]*V[i] + 1e-30) + cur[i] = psc; pin = 2*psc - pin + else: # outward: incoming = inner edge + pin = 0.0 for i in range(I): - c_ang = (A[i+1]-A[i])/w[m] - num = src[i] + mu[m]*(A[i+1]+A[i])*pin + c_ang*(a_hi+a_lo)*psi_edge_prev[i] - den = 2*mu[m]*A[i+1] + 2*c_ang*a_hi + sigt[i]*V[i] + 1e-30 - ps = num/den; cur[i] = ps; pin = 2*ps - pin + c = (A[i+1]-A[i])/w[m] + if scheme == 'step': + psc = (src[i] + mu[m]*A[i]*pin + c*a_lo*psi_edge[i]) / (mu[m]*A[i+1] + c*a_hi + sigt[i]*V[i] + 1e-30) + cur[i] = psc; pin = psc + else: + psc = (src[i] + mu[m]*(A[i+1]+A[i])*pin + c*(a_hi+a_lo)*psi_edge[i]) / (2*mu[m]*A[i+1] + 2*c*a_hi + sigt[i]*V[i] + 1e-30) + cur[i] = psc; pin = 2*psc - pin psi[:, m] = cur - psi_edge_prev = 2*cur - psi_edge_prev # angular-edge recursion psi_{m+1/2}=2psi_m-psi_{m-1/2} - return psi @ w # scalar flux + psi_edge = cur if scheme == 'step' else (2*cur - psi_edge) # angular-edge update + return psi @ w if __name__ == "__main__": - # validation: uniform absorber, point source at center -> phi = S exp(-st r)/(4pi r^2) 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] # point source in cell 0 - M = 16; mu, w = np.polynomial.legendre.leggauss(M) - alpha = np.zeros(M+1) - for m in range(M): alpha[m+1] = alpha[m] - mu[m]*w[m] - phi = solve_sphere(redge, sigt, Svol, mu, w, alpha) - ana = S*np.exp(-st*rc)/(4*np.pi*rc**2) - msk = (rc > 2) & (rc < 16) - err = np.abs(phi[msk]-ana[msk])/ana[msk] - print(f"spherical Sn validation (uniform absorber, point source):") - print(f" mean rel err vs analytic exp(-st r)/(4pi r^2): {100*err.mean():.1f}% max {100*err.max():.1f}%") - for rr in (4, 8, 12): - j = np.argmin(np.abs(rc-rr)); print(f" r={rr}: Sn {phi[j]:.4e} analytic {ana[j]:.4e} ratio {phi[j]/ana[j]:.3f}") + 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'): + 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():.1f}% max {100*err.max():.1f}% min phi {phi.min():.2e}") From 5062c6218216b725833164f1bec55949498a6507 Mon Sep 17 00:00:00 2001 From: shimwell Date: Sat, 27 Jun 2026 23:23:48 +0200 Subject: [PATCH 4/4] Add diamond+negative-flux-fixup scheme to spherical Sn (scheme='fixup') Set-to-zero & re-solve; validated best on smooth (0.74% vs diamond 0.89, step 1.17), positive flux. BUT on the resonant stack the set-to-zero clamp accumulates more error than step (Fe-56 2.06 vs step 1.08) -- many resonance negativities. Conclusion: no spherical scheme beats the robust 1D-slab transport (#114, Fe-56 0.60) on deep resonant materials; spherical only helps near-source. Recommend the #113+#114+NR-near-source hybrid. See FINDING_realgeom.md. --- FINDING_realgeom.md | 33 ++++++++++++++++++++ sphere_sn.py | 74 +++++++++++++++++++++++++++------------------ 2 files changed, 78 insertions(+), 29 deletions(-) diff --git a/FINDING_realgeom.md b/FINDING_realgeom.md index 78810cdd9c5..158932c4f94 100644 --- a/FINDING_realgeom.md +++ b/FINDING_realgeom.md @@ -60,3 +60,36 @@ accuracy AND positivity. 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. diff --git a/sphere_sn.py b/sphere_sn.py index 48c09d3c212..179151e06c0 100644 --- a/sphere_sn.py +++ b/sphere_sn.py @@ -1,50 +1,66 @@ -"""1D spherical Sn fixed-source solver. Two schemes: - scheme='diamond' : weighted-diamond (2nd order, can go negative on resonant problems) - scheme='step' : step characteristic (upwind in space AND angle). Since the curvature - coefficients alpha_{m+1/2} are >= 0 (tent: 0 -> peak -> 0), step is - UNCONDITIONALLY POSITIVE -- robust for resonant/optically-thick cells. -Validated against the analytic point-source-in-uniform-absorber flux S*exp(-st r)/(4 pi r^2).""" +"""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 solve_sphere(redge, sigt, Svol, mu, w, alpha, scheme='step'): +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 sweep + # 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 scheme == 'step': # psi_cell = psi_out (inner edge) - psc = (src[i] + A[i+1]*pin) / (A[i] + sigt[i]*V[i] + 1e-30) - psi_edge[i] = psc; pin = psc + if step: + psc = (src[i] + A[i+1]*pin)/(A[i] + sigt[i]*V[i] + 1e-30); pout = psc else: - psc = (src[i] + (A[i+1]+A[i])*pin) / (2*A[i] + sigt[i]*V[i] + 1e-30) - psi_edge[i] = psc; pin = 2*psc - pin + 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) - if mu[m] < 0: # inward: incoming = outer edge + 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 scheme == 'step': - psc = (src[i] - mu[m]*A[i+1]*pin + c*a_lo*psi_edge[i]) / (-mu[m]*A[i] + c*a_hi + sigt[i]*V[i] + 1e-30) - cur[i] = psc; pin = psc + 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 = (src[i] - mu[m]*(A[i+1]+A[i])*pin + c*(a_hi+a_lo)*psi_edge[i]) / (-2*mu[m]*A[i] + 2*c*a_hi + sigt[i]*V[i] + 1e-30) - cur[i] = psc; pin = 2*psc - pin - else: # outward: incoming = inner edge + 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 scheme == 'step': - psc = (src[i] + mu[m]*A[i]*pin + c*a_lo*psi_edge[i]) / (mu[m]*A[i+1] + c*a_hi + sigt[i]*V[i] + 1e-30) - cur[i] = psc; pin = psc + 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 = (src[i] + mu[m]*(A[i+1]+A[i])*pin + c*(a_hi+a_lo)*psi_edge[i]) / (2*mu[m]*A[i+1] + 2*c*a_hi + sigt[i]*V[i] + 1e-30) - cur[i] = psc; pin = 2*psc - pin + 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 - psi_edge = cur if scheme == 'step' else (2*cur - psi_edge) # angular-edge update return psi @ w if __name__ == "__main__": @@ -56,7 +72,7 @@ def solve_sphere(redge, sigt, Svol, mu, w, alpha, scheme='step'): 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'): + 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():.1f}% max {100*err.max():.1f}% min phi {phi.min():.2e}") + print(f" {sch:8}: mean {100*err.mean():.2f}% max {100*err.max():.2f}% min phi {phi.min():.2e}")