PE toroidal torque changes sign under grid refinement while its magnitude stays stable to ~1%

[d-burg]

Summary

perturbed_equilibrium/energies/toroidal_torque on the DIII-D-like ideal example changes sign under small, benign changes to the discretization, while its magnitude stays stable to about 1%. Radial refinement, poloidal refinement, and (independently reported) a dependency update all flip it. A physical torque should not behave this way: the sign is being set by something arbitrary rather than by the physics.

Measured after #350, so the coil-forcing fixes in that PR (b2ec04c3, 8dab0ed9) do not resolve it.

Evidence

All runs use examples/DIIID-like_ideal_example unchanged except for the grid controls named.

Radial refinement on a fixed ldp grid, post-#350 (grid_type="ldp", mpsi as shown; mpsi=0 is the deck's auto grid):

grid	torque
auto (561 knots)	+0.050865
ldp 512	+0.050892
ldp 1024	−0.050645
ldp 2048	+0.050876

Same scan before #350:

mpsi	128	256	512	1024	2048
torque	−0.0507289	−0.0506007	+0.0508128	−0.0510857	+0.0509140

Poloidal refinement at fixed mpsi = 1024 (pre-#350) — changing only mtheta:

mtheta	torque
256	−0.0510857
512	+0.0513171

Across every configuration above the magnitude sits in 0.0506–0.0513, a spread of ~1.4%, while the sign takes both values.

It is deterministic, not noise

Re-running mpsi = 1024 and 2048 with the Euler-Lagrange tolerance tightened from 1e-10 to 1e-12 reproduces each value — including its sign — to 1.8e-10 or better:

quantity	mpsi	el_tol=1e-10	el_tol=1e-12	rel diff
pe_torque	1024	−0.051085748	−0.051085748	4.6e-11
pe_torque	2048	+0.0509139592	+0.0509139592	1.8e-10

So the sign is a reproducible function of the grid, not run-to-run instability, and not integration error.

Independent report of the same symptom

In #338, @matt-pharr reported a set of package updates that "flip the sign of the torque in the d3d ideal example". That is a third, unrelated trigger producing the same signature — sign inverts, magnitude survives.

Why this looks like a convention rather than a convergence problem

A quantity that is converging would show its magnitude settling and its sign fixed. Here the magnitude is already stable to ~1% across an 8× radial refinement, a 2× poloidal refinement, and a code change — while the sign is not stable under any of them. That pattern fits a sign carried by something with an arbitrary choice in it (an eigenvector phase or normalization convention, or a branch selection somewhere in the torque assembly) rather than a physical quantity that has not yet converged.

Worth checking whether the torque is computed from an expression that is genuinely phase-invariant. Quantities built as Im(a* b) are invariant to a global phase on the mode; anything effectively linear in the eigenvector's sign is not.

Why it matters now

The regression harness tracks this quantity (pe_toroidal_torque in regression-harness/cases/diiid_n1.toml, noise_threshold = 1e-12). Any run that lands on the other sign reports a ~200% regression — which is exactly the 200.47% torque "regression" quoted in #338's original table. Until this is resolved the quantity cannot carry a golden value, and it will keep generating false regressions whose cause is invisible to whoever is reading the report.

Reproducing

# in a copy of examples/DIIID-like_ideal_example, set Equilibrium.grid_type="ldp"
# and Equilibrium.mpsi to each of 512, 1024, 2048, then:
julia --project=. -e 'using GeneralizedPerturbedEquilibrium; GeneralizedPerturbedEquilibrium.main(["<deck>"])'
# read perturbed_equilibrium/energies/toroidal_torque from gpec.h5

[matt-pharr]

@d-burg the toroidal torque is an error, it should technically be 0 for ideal cases (so it is very small at 0.05)

[d-burg]

@matt-pharr thanks — that reframes this usefully, and the definition agrees with you: toroidal_torque = -2·n·Im(⟨Φ_tot, Λ⁻¹·Φ_tot⟩/4) is the imaginary part of a quadratic form, identically zero for Hermitian Λ. So the nonzero value is measuring non-Hermiticity, and an arbitrary sign is exactly what noise about zero would look like.

I went looking for the discretization knob that controls it, expecting it to shrink under refinement. It doesn't shrink under any of them. All runs on examples/DIIID-like_ideal_example, one knob at a time:

axis	values	pe_torque
mpsi (ldp)	128 → 2048	−0.050729, −0.050601, +0.050813, −0.051086, +0.050914
mtheta @ mpsi=1024	256 → 512	−0.051086 → +0.051317
mthvac (auto grid)	512 → 1024 → 2048	+0.05086535, +0.05086513, +0.05086502
eulerlagrange_tolerance	1e-10 → 1e-12	identical to 1.8e-10

The vacuum-grid row is the one I expected to move, given the Wᵛ is non-Hermitian above tolerance 5.6e-4 > 1e-4 … Increase vacuum grid resolution warning — a 4× refinement changes |torque| by 4e-6 relative, i.e. not at all.

So across an 8× radial refinement, a 2× poloidal refinement, a 4× vacuum refinement and a 100× tighter ODE tolerance, the magnitude sits at 0.0509 ± 1%. The sign varies with the radial/poloidal grid but not with the vacuum grid or the tolerance. On a given grid it is bit-reproducible (two BLAS thread counts agree to 1.5e-14).

That is the part I can't reconcile with "numerical error": a discretization residual should decrease when the discretization is refined, and this one is flat. For scale it is ~4–6% of the energy quantities on the same run (ep ≈ −1.35, ev ≈ 2.15, et ≈ 0.80). Is a systematic, refinement-insensitive residual at that level what you'd expect from the Λ⁻¹ assembly, or does it suggest something is contributing an anti-Hermitian part that isn't a discretization artifact?

Practical consequence for the golden-value work either way: the torque will be pinned as a bound (|T| < measured atol, target 0) rather than as a value, since pinning a number that should be zero would encode the bug as truth.

[matt-pharr]

@d-burg the working hypothesis was that it is due to some small amount of contamination of the small-large solutions near the rational surfaces. Large solutions in the vicinity of the rational should be 0 in ideal cases, and this rule failing can cause an island which gives rise to a torque and non-hermitian-ness in W.

This is something that we have inherited from the fortran too, you should try running a fortran case.

[d-burg]

@matt-pharr thanks — the data suggest the ca route can't answer this. sing_get_ca solves
ua·ca = u, and since EulerLagrange.jl:639 introduces the small resonant solution at
ua[:, ipert_res + numpert_total, :],
ca[:,:,1] holds the large-solution coefficients — so the resonant row of ca_left/ca_right is
exactly the "large solution should be 0" test. Measured on the DIII-D-like ideal deck, the resonant
large/small ratio is 1e-4 to 1e-9. That looks like residual contamination until you check the
precision available: the ca blocks span 45–58 orders of magnitude, and against that roundoff
floor two of the five surfaces sit below it and two more are within 3×:

surface	‖ca_large‖	roundoff floor
1	2.88e30	1.7e29
2	3.54e35	2.7e35
3	2.14e39	7.8e38
4	1.26e38	2.5e40
5	2.78e39	7.7e41

The forward integrator's columns grow exponentially from the axis, so the asymptotic decomposition
has no precision left to resolve a small residual. Consistent with zero, but genuinely unresolvable —
I don't think this diagnostic can confirm or refute the hypothesis, and it's worth someone not
repeating it.

A fifth axis, which the earlier four missed. EulerLagrange.jl:533 stops the sweep at
psi_end = sing.psifac - singfac_min/|n·q1|, and line 630 crosses with a fixed 2-step trapezoidal
extrapolation u .+= (du1 + du2)·dpsi. singfac_min is the only knob setting that jump width — and
mpsi, mtheta, mthvac and eulerlagrange_tolerance all leave it untouched, which would explain
why the torque ignored every one of them.

Scanning it over three decades (deck otherwise byte-identical):

singfac_min	pe_torque	Δ vs previous	et[1]	|torque|/|et|
1e-3	+5.1988937e-02	—	0.822767	6.319%
1e-4	+5.0865351e-02	−2.16%	0.801371	6.347%
1e-5	+5.0767211e-02	−0.193%	0.798981	6.354%
1e-6	+5.0758427e-02	−0.017%	0.798726	6.355%

The torque converges rather than vanishing. Each 10× narrower jump cuts the change ~11× — clean
first-order convergence in singfac_min — to a Richardson-extrapolated limit of 5.0758e-2. If the
crossing extrapolation were manufacturing the torque, shrinking dpsi by 1000× would have killed it;
instead it settles onto a finite nonzero value, at a fixed 6.355% of the free-boundary energy (et[1]
converges the same way, to ≈0.79870).

So across all five axes — mpsi (8×), mtheta (2×), mthvac (4×), eulerlagrange_tolerance (100×)
and now singfac_min (1000×) — the torque is a converged property of the formulation, not a
resolution or jump artifact. That is the part I find hard to reconcile with a discretization error,
and it does point at something structural in the way the solution is continued across the rationals,
which is your reading. It just isn't reachable by refining anything.

On your Fortran suggestion — now the obvious next step, since refinement is exhausted: upstream
develop now ships TkMkr_D3Dlike_Hmode.geqdsk in docs/examples/run_kinetic_example/,
byte-identical to our deck's (sha256 9449e54b8488…), so the comparison can be
same-equilibrium rather than confounded by different inputs. If Fortran shows the same converged
~6.4%-of-energy residual, that settles it as inherited rather than a porting artifact.

Golden treatment is unchanged either way: pinned as a bound (|T| < atol, target 0), not a value.

[d-burg]

Ran the Fortran comparison. It does not reproduce the torque — this looks like a Julia-side
defect in Im(py), not something inherited.

Fortran GPEC rebuilt from the exact checkout (v1.5.7-611-ge68d7ac2, gfortran 15.2), example inputs
unmodified. The decks match closely: byte-identical TkMkr_D3Dlike_Hmode.geqdsk, same d3d C-coil
currents [1000, 500, −500, −1000, −500, 500], hamada, nn=1, singfac_min=1e-4.

Both codes compute the same quadratic form — Julia Response.jl:98 is a direct port of
gpout.f:1353, py = ⟨Φ_resp, L_plasma⁻¹·Φ_resp⟩/4, with the energy as Re(py) and the torque as
−2n·Im(py):

	Fortran	Julia
Re(py) — energy to perturb plasma	3.423	3.423266	0.008%
energy to perturb vacuum	5.955	5.839	1.9%
amplification factor	1.740	1.706	2.0%
−2n·Im(py) — torque	5.564e-05	5.0758e-02	912×

The real part of the same complex number agrees to 8 parts in 100,000 — so the whole chain feeding
it (equilibrium, DCON solve, coil forcing, inductances, response flux) is faithfully ported. The
imaginary part is off by ~900×. Fortran's torque sits at 1.6e-5 of its energy scale, i.e. zero to
numerical precision, exactly as the Hermitian argument demands; Julia's sits at 1.5e-2.

An older Fortran build (v1.5.5-323) gave 8.1e-5, so this is stable across ~290 commits of Fortran
history rather than a recent change.

That also resolves why five refinement axes never touched it: it was never a discretization error.
Something in the Julia Λ⁻¹ / inductance assembly is producing an anti-Hermitian part the Fortran
does not have. Worth noting the free-boundary energy carries a small imaginary part in Fortran too
(DCON reports imaginary = 1.4e-04 against real = 7.5e-01), so the mechanism exists there — it is
just ~1000× smaller.

Correction to my previous comment: I normalized the torque by et[1], the DCON free-boundary
eigenvalue, giving 6.355%. The physically matched denominator is Re(py) from the same quadratic
form, which gives 1.48%. The singfac_min convergence result is unaffected — only the normalization
was wrong.

Golden treatment unchanged, and now clearly right: |T| < atol with target 0, since Fortran shows
the converged value should be ~1e-5 relative, not 1.5e-2.

Happy to open a separate issue for the Im(py) assembly if that is easier to track than this thread.

[d-burg]

Split the root cause out to #401 — Im(py) is ~900× larger than Fortran GPEC's, with Re(py) agreeing to 0.008%. This issue's original symptom (sign flips under grid refinement) is explained by that: it was never a discretization effect, so there is nothing to converge.

Most specific divergence found so far is that Fortran Hermitian-symmetrizes the surface inductance (gpresp.f:284) and the Julia port does not (ResponseMatrices.jl:306) — untested as a fix, detailed in #401.

[matt-pharr]

Is claude missing an hermitianization step that fortran is doing? Fortran hermitianizes a number of its matrices artificially.

Regardless, if it is a difference between fortran and julia, then it is debuggable. Have you tried checking out #393 and running the three different integrators and looking at W then? galerkin doesn't have a W matrix implemented yet but it is not too complicated and there is no time like the present :P