ForceFreeStates - FEATURE! - Cap the integration domain where resistive layers overlap - #298
ForceFreeStates - FEATURE! - Cap the integration domain where resistive layers overlap#298amlakbekaluachule wants to merge 6 commits into
Conversation
amlakbekaluachule
assigned
amlakbekaluachule and
d-burg and unassigned
amlakbekaluachule
amlakbekaluachule
added
refactor
|
@logan-nc this PR does two things:
Let me know if anything looks off. |
|
I synced this branch with Worth knowing why the retarget was needed: The merge hit two conflicts, both caused by Two things worth a look before review: 1.
|
|
@d-burg whats the status of this? Should it be reviewed now? Or should it still be marked as draft? |
|
Fixed both flagged issues — SetResistiveWidthPsihigh.jl is now included/exported, and δ_FKR/δ_visco are documented in the LayerWidths docstring. |
|
@claude review this |
ReviewReminder (per repo policy): this PR must not be merged without third-party human review, regardless of any automated approval. Also — please run the regression harness ( Correctness / design concerns
Physics traceability
Style / consistency
Documentation / tests
None of the above touch the pre-existing |
StatusRewritten from the original 172-line contribution, which had zero call sites and compared layer widths in metres against positions in normalized flux. That file is gone. Feature — Default is off because enabling it crosses a rational surface on a diverted plasma: on Nine defects fixed along the way, all found by trying to make the feature trustworthy:
Verification — 122 + 67 + 59 + 32 + 14 tests; docs clean; the metres→flux conversion triple-checked three independent ways (the binding pair is 1.821 mm apart with half-widths summing to 2.338 mm, argued in pure metres); the back-off gate validated against ForceFreeStates energies as ground truth, and it flips one step before they degrade. Still open
|
|
@d-burg — ready for review, but please take #399 first. #399 is the one-line That ordering also separates the regression movement. Worth your attention specifically:
Everything except the opt-in flag is a bug fix that stands on its own — the field-line ODE |
logan-nc
changed the title
|
@d-burg, as claude says this has hit a point where it needs your expert review. Workflows are in place to decide the truncation from the resonant layer widths. The issue is that it is reporting big widths and wanting to truncate way in (~0.9892 for DIII-D). I find this suspect, since Fitzpatrick had been finding things like to be very far out past 0.9999. In the process of adding the truncation logic, claude did convince me the units of the layer widths had been wrong in the original overlap-checking-loop (comparing layer widths of m to rational spacing in psin) and a conversion from m to psin was added... so maybe double check that? Do you think we need to try and reproduce the exact case from Fitzpatricks paper on this as a sanity check? |
|
the 1. Blocking: the trial/back-off loop cannot compile on current
|
| (Q̂*, P̂) | Re GPEC | fig 6.2 | Im GPEC | fig 6.3 |
|---|---|---|---|---|
| (0.25, 0.25) | 0.818 | ~0.8–0.9 | −0.222 | ~−0.15…−0.3 |
| (0.5, 0.5) | 0.943 | ~0.9–1.0 | −0.354 | ~−0.3…−0.45 |
| (2.0, 0.5) | 1.351 | ~1.3–1.5 | −1.417 | ~−1.2…−1.35 |
| (4.0, 0.5) | 2.008 | ~1.9–2.1 | −2.225 | beyond −1.8 |
| (0.5, 4.0) | 1.537 | ~1.5 | −0.119 | ~−0.15 |
| (4.0, 4.0) | 1.350 | ~1.3–1.4 | −1.466 | ~−1.35…−1.5 |
(τ = 1 assumed — the caption does not state it; the agreement is the evidence.)
The widths on the actual DIII-D surfaces are also not anomalous. Table 6.2 of the book quotes
|δ_s| = 4.69 mm (low-field reactor) and 2.95 mm (high-field); ours span 1.8–4.4 mm:
| m | ψ | r_s [m] | |δ_s| [m] | δ_visco [m] | S | δ_visco/δ_s |
|---|---|---|---|---|---|---|
| 2 | 0.518 | 0.362 | 3.91e−3 | 2.81e−3 | 1.9e7 | 0.720 |
| 3 | 0.770 | 0.458 | 4.13e−3 | 2.85e−3 | 3.6e7 | 0.689 |
| 4 | 0.893 | 0.504 | 3.85e−3 | 2.62e−3 | 5.5e7 | 0.680 |
| 5 | 0.968 | 0.531 | 4.42e−3 | 2.74e−3 | 3.3e7 | 0.620 |
| 6 | 0.993 | 0.540 | 2.85e−3 | 1.76e−3 | 5.4e7 | 0.615 |
| 7 | 0.999 | 0.542 | 1.82e−3 | 1.20e−3 | 1.4e8 | 0.661 |
So where does 0.9892 come from? Between m=6 (ψ=0.99328) and m=7 (ψ=0.99851) the spacing is
Δψ = 0.0052, while dr/dψ ≈ 0.35 m there, so the two layers occupy δψ ≈ 0.0081 and 0.0052. The
half-widths sum to 0.0067 > 0.0052 — overlap — and cutting at ψ(6) − w(6)/2 gives 0.9893, the
number you are seeing. Nothing is wrong with the width; the surfaces are simply 5 mΨ apart out there.
Two notes on delta_visco. It is a consistent 0.62–0.72 × δ_s across the whole domain, so it is
the narrower width, overlaps later, yields a larger psihigh, and therefore minimum(cands) never
selects it — the viscous criterion is effectively dead code on this deck. And since it is algebraic
(δ_norm·P_perp^(1/6)) rather than from the Riccati solve, it is the weaker of the two anyway. Worth
deciding whether it earns its place. (delta_FKR is a non-issue: it was renamed delta_norm and is
not used as a criterion — the primary width is already the validated Riccati δ_s.)
2b. Why the criterion is nonetheless close to self-fulfilling
This is where I think your suspicion about ~0.9892 is right, and I do not think the metres→ψ_N
conversion is the culprit.
Past the last real surface the edge law is q ~ -A·ln(1-ψ), so surfaces sit at
ψ_m = 1 - exp(-m/(nA)) and their spacing shrinks geometrically as ψ → 1. Meanwhile the layer
width does not shrink: toward the edge T_e falls, η rises, S drops, and δ_s ∝ r_s S^(-1/3)
grows. A shrinking gap against a growing width means overlap is guaranteed at some radius for any
equilibrium — it is a property of extrapolating toward a separatrix, not a measurement of this
plasma.
So "where do layers first overlap" is set by where the geometric spacing crosses the width, and both
inputs there are extrapolations: A from edge_q_law fitted on the outer knots, and the kinetic
profiles continued past the grid. That is consistent with what you are seeing — the criterion always
wants to cut in aggressively, and the psiedge floor ends up binding on every shipped deck.
Fitzpatrick's ψ > 0.9999 is not really the same quantity: TJ places its outermost surface of
interest, it does not ask where layers first touch under an extrapolated q. I would expect these to
disagree by construction.
3. The viscous criterion fails open
da_dpsi (LayerOverlap.jl:257) and delta_s_m (:277) are both finite-checked, each recording a
note and continue-ing. delta_visco is not — :290 pushes lw.delta_visco / abs(da_dpsi)
unguarded. A NaN there makes every comparison in _first_overlap_limit false, so the viscous scan
silently returns "no overlap" with no note, and recommended = minimum(cands) quietly loses its
more conservative half. The bias is outward, which is the unsafe direction, and it is the one path in
this file that does not follow its own record-and-skip discipline.
4. The scan's cost is bounded only by m_max = 2000 layer solves
Every surface in found runs build_slayer_inputs + slayer_layer_thickness, and the latter is a
Riccati solve. The extrapolation terminates on psi_cap, but how many surfaces that takes depends on
A; a slowly-rising edge law reaches the cap late. Worst case is 2000 Riccati solves — and since the
scan runs whenever kinetic profiles are readable, regardless of the flag, that cost lands on runs
that never asked for the feature. Bounding the extrapolated-surface count (or the work) rather than
only the m index would make it predictable. Worth quoting the added runtime on a default deck.
On reproducing the Fitzpatrick case
Done — see §2. It was cheaper than expected because the solver already implements the book's
equations; no new physics was needed, only driving it at prescribed normalized parameters. I would
not chase the width formula any further: it matches the published figures and the published
magnitudes.
I have also pushed two fixes to this branch (0c2b3ac): delta_visco now gets the same
finite-check-and-note treatment delta_s_m and da_dpsi already had (it was the one path that could
fail open), and the layer-solve count is bounded with a note when the bound bites. 122/122 tests pass.
One coverage note on your own observation that the floor governs every DIII-D deck: as it stands the
physics path is never exercised end-to-end by a shipped deck or by runtests_layer_overlap.jl, which
tests the pieces (conversion, q = m/n, the cut rule, the limited-plasma refusal) rather than a run
where the criterion actually sets the domain. A deck that reaches the physics path — or a test that
drives _choose_psihigh with a synthetic scan — would close that.
Also: this is dirty against develop (4 conflicts). In test/runtests_slayer_runner.jl the
cherry-picked 4304caad collides with #399 as merged — #399's source-text test was replaced with a
behavioural one in runtests_slayer_inputs.jl, so take develop's side and drop the cherry-pick.
|
Clean out any approximate visco layer widths from the overlap criteria - those and any other limit approximations can be output in the The Fitzpatrick model this is based off of is https://iopscience.iop.org/article/10.1088/1741-4326/ae4fdd/meta (also at https://arxiv.org/pdf/2511.07666 if you hit a paywall), NOT the book (which is a fine benchmark for the layer calculations but doesn't have the truncation logic in it). Figure 10 shows an example for typical JET parameters that finds 0.998 and Fig. 11 for n=4 finds 0.9938. @d-burg can you have your AI add this paper to the resources then cite it in annotations and the docs where appropriate? Might as well also have it double check it's own logic above based on the text... This feature is not much help if it regularly suggests truncating deep into the pedestal 😛 Maybe it gets better (further out) for an ITER case? |
|
@logan-nc ahh sounds good, I shall dig into this |
| magnetic perturbation" | ||
| - Location: `docs/resources/2025-Fitzpatrick-Response_of_a_magnetically_diverted_tokamak_plasma_to_a_resonant_magnetic_perturbation.pdf` | ||
| - Published: Nucl. Fusion (2025), doi 10.1088/1741-4326/ae4fdd; arXiv:2511.07666 | ||
| - Describes: **the source of the edge-truncation criterion**. q diverges logarithmically at the |
There was a problem hiding this comment.
@d-burg this is way longer than the other one-liner descriptions in this md. One liner goes here. The actual docs webpage should describe the layer widths and truncation logics and cite this paper using cross references.
| separatrix, so rational surfaces accumulate without bound; the paper shows their spacing | ||
| shrinks *faster* than the resistive layer width, so the layers merge in a narrow region | ||
| straddling the separatrix. Only surfaces with `0 < Ψ < 1 - ε_c` need be retained, the rest | ||
| treated as vacuum — which is exactly the empirical truncation existing response codes |
There was a problem hiding this comment.
"which is exactly" is sooooo AI 🙄
also, no other existing response code actually uses this that I know of
| semi-collisional (Eq. 12), and resistive-inertial (Eq. 13). The `P^(1/6)` factor of Eq. (11) | ||
| is the `delta_visco` scale in `LayerWidths`. | ||
|
|
||
| - **Fitzpatrick (2025)**: "Response of a magnetically diverted tokamak plasma to a resonant |
There was a problem hiding this comment.
In-line annotations of the layer widths and/or truncation logic in src should be citing this as per the claude.md instructions
|
@logan-nc read the paper — agreed it changes the picture, and with its width the truncation lands Sect. 5.6 shows the strong shear near the separatrix forces every resonant layer into the One trap for the implementer: the paper's S = τ_R/τ_A uses the shear-free τ_A of its Eq. (74), while Both truncation points on the DIII-D SLAYER deck, n=1, in-grid plus edge-q-law surfaces:
The paper's references: ε_c = 1.5e−3 (JET n=1, fig 9) and 4.8e−3 (n=4, fig 10); our 2.9e−3 sits in On where the cut should land: the paper's prescription (p. 22) is to solve the response equations Was there a reason the feature targeted Your references.md review comments are addressed on the branch (4201f67): the entry is trimmed to |
|
@logan-nc prototyped the psilim version. Three things changed since my last comment, one of The formula I quoted was right, my implementation of it wasn't. I reintroduced the Second, the radial label matters and we had it wrong. Eq. (100) is not covariant — the width
Worth saying plainly that the flux label is not a circular-plasma construct — it is defined Third, and the reason I would now trust the criterion: it reproduces the paper. Running the
n = 1 lands on it; n = 4 is 34% high in eps_c, which I would put on the digitised pedestal. The So, to your original worry: the 0.98 truncations were not real. They came from the ODE width Architecture, per your sign-off: the overlap point is an upper bound on psilim, nothing Still open, and I would not gate this PR on it: the two width channels now disagree by about two The radial-label plumbing landed separately as Regression harness: not run yet on this branch — I am re-running it. I will post the table separately. Expected |
d-burg
changed the title
github-actions
Bot
added
the
changed-results
06be788 to
baf6746
Compare
…e layer width
Adds three algebraic scales to `LayerWidths` alongside the Riccati `delta_s`:
`delta_norm` (r_s S^-1/3), `delta_visco` (the P^1/6 viscous-resistive broadening) and
`delta_dr`, the diffusive-resistive width of Fitzpatrick, Nucl. Fusion 2025
(doi 10.1088/1741-4326/ae4fdd) Eq. (100), which sets the edge layer-overlap criterion.
The mapping of Eq. (100) into SLAYER variables is not literal. The paper's tau_A
(its Eq. 74) carries no shear, whereas SLAYER's tau_h divides by n*s, so
lu = tau_R/tau_h = (n|s|) * (tau_R/tau_A). Substituting that cancels Eq. (100)'s
explicit (n|s|)^(-1/2) exactly, leaving
delta_dr = r_s^(3/2) * P_perp^(1/4) / ( lu^(1/2) * d_beta^(1/2) )
with no shear factor. Carrying the (n|s|)^(-1/2) on top of SLAYER's lu double-counts
the shear and understates the width by sqrt(n|s|) -- a factor of 6.6 at q=7 on the
DIII-D-like deck, enough to report no overlap anywhere.
Validated against the paper's own Sect. 5.8 JET model, which this reproduces at
Psi = 0.9985 for n = 1 (paper: 0.9985) and 0.9936 for n = 4 (paper: 0.9952).
Marked `!` because `LayerWidths` gains fields.
…er-overlap scan Factors the diverging edge model q ~ -A*ln(1-psi) out of the grid-refinement density floor into `edge_q_law`, so one fit serves both that floor and the resistive-layer overlap scan's search for rational surfaces beyond the equilibrium grid. A cubic in psi (on q or on iota) saturates instead of diverging and is not usable for the extrapolation; measured on the DIII-D-like deck, extrapolating from a grid ending at psi=0.97 to the true q=8 surface at psi=0.99976, the log law lands within 5% while a cubic iota undershoots by 18%. `edge_q_law` returns nothing when the diverging model does not describe the equilibrium, judged by knot count, sign of A, and the log fit beating a plain linear one. That test now gates the edge density floor, which was previously applied unconditionally: a limited plasma with finite edge q is no longer packed as if q blew up. Marked `!` because that changes the auto grid, and so results, for limited decks (Solovev, LAR); diverted decks are unaffected. Also adds an optional `psihigh` to `refined_psi_grid`, which truncates the measured knot density so a pass-1 equilibrium can supply the density for a reduced domain.
Locates the q = m/n rational surfaces, computes each one's resistive layer width, and reports the outermost psi at which adjacent layers are still separated -- the outermost domain the matched-asymptotic treatment supports, per Fitzpatrick, Nucl. Fusion 2025 Sect. 5.9. Surfaces beyond the equilibrium grid are located on the shared separatrix edge q-law, so the scan can report an overlap point lying outside psihigh; that is recorded rather than clipped, since it is the bookkeeping a user needs to see whether the domain was constrained by layer physics at all. Widths are compared in normalized flux, converted from metres by the Jacobian of whichever radial label `rs_method` selects -- the same label that sets the shear, tau_R and d_beta/r, because Eq. (100) is not covariant (delta ~ J^(1/2) under a change of radial variable) and is only valid in the coordinate it was derived in. On the DIII-D-like deck the label is worth one rational surface: the cut sits at psi=0.99752 under `:midplane` and 0.99939 under `:flux`, Fitzpatrick's Eq. (30) label. Reports both the |delta_s| and Eq. (100) channels separately so the two can be compared; they disagree by about two surfaces, which is the subject of issue #411.
…ve-layer overlap `sing_lim!` gains an optional `psilim_cap`. Where the resistive layers of adjacent rational surfaces overlap, no surface retains a well-separated inner region and matched asymptotics is not defined, so that location is an upper bound on the useful domain. It is applied as a cap on `qlim` rather than by assigning `psilim` directly, which keeps it a bound rather than a setpoint: the existing `dmlim` / `qhigh` truncation still selects the final surface from inside it, and a cap lying beyond `psihigh` is inert by construction with no special case -- it simply loses the `min`. The domain can therefore only shrink, never grow, so no equilibrium ever has to be re-formed. Measured on the shipped DIII-D decks with the Eq. (100) width in the `:flux` label, the cap changes no retained surface on any of them: it is beyond `psihigh` for the ideal (0.995), riccati (0.995) and gal_resistive (0.993) decks, and on the SLAYER deck (0.9995) it binds by 1.1e-4 in psi without crossing a rational. The scan needs kinetic profiles, which two of those decks do not carry, so it is inert there regardless.
…d Fitzpatrick JET case Rebuilds the model JET equilibrium of Fitzpatrick, Nucl. Fusion 2025 Sect. 5.8 from the paper's own definitions -- the Eqs. (32)-(35) safety factor with alpha-/alpha+ solved self-consistently against Eq. (36), and mtanh edge profiles anchored to its Fig. 8 -- and checks where adjacent resistive layers first overlap against the values the paper reports. The widths are driven through `slayer_parameters` and `slayer_layer_thickness` rather than recomputed here, so the test covers the shipped tau_R / tau_A / d_beta chain and the `delta_dr` formula itself. That distinction matters: the bug this test would have caught -- carrying Eq. (100)'s explicit (n|s|)^(-1/2) on top of a `lu` that already contains it -- lives in exactly that chain, and a test that reimplemented the width would have passed while the shipped one was wrong. Measured against the paper: Psi = 0.9985111 for n = 1 (paper 0.9985, deviation 1.1e-5) and 0.9949738 for n = 4 (paper 0.9952, deviation 2.3e-4). Bounds are set at 1e-4 and 1e-3, leaving roughly 9x and 4x margin; the model is quadrature and root-finding with no BLAS, so platform spread is orders of magnitude below either. The test logs both deviations so a future failure reports how far it moved rather than only that it moved. This is the suite's second check pinned to an external published reference rather than to GPEC's own history.
…'s own reproducibility The `diiid_slayer_n1` γ_Hz threshold of 1e-1 sits below what the threaded dispersion root search reproduces: re-running one commit against itself with fresh samples moved γ by 0.145 Hz at surface 3/1, 0.122 Hz at 2/1 and 0.076 Hz at 4/1, in discrete steps, with every layer input bit-identical. At 1e-1 the case therefore flags a null change as CHANGED roughly half the time. 2.5e-1 clears the largest measured step by 1.7x. This is a stopgap: the underlying fix is a reproducible root search, after which the threshold can come back down.
|
@d-burg please hold for me to find time to skim the full diff. Quick thoughts from the comments:
Great. I like this simplification a lot. Open physics question: Do we want to always cut inside this stochastic boundary (cut off some valid region) or always get at least to this boundary (include a little invalid region in the name of making sure we capture all the valid region)? The former is numerically easier (we don't have to go out as far), but perhaps the later in more physical (debatable - lets debate it here a bit)?
Please explicitly confirm non of the r-vs-rho and other normalization bugs need to be extended to delta_s. Are we tripple sure this one is being implemented correctly for getting a Jacobian independent psi width (note that this is easily testable by running with Hamada and PEST)?
Do we want this to be true? If the feature is the future and how we want users to default use the code, then shouldn't we be shipping examples that use it (i.e. psihigh 0.9995 for all DIII-D examples)? The DIII-D like example inherited the dmlim truncation behavior only to make it easy to compare to fortran, but we've moved past that I think (and the bechmark_with_fortran script can easily account for it by defaulting to a dmlim type comparison). @jhalpern30 and @matt-pharr may want to add their two cents. |
3 participants
Release note
grid_type = "auto"(harness @ 5fe6b42)LayerWidthsgainsdelta_norm,delta_viscoanddelta_dr; code constructing it positionally must supply them.sing_lim!andrefined_psi_gridgain optional keywords and are otherwise unchanged.GPEC can now locate the point past which adjacent rational surfaces' resistive layers overlap — the outer limit of the matched-asymptotic treatment, after Fitzpatrick, Nucl. Fusion 2025 Sect. 5.9 — and use it as an upper bound on the stability integration domain
psilim. The overlap point is recorded whether or not it binds, so a run always shows whether layer physics constrained the domain.psihighis untouched and remains user-set.Regression report
Four cases at the rebased head,
baf674659(base, #403) vs5fe6b4275, both sides run fresh with--forceon identical pinned environments (julia 1.11.6, manifest 7e5c34ad).The one
CHANGEDline is solver jitter, not this branch, and it is measured rather than asserted: re-runningbaf674659against itself with fresh samples moved γ at the same surface (3/1) by the same 0.1452 Hz, and5fe6b4275against itself by 0.1217 Hz at surface 2/1. The SLAYER root search is nondeterministic at the ~0.1 Hz (0.01 %) level per surface, in discrete steps, so the case's γ threshold sits below the solver's own reproducibility; that is a harness/threshold question (#418 and follow-up), not a change here. Every layer input feeding the root is bit-identical.Coverage gap, stated rather than hidden: the
edge_q_lawgate on the auto-grid edge floor is exercised bydiiid_n1anddiiid_slayer_n1(both auto grids, both unchanged) but by no limited deck, because every Solovev/LAR case usesgrid_type = "ldp"with fixedmpsi. That commit carries!as a precaution. A Solovev auto-grid harness case would close this.Notes for reviewers
This branch replaces the earlier design on this PR; the thread above records why. The architecture is the one agreed with @logan-nc: the overlap point is an upper bound on
psilim, applied as a cap onqlimat the top ofsing_lim!, sodmlim/qhighstill select the final surface from inside it and a cap beyondpsihighis inert with no special case. Nothing re-forms the equilibrium, which removes the trial/back-off loop,EdgeQuality, and the immutable-config collision from the previous version (~1800 insertions → ~700).Stacked on #403, which supplies
radial_label; base is set to that branch so the diff shows only these five commits, and will be retargeted todeveloponce #403 merges.Because the base is not
develop, the test workflow does not trigger on this PR while it is stacked. The full suite was run locally instead:julia -t 4 --project=. test/runtests.jlon treeb3d65ecf4, 25 testset groups, all passing, exit 0. CI will re-run it when the base is retargeted.Three things worth a close look:
LayerThickness.jl). The paper's τ_A carries no shear, whereas SLAYER'stau_hdivides byn·s, solu = (n|s|)·S_paperand Eq. (100)'s explicit(n|s|)^(-1/2)cancels exactly. The first implementation kept it and understated every width by√(n|s|).(q/|dq/dr|)^(1/2)— so it is only valid in the coordinate it was derived in, the toroidal-flux label of the paper's Eq. (30) (rs_method = :flux). On the DIII-D-like deck the label is worth one rational surface: cut at ψ = 0.99752 under:midplane, 0.99939 under:flux. The flux label reduces to the geometric minor radius on the circular LAR deck (r/a = 0.99–1.05) and departs from it by 20–45 % on the shaped deck.runtests_layer_overlap_fitzpatrick.jl) rebuilds the paper's Sect. 5.8 JET equilibrium and drives GPEC's ownslayer_parameters→slayer_layer_thickness, reproducing the paper's overlap boundary at Ψ = 0.9985111 (paper: 0.9985) for n = 1 and 0.9949738 (paper: 0.9952) for n = 4. It covers the exact chain where the bug in (1) lived.On the shipped DIII-D decks the cap binds nowhere: it sits beyond
psihighfor the ideal, riccati and gal_resistive decks, and on the SLAYER deck it binds by 1.1e-4 in ψ without crossing a rational. Two of those decks carry no kinetic profiles, so the scan cannot run on them at all. The two width channels (|δ_s| Riccati vs Eq. 100) disagree by about two surfaces; both are reported, and that disagreement is the subject of the ODE-vs-Eq.(100) reduction issue rather than this PR.