Tearing - BUGFIX! - Convert Δ' to the r_s reference length before slab-layer matching - #403
Tearing - BUGFIX! - Convert Δ' to the r_s reference length before slab-layer matching#403d-burg wants to merge 3 commits into
Conversation
Regression harness report —
|
Shaped-case check: shipped DIII-D-like SLAYER example (auto grid,
|
| q | k_ref | μ | K^(2μ) | Δ' raw | Δ' converted |
|---|---|---|---|---|---|
| 2 | 0.9125 | 0.540 | 0.9059 | 9.0969 | 8.2405 |
| 3 | 1.2245 | 0.529 | 1.2392 | −6.1961 | −7.6782 |
| 4 | 1.3597 | 0.558 | 1.4092 | −16.375 | −23.076 |
| 5 | 1.5051 | 0.831 | 1.9735 | −2475.3 | −4885.0 |
| 6 | 1.5502 | 0.521 | 1.5787 | −194.65 | −307.30 |
| 7 | 1.5541 | 0.504 | 1.5599 | 242.23 | 377.86 |
Notes for review:
- q=2 (the 2/1): −9.4 % in Δ'; with Δ_crit = 14.42 the effective drive moves −5.32 → −6.18. No sign changes on any surface (the q=2 flip threshold is 1.585).
- Shaping puts K > 1 outboard of q≈2.7, so the conversion amplifies |Δ'| at q ≥ 3 — this is why γ moved harder on outer surfaces in the regression report.
- q=5 has μ = 0.83 (strongly non-slab, D_I = −0.69); its ×1.97 factor is consistent bookkeeping, but the slab-layer match itself is a poor approximation at that surface regardless of convention.
🤖 Generated with Claude Code
Benchmark figures: Δ' alignment with TJ, before and after the K^(2μ) conversionFrozen TJ circular benchmark scans (q_edge = 3.6): TJ = hollow black rings (its own r_s-referenced Δ̂), julia_GPEC raw = dashed line (ψ_N-referenced BVP Δ', the pre-PR pairing), julia_GPEC ×K^(2μ) and Fortran STRIDE ×K^(2μ) = filled markers. K and μ are computed per point/surface from TJ's own equilibrium files — nothing fitted. Grey vertical band = pole region (each code approaches its own ideal-kink limit there; TJ's β-limit sits at p_f = 0.1675 vs 0.1851 for both STRIDE codes, so the comparison is not meaningful past the band edge).
Same data as percentage difference from TJ (grey horizontal band = ±3 % combined grid uncertainty of the two STRIDE codes):
Reading guide:
Figures live on the non-merging 🤖 Generated with Claude Code |
Same comparison on fixed grids: julia
|
| series | β roughness (med/max) | ε roughness (med/max) |
|---|---|---|
| julia auto@1e-3 | 1.38 % / 7.94 % | 0.91 % / 1.89 % |
| julia ldp@256 | 0.71 % / 2.45 % | 0.16 % / 0.60 % |
| fortran ldp@128 | 0.34 % / 1.75 % | 0.20 % / 1.57 % |
| fortran ldp@256 | 3.64 % / 6.83 % | 2.12 % / 10.24 % |
(The noise is a property of the auto grid's per-equilibrium adaptivity, not of the solver — julia's own fixed-grid series is the smoothest of the four. Fortran mpsi=256 is excluded from the figures because ldp knot oversampling above mpsi≈256 injects solver noise.)
With the fixed grid, the corrected 2/1 medians are 1.4 % (β) and 1.4 % (ε) from TJ — spreads [−6.1 %, +0.2 %] and [−3.7 %, −1.0 %] — and the two STRIDE implementations are nearly indistinguishable point-for-point after the conversion. Median offsets are grid-robust (16.7–18.5 % raw across all four grid choices), so the K^(2μ) result is unchanged; only the scatter tightens.
Generated with JULIA_SRC=ldp via the same scripts (see PR body links).
🤖 Generated with Claude Code
|
Cross-reference: the "unexplained low-ε 3/1 drift" flagged in this PR's benchmark comments is now root-caused and fixed — it was a Riccati shooting error-control bug (#412, draft PR #413), not toroidal physics. With #413's branch the ε-scan 3/1 median |difference from TJ| after the K^(2μ) conversion drops 5.7 % → 2.3 % (see the scan-scale validation on #413), while the 2/1 panels here are unaffected (Δ'_21 moves ≤ 0.11 % across all 196 scan points). The two PRs are independent; merging #413 first and rebasing this branch is the suggested order. 🤖 Generated with Claude Code |
0cfee9d to
8315b5a
Compare
Rebased onto develop (
|
Radial-label unification (commit 06be788)A parallel investigation into a resistive-layer-overlap truncation criterion (psihigh cap, PR to follow; #298 context) independently ran into the same radial-label question this PR addresses, and the two efforts collided in
Findings on record from the cross-validation of the two investigations
Tests: LayerInputs 56 (incl. a five-label analytic-vs-FD derivative testset), Runner 79/79, params green. Regression harness refresh of 8315b5a -> 06be788 on NO MERGE WITHOUT HUMAN REVIEW — this PR remains a draft and requires a named human reviewer and assignee before it can be considered for merge. This is non-negotiable. |
Regression refresh: 8315b5a -> 06be788 (
|
| Quantity | Diff | Status |
|---|---|---|
| surface indices, m, n, minor radius rs | 0.0e+00 | OK |
| r-based shear | 6.0e-08 (0.00%) | changed |
| Lundquist S | 7.0e-02 abs (~1e-9 rel) | changed |
| D_norm | 3.1e-10 (0.00%) | changed |
| P_perp, tauk, iota_e | 0 / 5.5e-15 / 1.1e-16 | OK |
| Q_root [2/1,3/1,4/1] | 5.8e-06 | OK |
| omega_Hz | 1.9e-01 | OK |
| gamma_Hz | 1.5e-01 (0.01% max element) | changed |
| no_root flags, enabled | 0.0e+00 | OK |
Interpretation: rs is bit-identical (the default :midplane radius formula is untouched); the 1e-8-1e-9 relative movement in shear/S/D_norm is precisely the h^2 ~ 1e-10 truncation error of the removed finite-difference stencil, i.e. the analytic derivative correcting the FD approximation, and it amplifies to at most 0.01 percent in the extracted growth rates through root-finding. The refactor is inert at defaults to within the old stencil's own numerical noise. No tolerance or acceptance criterion was changed anywhere in this commit.
Post-hoc annotation: the omega_Hz (1.9e-01) and gamma_Hz (0.01%) lines in this table were later found to be contaminated by a harness cached-sample artifact (see the forced-rerun verification comment below) and are upper bounds, not measured effects; the true effect of the derivative swap is at or below the ~1e-5 Hz sample-noise ceiling. The interpretation paragraph's conclusion (inert at defaults) stands, conservatively.
d-burg
changed the title
github-actions
Bot
added
bugfix
|
Forced-rerun verification of the regression stamp. A harness caching subtlety was found today: the results DB keeps one cached sample per commit ( To rule it out for this PR, the stamped One earlier comment in this thread needs a precision note: the intermediate |
…b-layer matching The SLAYER dispersion relation paired a ψ_N-referenced BVP Δ' with an r_s-referenced layer side (Δ(Q), S=τ_R/τ_H on r_s, and the critical-Δ, whose Ŵ_d is W/r_s). Apply the Frobenius reference-length transform Δ̂_ij = K_i^(1/2+μ_i)·Δ'_ij·K_j^(μ_j−1/2), K = r_s·(dψ_N/dr)|_s, μ = √(−D_I), at the matching point. Verified parameter-free against the TJ circular benchmarks (median residual +1% over 23 points; absolute 2/1 agreement 17.5%→2.5% (β) and 17.6%→1.1% (ε)). BVP Δ' outputs are unchanged; only the SLAYER matching (γ, Δ_eff) moves. GGJ untouched. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…el options to build_slayer_inputs The cylindrical layer model needs a single minor radius; on a shaped torus that label is ambiguous. Add rs_method options :halfwidth (midplane half-chord, the shift-free stand-in for the circular-theory flux label; reproduces it to ~1% on circular benchmarks) and :volume (cylinder-equivalent √(V/2π²R₀), the Rutherford-literature convention) alongside the existing :midplane default and :fsa. Every label feeds S, the r-based shear, W_d, and k_ref together, so each choice is self-consistent by construction. Programmatic API only — not exposed via TOML; the default and all TOML-driven results are unchanged. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…nd unify the label plumbing Factor the rs_method radial-label options out of build_slayer_inputs into one exported radial_label(equil; rs_method, theta) returning (r, dr/dpsi) closures, and add a fifth label :flux — Fitzpatrick's toroidal-flux surface label (Nucl. Fusion 2025, Eq. 30), r = sqrt(2 psi_t / B0) with the g = F/(B0 R0) correction carried. All five labels now use analytic psi-derivatives from the interpolants themselves; the finite-difference stencils (which silently clamped within 1e-4 of the flux boundaries) are removed. The :fsa branch gains the max(r^2, 0) guard needed on extrapolated surfaces. k_ref and the K^(2mu) Delta-prime conversion consume the same closures, so the label choice drives S, the r-based shear, W_d and the Delta-prime reference length together for every option. Default :midplane behavior is unchanged up to the analytic-derivative refinement. Shared infrastructure for the resistive layer-overlap psihigh cap branch, which evaluates its criterion in :flux. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
06be788 to
baf6746
Compare
Closing summary: per-surface effect on the reference DIII-D-like equilibriumDirect base-vs-head runs of
What to read from it
Other changes in this PR worth noting for review
NO MERGE WITHOUT HUMAN REVIEW — this PR requires approval from its named reviewers (logan-nc, matt-pharr). Non-negotiable. |
1 participant
Release note
k_ref,mu_mercier, anddelta_prime_conversiondatasets documenting the transform, and theDelta_prime_matrixlong_name notes the r_s reference.SLAYER's dispersion relation compared the outer-region Δ' (referenced to unit ψ_N) against layer quantities referenced to the minor radius r_s, overstating or understating growth rates surface-by-surface. Δ' is now converted to the r_s reference length via the Mercier-exponent K^(2μ) transform before slab-layer matching, and the radial label defining r_s (five options, including Fitzpatrick's toroidal-flux label) is selectable programmatically, with the historical midplane default unchanged.
Regression report
regress --cases diiid_slayer_n1,diiid_n1_riccati --refs 8f4be5fe,baf674659 --force(develop tip vs branch head, both sides force-re-run so neither comes from the results cache; same pinned manifest, julia 1.11.6):The moved quantities are exactly the intent: the dispersion roots move because Δ' now enters the layer matching in the r_s reference. The 1e-8–1e-9 relative movement in shear/S/D_norm is the truncation error of the finite-difference stencil this branch replaces with analytic interpolant derivatives, not a physics change.
diiid_n1_riccatiis bit-identical throughout, including the Δ' BVP diagonal and raw matrix — this PR converts Δ' downstream of the outer-region solve and does not perturb it. No tolerance or acceptance criterion was altered.Summary
The SLAYER tearing dispersion relation
Δ' − Δ_crit = S^{1/3}Δ̂(Q)mixed two reference-length conventions: the outer BVP Δ' is referenced to unit Δψ_N (STRIDE convention), while every term on the layer side — the Fitzpatrick slab Δ(Q), the Lundquist numberS = (μ₀r_s²/η)/τ_H, and therfitzp/lar/toroidalcritical-Δ — is referenced to unitx̂ = (r−r_s)/r_s. This PR converts Δ' to the r_s reference at the matching point:whose diagonal is
K^(2μ)·Δ'_kk, reducing to the textbookΔ̂ = r_s·Δ'_physat the slab point D_I = −¼. K and μ are computed per surface from the equilibrium (surface_minor_radius/surface_da_dpsi, GGJ E+F+H−¼ fromResistGeometry).Why this is correct
A_L|x|^(1/2−μ) + A_S|x|^(1/2+μ)(Glasser, Wang & Park 2016, Eq. 26); rescalingx_ψ = K·x̂mapsÂ_L = A_L K^(1/2−μ),Â_S = A_S K^(1/2+μ), giving the transform above.Ŵ_d = W_d/r_s(proven by its source multiplying the reported width by r_s to convert to a-units), its Δ carriesr_s^(1+2ν)by definition (its docs' pressure-flattening derivation,1+2ν = 2μ), and its dispersion assembly(Δ − Δ_c)/S^{1/3}is identical in form to ours. Thedc_tmpformulas inLayerParameters.jlare a line-for-line port of that r_s-referenced formula set.K^(2μ)computed from TJ's own equilibrium files reproduces the measured Δ̂_TJ/Δ'_GPEC ratio to +1 % median over 23 pole-free points; in absolute values the 2/1 median |difference from TJ| drops 17.5 % → 2.5 % (β) and 17.6 % → 1.1 % (ε), inside the combined grid band. The factor is 0.84 at the 2/1 but 1.20 at the 3/1, and raw GPEC was correspondingly +18 % high / −28 % low — one formula corrects both directions.What changes numerically
SingularSurfaces/Delta_prime_matrixstays ψ_N-referenced; STRIDE benchmarking is unaffected.Tearing/PerSurface/Delta_prime_matrixnow stores the r_s-referenced matrix actually used in the matching (long_name updated), plus new per-surface datasetsk_ref,mu_mercier,delta_prime_conversion.rescale_deltahandles inner→outer units natively).SLAYERParametersdefault tok_ref = 1→ identity conversion, so existing unit tests andrun_slayer_from_inputsworkflows are unaffected unless the caller opts in.Changes
SLAYERParameters: new fieldsk_ref,mu_mercier(documented; default 1.0/0.5 = identity).build_slayer_inputs: fills both from the equilibrium; Mercier-unstable surfaces (D_I ≥ 0) get μ → 0 continuously (Δ' left raw) — documented in-line.run_slayer_from_inputs: appliesdelta_prime_to_rs_reference(new, documented) for the slab-layer path.surface_couplingdocstring corrected:S^{1/3}maps the inner Δ to the r_s-referenced outer convention, not "ψ-units" — the previous wording was the error that concealed this mismatch.Tearing/PerSurfacedatasets + annotations;SLAYERResult.dp_matrixdocstring notes the convention.K^(2μ), off-diagonal split factors, slab-point reduction toK, identity default).Validation
developvs this branch: only the SLAYER dispersion roots moved (γ up to 42 % max-element; ω 0.12 %); all layer inputs and all 17 Riccati/BVP quantities (incl.SingularSurfaces/Delta_prime_matrix) bit-identical — full reportBenchmark figures and the full theory audit (TJ source citations, ratio test, corrected Δ' comparison plots) are recorded in the maintainer's private benchmark repo: issue · figures + scripts (maintainer access).
Important
This changes every SLAYER γ (≈9 % shift in Δ_eff at q=2 on the shipped example). Do not merge without human review of the derivation and the regression report.
🤖 Generated with Claude Code