Skip to content

Tearing - BUGFIX! - Correct the coupled SLAYER determinant and Doppler-shift it by the kinetic-file E×B rotation - #463

Open
d-burg wants to merge 30 commits into
developfrom
bugfix/coupled-slayer-tauk-rescale-and-rotation
Open

d-burg wants to merge 30 commits into
developfrom
bugfix/coupled-slayer-tauk-rescale-and-rotation

Conversation

@d-burg

@d-burg d-burg commented Sep 16, 2026 •

Copy link
Copy Markdown
Collaborator

Release note

  • Audience: users
  • Numerical impact: SLAYER roots are now reported in the lab frame in both coupling modes. Uncoupled: each surface's Roots/Q_root and Roots/omega move by τ_k·n·Ω_E and n·Ω_E (on diiid_slayer_n1 ω goes −10.5 k → +31.0 k rad/s at the 2/1, −9.7 k → +21.0 k at the 3/1, −11.1 k → −2.8 k at the 4/1); γ is unchanged. Coupled (coupling_mode = "coupled"): the determinant is now Doppler-shifted per surface and rescaled in the right direction; on diiid_slayer_n1_coupled ω moves −10.5 k → +31.0 k rad/s and γ by −0.3 % (220.25 → 219.67 s⁻¹), inside the coupled root extraction's accuracy. Every layer input is bit-identical. (harness @ 546d9fb)
  • Migration: Roots/Q_root and Roots/omega are lab-frame in both modes; an uncoupled root's previous, E×B-frame value is ω − n·Ω_E from PerSurface/omega_E and PerSurface/n, and Diagnostics/* stay in each surface's E×B frame. A kinetic file without omega_E gives Ω_E = 0 (lab = E×B frame), records NaN in PerSurface/omega_E, and warns. In coupled mode widen Q_re_range if the run warns that Doppler offsets fall outside it. Interface removals: slayer_parameters(omega=…), multi_surface_coupling_full(rotation=…, ntor=…), KineticProfiles(omega_e=…, omega_i=…) and build_slayer_inputs(compute_omega_star=…); drop those arguments. KineticProfiles.omega is renamed omega_E.

The coupled tearing determinant evaluated every non-reference surface at the wrong Q (off by (τ_k/τ_ref)²), so its roots depended on which surface was chosen as reference. It also ignored the E×B rotation that the kinetic file already supplies. Both are fixed, and the root extraction is now reproducible run to run.

image

Left: on develop the coupled root sits at the 2/1's own-frame frequency (−10.5 krad/s). With this PR it moves by n·Ω_E at the 2/1 (41.5 krad/s) into the lab frame, landing where shifting the uncoupled 2/1 root predicts (30,993 against 30,996 rad/s); γ changes by under 1 % (221.9 → 219.7 s⁻¹ in that run). Right: the E×B Doppler shift n·Ω_E now applied at each surface, from the kinetic file. Both from the harness run at c4c95c4. With uncoupled roots now also reported in the lab frame, that comparison is in the raw output: the uncoupled 2/1 is at 30,996 rad/s and the coupled root at 30,993.

What changed

change where
Surface k's layer is evaluated at Q·τ_k/τ_ref, not Q·τ_ref/τ_k Coupled.jl, CoupledFullMatch.jl (one line each)
Real per-surface Doppler offset q_shift = −τ_k·n·Ω_E, from the kinetic file's omega_E, applied in coupled mode only SurfaceCoupling.jl, run_slayer.jl, TearingProblem.jl
The full 4m×4m form's + i·ntor·rotation[k] is removed: it shifted γ instead of ω and added s⁻¹ to a dimensionless Q CoupledFullMatch.jl
The AMR Delaunay triangulation is seeded GrowthRateExtraction.jl (3 lines)
PerSurface/omega_E and PerSurface/q_shift are written to the output; omega_E is NaN when the kinetic file has none Result.jl, HDF5Output.jl
Uncoupled roots are reported in the lab frame, Q + τ_k·n·Ω_E; each layer is still scanned and solved in its own E×B frame run_slayer.jl
Always-zero rotation and diamagnetic arguments are removed KineticProfiles.jl, LayerInputs.jl, LayerParameters.jl and eight test files
New harness case diiid_slayer_n1_coupled regression-harness/cases/

The source change is +169/−119 lines. Most of the 24 files are tests that only drop the removed arguments.

What moves

regress --refs 9a51446df,546d9fbec --force on feynman (SLURM, 8 threads), all 20 cases: develop @ 9a51446 against this branch @ 546d9fb, which has develop merged. Julia 1.11.9, manifest pinned.

