Repository navigation
PE toroidal torque changes sign under grid refinement while its magnitude stays stable to ~1% #370
Description
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
Activity
matt-pharr commentedon Aug 14, 2026
@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 commentedon Aug 16, 2026
@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:pe_torquempsi(ldp)mtheta@ mpsi=1024mthvac(auto grid)eulerlagrange_toleranceThe 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 resolutionwarning — 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 commentedon Aug 17, 2026
@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 commentedon Aug 17, 2026
@matt-pharr thanks — the data suggest the
caroute can't answer this.sing_get_casolvesua·ca = u, and sinceEulerLagrange.jl:639introduces the small resonant solution atua[:, ipert_res + numpert_total, :],ca[:,:,1]holds the large-solution coefficients — so the resonant row ofca_left/ca_rightisexactly 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
cablocks span 45–58 orders of magnitude, and against that roundofffloor two of the five surfaces sit below it and two more are within 3×:
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:533stops the sweep atpsi_end = sing.psifac - singfac_min/|n·q1|, and line 630 crosses with a fixed 2-step trapezoidalextrapolation
u .+= (du1 + du2)·dpsi.singfac_minis the only knob setting that jump width — andmpsi,mtheta,mthvacandeulerlagrange_toleranceall leave it untouched, which would explainwhy the torque ignored every one of them.
Scanning it over three decades (deck otherwise byte-identical):
singfac_minpe_torqueet[1]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 thecrossing extrapolation were manufacturing the torque, shrinking
dpsiby 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 aresolution 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
developnow shipsTkMkr_D3Dlike_Hmode.geqdskindocs/examples/run_kinetic_example/,byte-identical to our deck's (sha256
9449e54b8488…), so the comparison can besame-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 commentedon Aug 17, 2026
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 inputsunmodified. The decks match closely: byte-identical
TkMkr_D3Dlike_Hmode.geqdsk, samed3dC-coilcurrents
[1000, 500, −500, −1000, −500, 500],hamada,nn=1,singfac_min=1e-4.Both codes compute the same quadratic form — Julia
Response.jl:98is a direct port ofgpout.f:1353,py = ⟨Φ_resp, L_plasma⁻¹·Φ_resp⟩/4, with the energy asRe(py)and the torque as−2n·Im(py):Re(py)— energy to perturb plasma−2n·Im(py)— torqueThe 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 Fortranhistory 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 Fortrandoes not have. Worth noting the free-boundary energy carries a small imaginary part in Fortran too
(DCON reports
imaginary = 1.4e-04againstreal = 7.5e-01), so the mechanism exists there — it isjust ~1000× smaller.
Correction to my previous comment: I normalized the torque by
et[1], the DCON free-boundaryeigenvalue, giving 6.355%. The physically matched denominator is
Re(py)from the same quadraticform, which gives 1.48%. The
singfac_minconvergence result is unaffected — only the normalizationwas wrong.
Golden treatment unchanged, and now clearly right:
|T| < atolwith target 0, since Fortran showsthe 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 commentedon Aug 17, 2026
Split the root cause out to #401 —
Im(py)is ~900× larger than Fortran GPEC's, withRe(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 commentedon Aug 17, 2026
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