From 5edaaa7b9327db0e2e851f2bf88d090963641117 Mon Sep 17 00:00:00 2001 From: Jonathan Brodrick Date: Thu, 23 Jul 2026 16:02:04 +0100 Subject: [PATCH] fix(vfp1d): SNB heat-flux diagnostic on a global absolute-energy grid snb_heat_flux discretized the energy groups on a per-cell beta = v^2/(2T) grid, so a group's absolute energy (and hence its mean free path lambda ~ eps^2) varied from cell to cell. Because the SNB multigroup diffusion transports each group at fixed kinetic energy (energy is conserved as electrons stream, which is what lets the groups decouple in energy and couple only in space), the per-cell beta grid made the spatial stencil couple different absolute energies between neighbouring cells -- a spurious energy advection -- and dropped the suprathermal carriers that stream from hot to cold and set the preheat. The peak flux was only mildly affected but the cold-side preheat was under-predicted by ~25% (mild gradients) to ~5x (steep gradients). Rewrite to a global absolute-energy grid eps = v^2/2 (uniform in space, spanning [0, beta_max * max(T)]) with the local Spitzer-Harm group weight W_g(x) = int beta^4 e^-beta/24 dbeta over beta = eps/T(x). This matches the validated reference SNB implementation. Defaults ngroups 32 -> 300 and beta_max 30 -> 20 (the global-energy grid needs a few hundred groups since cold cells sample beta coarsely; beta_max ~ 20 keeps the beta ~ 16 flux carriers on-grid at the hottest cell). Verified against the reference on the Z=7 Gaussian problem: q_snb/q_sh reproduces the reference to ~1e-5 in local, mild, and steep regimes. (The electric-field / l_s stopping term in the effective diffusion mfp is a separate, second-order refinement and is intentionally still omitted.) Co-Authored-By: Claude Opus 4.8 (1M context) --- adept/vfp1d/heat_flux.py | 51 +++++++++++++++++--------- configs/vfp-1d/spherical-heatflow.yaml | 4 +- docs/source/solvers/vfp1d/config.md | 10 +++-- 3 files changed, 41 insertions(+), 24 deletions(-) diff --git a/adept/vfp1d/heat_flux.py b/adept/vfp1d/heat_flux.py index c1e8a548..8c377575 100644 --- a/adept/vfp1d/heat_flux.py +++ b/adept/vfp1d/heat_flux.py @@ -27,7 +27,7 @@ from adept.vfp1d.storage import calc_EH -DEFAULT_SNB_CFG = {"ngroups": 32, "beta_max": 30.0, "r": 2.0} +DEFAULT_SNB_CFG = {"ngroups": 300, "beta_max": 20.0, "r": 2.0} def _geometry_factors(x: np.ndarray, dx: float, geometry: str) -> tuple[np.ndarray, np.ndarray]: @@ -103,8 +103,9 @@ def snb_heat_flux( """Multigroup SNB heat flux at cell edges (nx+1,). Implements the "separated" SNB variant recommended by Brodrick et al. (2017): - for each energy group g (beta = v^2 / 2T, group weights from the Spitzer-Härm - spectrum W(beta) = beta^4 exp(-beta) / 24), + for each energy group g on a *global absolute-energy grid* eps = v^2/2 (with + group weights from the Spitzer-Härm spectrum W(beta) = beta^4 exp(-beta) / 24, + beta = eps/T the local position in the Maxwellian spectrum), H_g / (lambda_g^ee / r) - div( (xi * lambda_g^ei / 3) grad H_g ) = -div(W_g q_SH) @@ -112,6 +113,15 @@ def snb_heat_flux( q_SNB = q_SH - sum_g (xi * lambda_g^ei / 3) grad H_g. + The energy groups are bands of fixed *kinetic energy* eps = v^2/2 (uniform in + space), NOT of the local beta = eps/T. This is required for correctness: the + multigroup diffusion transports each group at fixed energy (energy is conserved + as electrons stream, so the groups decouple in energy and couple only in space). + Gridding on the local beta would make a group's absolute energy - and hence its + mfp lambda ~ eps^2 - vary from cell to cell, so the spatial stencil would couple + different energies between neighbouring cells (a spurious energy advection) and + would drop the suprathermal carriers that stream hot -> cold and set the preheat. + The mean free paths use the same normalized collision rates as the kinetic solver: nu_ei(v) = nuee_coeff * logLam_ratio * Z^2 * ni / v^3 (the l=1 rate in FLMCollisions) and nu_ee(v) = nuee_coeff * n / v^3, so lambda = v^4 / (...). @@ -133,34 +143,39 @@ def snb_heat_flux( area, vol = _geometry_factors(x, dx, geometry) q_sh = spitzer_harm_heat_flux(x, dx, n, T, Z, cfg) - # group edges in beta = v^2/(2T); logarithmic spacing resolves the low-energy - # groups while still covering the suprathermal tail that carries q_SH - beta_edges = np.geomspace(beta_max / 400.0, beta_max, ngroups + 1) - beta_edges[0] = 0.0 + # Global absolute-energy grid eps = v^2/2 (code units), uniform in space and + # spanning eps in [0, beta_max * max(T)]. The flux-carrying groups sit at + # beta = eps/T ~ 16 (v ~ 4 v_th, the v^9 heat-flux moment), so beta_max ~ 20 + # keeps them on the grid even at the hottest cell. (Cold cells then reach + # beta = beta_max * max(T)/T_local >> beta_max, retaining the suprathermal + # carriers that stream in from the hot region and set the preheat.) + eps_edges = np.linspace(0.0, beta_max * float(np.max(T)), ngroups + 1) + eps_c = 0.5 * (eps_edges[1:] + eps_edges[:-1]) # group midpoints (nx-independent) def _cdf(b): # int_0^b beta^4 exp(-beta)/24 dbeta = 1 - exp(-b)(1 + b + b^2/2 + b^3/6 + b^4/24) return 1.0 - np.exp(-b) * (1.0 + b + b**2 / 2.0 + b**3 / 6.0 + b**4 / 24.0) - weights = np.diff(_cdf(beta_edges)) - beta_c = np.sqrt(beta_edges[:-1] * beta_edges[1:]) - beta_c[0] = 0.5 * beta_edges[1] - - xi = (Z + 0.24) / (Z + 4.2) + xi_e = _interp_c2e((Z + 0.24) / (Z + 4.2)) T_edge = _interp_c2e(T) + Z2ni_e = _interp_c2e(Z**2 * ni) q_nonlocal = np.zeros_like(q_sh) H_all = np.zeros((ngroups, x.size)) for g in range(ngroups): - v2_c = 2.0 * beta_c[g] * T # v_g^2 at centers - v2_e = 2.0 * beta_c[g] * T_edge - lam_ee_c = v2_c**2 / (nuee_coeff * n) - lam_ei_e = v2_e**2 / (nuee_coeff * logLam_ratio * _interp_c2e(Z**2 * ni)) - D_e = _interp_c2e(xi) * lam_ei_e / 3.0 + v2 = 2.0 * eps_c[g] # v_g^2, fixed in space (energy is conserved during transport) + lam_ee_c = v2**2 / (nuee_coeff * n) # centers (nx,) + lam_ei_e = v2**2 / (nuee_coeff * logLam_ratio * Z2ni_e) # edges (nx+1,) + D_e = xi_e * lam_ei_e / 3.0 + + # local Spitzer-Härm group weight W_g(x) = int_{eps_g}^{eps_g+1} beta^4 e^-beta/24 dbeta, + # beta = eps / T(x): spatially varying because a fixed absolute-energy band samples a + # different part of the local Maxwellian spectrum as T changes. + w_e = _cdf(eps_edges[g + 1] / T_edge) - _cdf(eps_edges[g] / T_edge) # source: -div(W_g q_SH), finite volume with zero flux through boundary edges - q_g = weights[g] * q_sh + q_g = w_e * q_sh src = -np.diff(area * q_g) / vol # tridiagonal finite-volume diffusion operator diff --git a/configs/vfp-1d/spherical-heatflow.yaml b/configs/vfp-1d/spherical-heatflow.yaml index 1ff61b0d..8c3f3d07 100644 --- a/configs/vfp-1d/spherical-heatflow.yaml +++ b/configs/vfp-1d/spherical-heatflow.yaml @@ -98,6 +98,6 @@ terms: diagnostics: snb: - ngroups: 32 - beta_max: 30.0 + ngroups: 300 + beta_max: 20.0 r: 2.0 diff --git a/docs/source/solvers/vfp1d/config.md b/docs/source/solvers/vfp1d/config.md index a2651d4c..35b68cb6 100644 --- a/docs/source/solvers/vfp1d/config.md +++ b/docs/source/solvers/vfp1d/config.md @@ -370,15 +370,17 @@ Optional post-processing diagnostics configuration. The Spitzer-Härm / SNB / ki | Field | Type | Default | Description | |-------|------|---------|-------------| -| `ngroups` | int | `32` | Number of energy groups in $\beta = v^2/(2T)$ | -| `beta_max` | float | `30.0` | Upper edge of the group grid | +| `ngroups` | int | `300` | Number of energy groups on the absolute-energy grid $\varepsilon = v^2/2$ spanning $[0,\ \beta_\text{max}\max(T)]$ | +| `beta_max` | float | `20.0` | $\beta = \varepsilon/T$ coverage at the hottest cell (sets the grid's top energy $\varepsilon_\text{max} = \beta_\text{max}\max(T)$) | | `r` | float | `2.0` | e-e absorption scaling of the group mfp, $\lambda_g^{ee}/r$ | +The energy groups are bands of fixed *kinetic energy* $\varepsilon = v^2/2$ (uniform in space), so the multigroup diffusion transports each group at fixed energy — required for the nonlocal preheat to be captured correctly. `ngroups` needs to be a few hundred because cold cells sample the grid coarsely in $\beta = \varepsilon/T$; the flux-carrying groups sit at $\beta \approx 16$, so `beta_max` $\gtrsim 20$ keeps them on the grid. + ```yaml diagnostics: snb: - ngroups: 32 - beta_max: 30.0 + ngroups: 300 + beta_max: 20.0 r: 2.0 ```