case result
diiid_slayer_n1 (uncoupled): ω 2/1 −10,526 → +30,996; 3/1 −9,676 → +21,048; 4/1 −11,097 → −2,768 rad/s, each by n·Ω_E
diiid_slayer_n1 (uncoupled): γ unchanged (≤ 0.03 s⁻¹, develop's triangulation noise, below the harness threshold); layer inputs bit-identical
diiid_slayer_n1_coupled: ω −10,522 → +30,993 rad/s
diiid_slayer_n1_coupled: γ 220.25 → 219.67 s⁻¹ (−0.3 %)
diiid_slayer_n1_coupled: q_shift new output, absent on develop
the other 18 cases unchanged (gal_resistive_pe extracts nothing on either ref, as on develop)

Since the previous report (e7ce73c on SDCC), this branch adds the sign-anchor test, the NaN omega_E and the lab-frame uncoupled roots. Only the last moves a number, and only uncoupled ω and Q_root. The coupled γ came out −1.0 % (feynman, c4c95c4), −0.3 % (SDCC, e7ce73c) and −0.3 % here; see How far to trust the coupled γ below.

For the reviewer

  • Why τ_k/τ_ref. With Q = τ_k·(ω + iγ), one physical eigenvalue reaches surface k as Q·τ_k/τ_ref. With a diagonal Δ′ the coupled roots must then reproduce each surface's uncoupled root whichever surface is the reference; the old direction got only the reference surface right. Pinned by "Coupled roots do not depend on the reference surface" in runtests_dispersion_rotation.jl.
  • Why the factor n. GPEC's omega_E is per unit n, so the mode sees n·Ω_E.
  • Why this sign. The layer's Q is the mode frequency in the local E×B frame, so Q_layer = τ_k·ω_lab − τ_k·n·Ω_E. A test on two real SLAYER surfaces checks that the rotating determinant at Q + τ_ref·n·Ω_E equals the static one at Q to within 1e-4 (measured agreement 3e-6, the layer solve's rounding noise); the opposite sign is off by 27–394 %.
  • How far to trust the coupled γ. The coupled root is the raw contour intersection, with no polishing. Re-extracting one set of scan samples under different triangulations gave γ = 190.6 and 191.9 s⁻¹, a 0.65 % spread (measured before develop's Sauter default). The seed makes it repeatable, but it still shifts between environments: this branch gives 219.67 s⁻¹ on feynman (Julia 1.11.9) and 220.05 s⁻¹ on SDCC (Julia 1.12.6). The ω shift is resolved; the γ shift (−1.0 % in one run, −0.3 % in the other) is not.
  • Fortran has the same inversion on its slayer_growthrate branch, so the earlier Julia/Fortran coupled comparison did not validate the rescaling direction.
  • Not in this PR. A per-surface override of the file rotation for rotation scans (omega_E_kHz) is split into Tearing - FEATURE - Add a per-surface omega_E_kHz override of the kinetic-file E×B rotation #519, stacked on this one.
Full regression report (SLAYER cases)
Regression Report: diiid_slayer_n1
============================================================================================
Ref 1: 9a51446df  @ 9a51446df (2026-10-22)
       env: julia 1.11.9, x86_64-linux-gnu, manifest f826ee26 (pinned), 8 threads/8 BLAS, jit build
Ref 2: 546d9fbec  @ 546d9fbec (2026-10-23)
       env: julia 1.11.9, x86_64-linux-gnu, manifest f826ee26 (pinned), 8 threads/8 BLAS, jit build
--------------------------------------------------------------------------------------------
Quantity                            9a51446df  546d9fbec  Diff                 Status
--------------------------------------------------------------------------------------------
SLAYER surface indices              [6 elem]   [6 elem]   0.0e+00              OK
SLAYER poloidal m                   [6 elem]   [6 elem]   0.0e+00              OK
SLAYER toroidal n                   [6 elem]   [6 elem]   0.0e+00              OK
SLAYER minor radius rs              [6 elem]   [6 elem]   0.0e+00              OK
SLAYER r-based shear                [6 elem]   [6 elem]   0.0e+00              OK
SLAYER Lundquist S                  [6 elem]   [6 elem]   0.0e+00              OK
SLAYER D_norm                       [6 elem]   [6 elem]   0.0e+00              OK
SLAYER P_perp                       [6 elem]   [6 elem]   0.0e+00              OK
SLAYER tauk                         [6 elem]   [6 elem]   0.0e+00              OK
SLAYER iota_e                       [6 elem]   [6 elem]   0.0e+00              OK
SLAYER Q_e                          [6 elem]   [6 elem]   0.0e+00              OK
SLAYER Q_i                          [6 elem]   [6 elem]   0.0e+00              OK
SLAYER Q_root [2/1,3/1,4/1]         [3 elem]   [3 elem]   2.463e+00 (394.38%)  ** CHANGED **
SLAYER ω_Hz [2/1,3/1,4/1]           [3 elem]   [3 elem]   4.152e+04 (374.18%)  ** CHANGED **
SLAYER γ_Hz [2/1,3/1,4/1]           [3 elem]   [3 elem]   2.5e-02              OK
SLAYER no_root flags [2/1,3/1,4/1]  [3 elem]   [3 elem]   0.0e+00              OK
SLAYER enabled flag                 1          1          0.0e+00              OK
Runtime (s)                         298.6s     300.2s                          --
============================================================================================
Summary: 2 changed, 15 unchanged

Regression Report: diiid_slayer_n1_coupled
=====================================================================================
Ref 1: 9a51446df  @ 9a51446df (2026-10-22)
       env: julia 1.11.9, x86_64-linux-gnu, manifest f826ee26 (pinned), 8 threads/8 BLAS, jit build
Ref 2: 546d9fbec  @ 546d9fbec (2026-10-23)
       env: julia 1.11.9, x86_64-linux-gnu, manifest f826ee26 (pinned), 8 threads/8 BLAS, jit build
-------------------------------------------------------------------------------------
Quantity                     9a51446df  546d9fbec  Diff                 Status
-------------------------------------------------------------------------------------
SLAYER tauk                  [6 elem]   [6 elem]   0.0e+00              OK
SLAYER q_shift               N/A        [6 elem]   N/A                  N/A
SLAYER coupled Q_root        [1 elem]   [1 elem]   2.463e+00 (394.48%)  ** CHANGED **
SLAYER coupled ω             [1 elem]   [1 elem]   4.151e+04 (394.56%)  ** CHANGED **
SLAYER coupled γ             [1 elem]   [1 elem]   5.753e-01 (0.26%)    ** CHANGED **
SLAYER coupled no_root flag  [1 elem]   [1 elem]   0.0e+00              OK
SLAYER enabled flag          1          1          0.0e+00              OK
Runtime (s)                  310.1s     311.8s                          --
=====================================================================================
Summary: 3 changed, 3 unchanged, 1 missing/N/A

The other 18 cases report every quantity unchanged.

The run-to-run scatter in the coupled γ, and its fix
  • An earlier version of this PR saw the coupled γ take two values run to run with identical source.
  • Δ′, the energies, all per-surface SLAYER inputs and the stored AMR scan samples were bit-identical in every run; only the extracted roots differed.
  • The cause was the triangulation of the scan samples. DelaunayTriangulation inserts points in a random order drawn from Julia's global RNG, and an AMR grid is full of cocircular points, so the order picks the diagonals. Re-extracting the same stored samples under 10 seeds gave exactly the two observed values.
  • The fix draws the order from a fixed-seed generator. A test asserts that extraction no longer touches the global RNG; it fails without the fix.
Other notes
  • Rotation enters the determinant in coupled mode only. Each uncoupled layer is solved in its own plasma frame, and its root is then relabelled into the lab frame; PerSurface/q_shift stays 0 there. PerSurface/omega_E is written in both modes.
  • GGJ surfaces are not rotated, because they carry no time normalization; the run warns.
  • n = 0 with nonzero rotation is an error, since the shift would silently vanish.
  • External callers that build KineticProfiles with omega_e=zeros(...), omega_i=zeros(...) must drop those two arguments.
  • ForceFreeStates.ResistiveMatch(; rotation) is unaffected: it uses the Laplace-variable eigenvalue, where the imaginary sign is correct.

NO MERGE WITHOUT HUMAN REVIEW. This PR needs a third-party human reviewer's approval of the head that lands.

🤖 Generated with Claude Code

@d-burg d-burg self-assigned this Sep 16, 2026
@github-actions github-actions Bot added bugfix Something was wrong and now is not changed-results Results move or an interface breaks - read before upgrading labels Sep 16, 2026
d-burg and others added 10 commits September 23, 2026 16:31
…e inner-layer Q argument

The coupled determinant evaluated every surface's inner layer at the same
lab-frame frequency, so there was no way to impose differential rotation
between rational surfaces and no way to test the coupled growth rate for
rotation invariance.

Add a real `q_shift` to `SurfaceCoupling`, applied to the layer's Q argument
in the scalar residual, the reduced m x m determinant, and the 4m x 4m
Pletzer-Dewar matching. Re(Q) is the lab-frame frequency in each surface's own
normalization, so the offset is real and Dopplers that surface away from the
common eigenvalue; zero reproduces the previous static behaviour exactly.

The runner builds the shifts from a new `[SLAYER] omega_shift_kHz` vector
(per-surface lab-frame offsets in kHz, core to edge), converting with each
surface's tau_k, and records the applied dRe(Q) at `PerSurface/q_shift`.

On the shipped DIII-D-like deck in coupled mode, a 3 kHz shift on the 3/1
alone moves the n=1 coupled growth rate from 205.6 to 233.4 1/s (+13.5%),
with omega essentially unchanged.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
(cherry picked from commit bdafb63dbe9334232bcead6357e7ed9190cb5835)
…he direct tauk rescale

Verified the layer normalization against Fitzpatrick's TJ
(Documentation/Layer.tex and Layer/Layer.cpp). TJ solves for the layer
eigenvalue as ghat = i(Q_E - omega*tau_k) with Q_E = tau_k*omega_E, and
reports omega = (Q_E - Im(ghat))/tau_k. Our Riccati bridge already documents
ghat = -i*Q, so Q_k = tau_k*(omega - omega_E_k).

Two consequences:

Rename `omega_shift_kHz` to `omega_E_kHz` and negate the conversion, so the
input is the surface's E*B rotation frequency and the layer responds to the
mode frequency in its own fluid frame, as in TJ. The knob was new in the
previous commit, so nothing downstream depended on the old spelling.

Add `tauk_rescale`, selecting the inter-surface Q normalization used by the
coupled determinant. The default `:legacy` keeps Q*tauk_ref/tauk_k, which
every coupled growth rate produced by this code so far has used. `:direct`
uses Q*tauk_k/tauk_ref, the direction implied both by our own extraction
convention omega = Re(Q)/tauk and by TJ. The two differ by (tauk_k/tauk_ref)^2
and move results, so the correction is opt-in rather than applied silently.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
(cherry picked from commit d1c4e6766ffe958a548a61deb47627728a25588a)
…in the coupled determinant

The coupled determinant evaluated surface k's inner layer at
Q_k = Q * (tauk_ref / tauk_k). That direction is inverted.

Q is defined as tauk * omega -- Fortran GPEC slayer/params.f sets
`Q = Qconv*omega` with `Qconv = lu**(1/3)*tau_h`, this code builds Q_e and
Q_i the same way, and the growth-rate extraction inverts it as
omega = Re(Q)/tauk. One shared physical eigenvalue therefore reaches surface
k by MULTIPLYING by tauk_k, giving Q_k = Q * (tauk_k / tauk_ref).

The inconsistency is visible inside a single layer solve: the scanned
frequency was scaled as 1/tauk_k while that surface's Q_e/Q_i scale as
tauk_k, yet the layer equations add them (ghat + i*Q_e). The two were
mis-scaled by (tauk_k/tauk_ref)^2 -- a factor of 18 at the 3/1 on the TJ
pc=0.1 equilibrium.

Confirmed three ways: the definitional argument above; Fitzpatrick's TJ,
where ghat_k = tauk_k * g_k (TJ Documentation/Layer.tex); and a rigid-rotation
invariance test, which the corrected direction satisfies to four significant
figures in omega while the old one misses it by 15x and swings gamma by 44%.

`tauk_rescale` now defaults to `:direct`. `:legacy` restores the old
behaviour for reproducing results generated before this correction. The
existing rescale test only restated the implementation, so it is rewritten to
pin the physics and to cover both directions.

NOTE: Fortran GPEC has the same inversion on its `slayer_growthrate` branch,
slayer/growthrates.f: `g_tmp = (g_in*sl_in%Qconv_arr(1))/tauk`. This Julia
code was a faithful port; the upstream fix is still outstanding.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
(cherry picked from commit d352ae24373a0052f44277b7d63cb460ab7dbaa5)
…rom slayer_parameters

slayer_parameters required an `omega` rotation keyword but never read it, and
build_slayer_inputs dutifully passed the kinetic file's rotation into it. A single
layer is solved in its own plasma frame, so rotation has no place in the per-surface
parameters; it belongs in the coupled determinant, where surfaces rotating at
different rates must share one lab-frame eigenvalue. Drop the argument.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…e imaginary rotation term in the full coupled determinant

The reduced determinant was corrected to evaluate surface k at Q·tauk_k/tauk_ref,
but the full 4m×4m Pletzer-Dewar form still used the inverted tauk_ref/tauk_k and
still added i·ntor·rotation[k] alongside the new real q_shift. That imaginary term
is Fortran's guess_modify carried over without the change of convention: with
Q = tauk·(ω + iγ) it moves γ rather than ω, and it adds a frequency in s⁻¹ to a
dimensionless Q.

Give the full form the same tauk_rescale switch (default :direct) and take the
Doppler offset solely from each SurfaceCoupling's real q_shift. The rotation and
ntor fields are removed. Tests pin the per-surface Q argument under both
rescalings and that the shift moves Re(Q) only.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…y it only to the coupled determinant

run_slayer loaded the kinetic file's omega_E and then discarded it; the only way
to rotate a surface was a hand-typed omega_E_kHz vector. That vector was also
applied to uncoupled scans, where a layer is already solved in its own plasma
frame and a shift only relabels the reported frequency, and it omitted the
toroidal mode number: the file's omega_E is per unit n, while the layer sees
n·Ω_E (TJ's ω_E = (m/r_s)·V_E).

- Sample Ω_E from the kinetic file at each rational surface by default; a
  non-empty omega_E_kHz overrides it in the same per-unit-n convention.
- Doppler-shift only the coupled determinant, by ΔRe(Q_k) = −τ_k·n·Ω_E,k.
  Uncoupled runs apply no shift. GGJ surfaces carry no time normalization and
  are not shifted.
- Record the shifts actually applied in SLAYERResult.q_shift and write those to
  PerSurface/q_shift, instead of recomputing from the control block.
- Warn when a surface's offset in the scanned Q falls outside Q_re_range, since
  the lab-frame root then silently becomes :no_root.

The sign needs no mapping onto TJ's conventions. On the DIII-D-like deck the
static uncoupled roots sit on ω_*e = −Q_e/τ_k (0.1%, 0.04%, 0.4% at the 2/1, 3/1,
4/1), and the file's omega_E is built from the same ion diamagnetic formula as
SLAYER's ω_*i, so ω, ω_* and Ω_E share one convention and the Doppler term is a
plain subtraction. In coupled mode the deck's root then lands at ω = 30993 rad/s
against the predicted n·Ω_E + ω_*e = 30986 at the 2/1, with γ = 195.0 s⁻¹ against
the uncoupled 194.8.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
No harness case exercised the coupled determinant, so neither the inter-surface
Q normalization nor the E×B Doppler shift was guarded. Reuse the n=1 SLAYER deck
in coupled mode, widening Re(Q) to [-4, 6] because the kinetic-file rotation
moves the root to Re(Q) ≈ 2.06, outside the deck's [-2, 2]. The run adds about
17 s of SLAYER time on top of the shared equilibrium and stability stages.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…s from KineticProfiles

KineticProfiles carried omega_e and omega_i, but nothing could put a real value
in them: run_slayer filled them with zeros, the kinetic file schema has no ω_*
column, and compute_omega_star defaulted to true with no caller in src/ ever
passing it, so build_slayer_inputs always derived ω_*e and ω_*i from the
density and temperature splines instead. The fields, the zeros and the
compute_omega_star branch are removed; the container keeps the profiles that
are actually read, including omega, which now carries the E×B rotation into the
coupled determinant.

The two tests that used the pass-through to pin Q_e and Q_i now assert the
derived values instead: the fixtures are linear in ψ at constant n_e, so
ω_*e = -700/psio and ω_*i = +600/psio exactly.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…and tidy the rotation plumbing

- Drop `tauk_rescale` entirely: `Q·tauk_k/tauk_ref` is the only normalization, so the
  `:legacy` reproduction shim (which could not reproduce pre-correction results once the
  kinetic file carries omega_E) goes, with its field, keyword, TOML key and validation.
- The example deck's rotation note gives the shift in the scanned plane (+τ_ref·n·Ω_E, the
  same for every surface) rather than each surface's own q_shift.
- `Roots/omega` documents its frame: E×B per surface when uncoupled, lab when coupled.
- GGJ surfaces take no q_shift (they have no time normalization); Ω_E is resolved once per run.
- Replace the documentation-file citations with the relation itself and trim the comments
  that narrated the old bug.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…d roots are reproducible

`_extract_growth_rates_amr` triangulated the AMR samples with DelaunayTriangulation's default
randomized insertion order, drawn from the global RNG, which Julia seeds from entropy at startup.
An AMR grid is full of cocircular points, where the Delaunay triangulation is not unique, so each
run could pick different diagonals and hence different contour crossings. Polished roots
converge past that; the unpolished coupled root did not. On the DIII-D-like coupled deck, the
identical scan samples gave γ = 190.646 or 191.888 s⁻¹ depending only on the seed (10 seeds,
5/5 split). This is the run-to-run scatter this PR reported. Threads were not the cause.

The triangulation now draws its insertion order from a fixed-seed Xoshiro, and a test asserts
that extraction leaves the global RNG untouched. The coupled root is now reproducible, but it
still carries the triangulation's ~0.7% interpolation error until coupled root polishing exists.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
d-burg and others added 7 commits September 25, 2026 11:04
…and rescale test

- Throw when a surface has E×B rotation but SLAYERParameters.n < 1, instead of
  silently zeroing the Doppler offset.
- Put the toroidal mode number back into the documented q_shift = −τ_k·n·Ω_E.
- State the omega_E_kHz units correctly (per unit n like the kinetic file's
  omega_E, but kHz rather than rad/s) and drop the duplicated inline comment.
- Define the coupled reference surface once and word the Doppler-window check
  as a heuristic; drop the per-run Doppler @info (q_shift is in the output).
- Drop deck-specific numbers from the coupled-root comment and example deck;
  note that coupled runs need a wider Q_re_range.
- Fix the known-coupled-root test to use the tauk_k/tauk_ref direction and a
  second surface with a different tauk, so it exercises the rescale.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…ga to omega_E

The field and keyword now carry the same name as the kinetic-file dataset they are read from, Ω_E per unit n in rad/s.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…otation per surface

omega_E_kHz is now a TOML table keyed by the surface's m/n, e.g. {"2/1" = 0.0, "3/1" = 3.0}; unlisted surfaces keep the kinetic file's Ω_E and an m/n matching no analysed surface is an error. The positional list form is removed.

The Ω_E actually resolved on each surface is carried as SLAYERResult.omega_E and written to Tearing/PerSurface/omega_E [rad/s, per unit n] in both coupling modes; PerSurface/q_shift stays the applied Doppler offset (zero when uncoupled).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
State the lab-to-E×B-frame relation behind q_shift = −τ_k·n·Ω_E in the _q_shift docstring, and test on two real SLAYER surfaces that the rotating coupled determinant at Q + τ_ref·n·Ω_E equals the static one at Q, so every root keeps its γ and its lab frequency moves by +n·Ω_E.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…s rounding noise

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@d-burg d-burg changed the title Tearing - BUGFIX! - Correct the coupled SLAYER determinant and Doppler-shift it by the kinetic-file E×B rotation Tearing - BUGFIX! - 🚨 Correct the coupled SLAYER determinant and Doppler-shift it by the kinetic-file E×B rotation Oct 5, 2026
d-burg and others added 4 commits October 5, 2026 15:50
…yer-tauk-rescale-and-rotation

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…tity in the coupled SLAYER case

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…nd convention statements

Rotation is tested once, through run_slayer_from_inputs, with the n = 2 and
wrong-length cases moved there. omega_E_kHz keys are validated in one place,
the run-time check against the analysed surfaces. The Doppler convention is
derived once, in _q_shift.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…e coupled-determinant fix

The fix takes E×B rotation from the kinetic file only. The per-surface
override is a rotation-scan feature and goes to its own pull request.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
d-burg and others added 2 commits October 19, 2026 15:00
Brings in the post-solve TearingProblem API and the TOML renames. The E×B rotation
per surface now feeds run_slayer_from_inputs from TearingProblem's solve, since the
run_slayer pipeline it used to come from is gone. The overlap-scan test from develop
drops the toroidal rotation and diamagnetic fields that this branch removed from
SLAYERParameters and KineticProfiles.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@d-burg
d-burg requested a review from matt-pharr October 21, 2026 01:57
@matt-pharr

Copy link
Copy Markdown
Collaborator

@d-burg you requested I review this today so I will assume it is ready for review. I will review it and leave it to you to edit after that

@matt-pharr
matt-pharr marked this pull request as ready for review October 21, 2026 15:34
@d-burg d-burg changed the title Tearing - BUGFIX! - 🚨 Correct the coupled SLAYER determinant and Doppler-shift it by the kinetic-file E×B rotation Tearing - BUGFIX! - Correct the coupled SLAYER determinant and Doppler-shift it by the kinetic-file E×B rotation Oct 22, 2026
@matt-pharr
matt-pharr added this pull request to stack #562 October 23, 2026 07:25
matt-pharr and others added 3 commits October 23, 2026 09:31
…-fluid zero

The rotating-equals-static determinant check holds identically for q_shift = −τ_k·n·Ω_E, so it
pins how _q_shift composes with the determinant but not the physical sign. Add a test that the
layer Δ vanishes at Re(Q) = +τ_k·ω_*e, tying the scanned Re(Q) to the ζ axis of GPEC's ω_*e and
the kinetic file's Ω_E. Reword the coupled harness case's threshold comment.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
_load_profiles reports whether the file carries omega_E, and the tearing solve passes no rotation
when it does not. No shift is applied either way, but result.omega_E and PerSurface/omega_E now
record NaN instead of zeros, so a missing profile is distinguishable from a measured zero rotation.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@matt-pharr

Copy link
Copy Markdown
Collaborator

I pushed two small commits on top of the review:

The Doppler sign test now has a physical anchor (test/runtests_dispersion_rotation.jl)

  • The "Real SLAYER surfaces" test checks that the rotating determinant at Q + τ_ref·n·Ω_E equals the static one at Q. With q_shift = −τ_k·n·Ω_E that holds identically. So it pins how _q_shift composes with the determinant, but it would also pass with the sign flipped throughout. Its comments now say so, and its inputs use GPEC-signed ω_* (ω_*e < 0, ω_*i > 0).
  • New testset, "Scanned Re(Q) runs along the ζ axis of Ω_E". The Doppler sign is right only if the scanned Re(Q) lies along the same ζ axis as the kinetic file's Ω_E and GPEC's ω_*e. The layer Δ vanishes where the mode co-rotates with the electron fluid. The test checks that this zero sits at Re(Q) = +τ_k·ω_*e, not −τ_k·ω_*e: the |Δ| ratio is about 0.015 at Im Q = 0.01. Reversing either Q_e = −τ_k·ω_*e or the i·conj(Q) bridge in solve_inner moves the zero to −τ_k·ω_*e and fails the test.

omega_E is NaN when the kinetic file has none (run_slayer.jl, TearingProblem.jl, Result.jl, HDF5Output.jl)

  • _load_profiles now reports whether the file carries omega_E, and the tearing solve passes no rotation when it doesn't. That still applies no shift (q_shift stays 0). But result.omega_E and PerSurface/omega_E now record NaN, the HDF5 not-computed sentinel, so a missing profile can be told apart from a measured zero rotation. The same holds for run_slayer_from_inputs called without omega_E. A new check in runtests_slayer_runner.jl covers it.

Harness case comment (regression-harness/cases/diiid_slayer_n1_coupled.toml)

  • The comment above the coupled-root thresholds now explains why they are tight: the seeded triangulation repeats exactly in one environment. It also says the 0.65 % triangulation spread applies across environments. The unfinished "ω and Q need the same measurement" line is gone, and the thresholds are unchanged.

No harness number moves. The NaN only replaces zeros in an output that applies no shift; everything else is tests and comments.

@matt-pharr matt-pharr left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@d-burg the physics checks out. The τ_k/τ_ref direction is right, and so is the Doppler sign: the two sign flips (Q_e = −τ_k·ω_*e and the i·conj(Q) bridge in solve_inner) cancel, so Re(Q) runs along GPEC's ζ axis, the same axis as Ω_E. The new sign-anchor test pins that. The harness at e7ce73c moves only the coupled-SLAYER results, as described.

One design question before I approve: which frame should Roots/omega and Roots/Q_root report in?

Right now it depends on the mode: uncoupled roots are in each surface's own E×B frame, and the coupled root is in the lab frame. I'd prefer always lab frame:

  • The coupled root has no single E×B frame, because each surface rotates at its own Ω_E. So the lab frame is the only frame both modes share. Reporting the coupled root in the reference surface's frame instead would bring back the reference-surface dependence this PR removes.
  • Uncoupled and coupled results become directly comparable. The figure shifts the uncoupled 2/1 root by n·Ω_E by hand (30,996 against 30,993 rad/s); with lab-frame output, the raw data shows that comparison directly.
  • It's what experiments measure. Mirnov coils see lab-frame frequencies.
  • For uncoupled roots it is only a relabel. Report Q_root + τ_k·n·Ω_E,k and the matching ω, so ω = Re(Q)/τ_k still holds. The scan and root-finding stay in each layer's own frame, and γ doesn't change. The E×B-frame value stays recoverable as ω − n·Ω_E from PerSurface/omega_E and PerSurface/n. When the kinetic file has no omega_E, take Ω_E = 0, as the coupled path already does.

The cost is that uncoupled ω moves, by n·Ω_E at each surface (+41.5 krad/s at the 2/1 in diiid_slayer_n1). So the release note would need updating, plus a fresh harness run.

Could you either:

  1. switch to always reporting in the lab frame, as above, or
  2. tell me you feel strongly that uncoupled roots should stay in their own E×B frame, and why? Then I'll approve as is, and we'll add the frame to the Roots/Q_root description, which doesn't mention it yet.

d-burg added a commit to PrincetonUniversity/GPEC that referenced this pull request Oct 23, 2026
…d determinant

dispersion_det evaluated surface k's layer at g*tau_1/tau_k. With
g = tau_1*(gamma + i*omega) the same physical eigenvalue reaches surface
k as g*tau_k/tau_1, which is the normalisation its Q_e and Q_i carry.
The old direction got only the reference surface right, so every other
surface's branch of the coupled determinant sat at the wrong frequency
and growth rate.

Check on the DIII-D-like example with a real, diagonal Delta' (where
the coupled roots must equal the uncoupled ones): the 3/1 and 4/1 roots
now land on Q_k*tau_1/tau_k to 7 digits. Before, the 3/1 branch sat at
Re(Q) = -0.514 instead of -1.062.

Same defect and fix as OpenFUSIONToolkit/GPEC#463.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
d-burg added a commit to PrincetonUniversity/GPEC that referenced this pull request Oct 23, 2026
…xB rotation

The coupled determinant evaluated every surface's layer at the same
frame frequency, as if the surfaces rotated together. The layer's Q is
the mode frequency in the local ExB frame, so with one lab-frame
eigenvalue surface k is evaluated at

   Q_k = Q*tau_k/tau_1 - tau_k*n*Omega_E,k

using the ExB rotation already read from the kinetic file. Applied in
coupled mode only (rotation_flag, default true): an uncoupled layer is
solved in its own frame, where the shift only relabels the frequency.
The coupled root is therefore reported in the lab frame.

The scan window gets a Re(Q) centre (scan_center). Left unset, a
coupled scan with rotation is centred on the Doppler-shifted branches,
and a warning is printed when they do not all fit in the window.
The uncoupled scan warns when -Q_e lies outside the window.
Q_shift is written to the NetCDF file.

Check with a real, diagonal Delta' on the DIII-D-like example: the
coupled roots equal the uncoupled roots shifted by tau_1*n*Omega_E,k
to 6 digits. With the full matrix the 2/1 root moves from -10.52 to
+31.00 krad/s, with the growth rate unchanged at 79.5 1/s.

Follows OpenFUSIONToolkit/GPEC#463.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@matt-pharr
matt-pharr force-pushed the bugfix/coupled-slayer-tauk-rescale-and-rotation branch from a5ab588 to 9e95d23 Compare October 23, 2026 17:33
…s the coupled root is

Uncoupled roots were reported in each surface's own E×B frame and the coupled root in the lab
frame. The coupled root has no single E×B frame, so the lab frame is the one both modes share.
Each uncoupled layer is still scanned and solved in its own frame; its Roots/Q_root and
Roots/omega now carry τ_k·n·Ω_E and n·Ω_E, so γ is unchanged and the E×B-frame value is
ω − n·Ω_E. Where the kinetic file has no omega_E the shift is zero, and the tearing solve warns.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@matt-pharr

Copy link
Copy Markdown
Collaborator

@d-burg sorry I messed up my worktrees and accidentally pushed something here, I force pushed back over it to undo it but please beware of if I messed anything up on this branch

@d-burg

d-burg commented Oct 26, 2026

Copy link
Copy Markdown
Collaborator Author

@matt-pharr Lab frame everywhere, as you suggested, in 546d9fbec:

  • Uncoupled roots. Each layer is still scanned and solved in its own E×B frame. Its root is reported as Q + τ_k·n·Ω_E and ω + n·Ω_E, so γ is unchanged and the E×B-frame value is ω − n·Ω_E. Diagnostics/* (valid roots, poles, filtered roots) stay in each surface's E×B frame, and the docstrings and HDF5 descriptions now say which frame each output is in.
  • No omega_E in the kinetic file. Ω_E is taken as 0, so lab and E×B frame coincide. PerSurface/omega_E stays NaN, as in your commit, and the tearing solve now warns that the roots assume Ω_E = 0.
  • Harness (feynman, all 20 cases, develop 9a51446df vs 546d9fbec): diiid_slayer_n1 ω moves by n·Ω_E (2/1 −10,526 → +30,996, 3/1 −9,676 → +21,048, 4/1 −11,097 → −2,768 rad/s) with γ unchanged. The coupled case is as before (ω +30,993, γ −0.3 %), and the other 18 cases are unchanged. The uncoupled 2/1 and the coupled root now agree in the raw output, 30,996 against 30,993.
  • New test. The uncoupled Q_root moves by exactly τ_k·n·Ω_E, while γ and the layer-frame extraction stay identical.

Release note and description are updated. Thanks for merging develop (cbfa58d5e); since that touches src/, I'm re-running the harness there and will update the stamp. Your accidental push didn't cost anything: the branch was back at 9e95d2395 when I built on it.

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

bugfix Something was wrong and now is not changed-results Results move or an interface breaks - read before upgrading

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants