PE - BUG - Im(py) is ~900x larger than Fortran GPEC, producing a spurious toroidal torque

Open
  • Is something wrong?
#401 2 comments 0 reactions 0 assignees View on GitHub

Maintainers usually reply within 3 days

Nobody has claimed this yet.

Assessment

Difficulty
4/5
Estimated time
3-5 days
Newbie friendliness
42/100
Issue type
Bug
Clarity
Clearly specified
Activity status
Quiet
Tech stack
julia
Domain
backend

Research direction

Compare Fortran gpec/gpresp.f (explicit Hermitian symmetrize of surf_indmats around line 284) with Julia calc_surface_inductance in ResponseMatrices.jl (~306) and how surface_inductance feeds permeability/reluctance and response_flux in Response.jl. Reproduce with examples/DIIID-like_ideal_example and docs/examples/run_ideal_example; done when Im(py)/toroidal_torque matches Fortran (~1e-5 relative) on the ideal deck without changing Re(py).

Written by the indexing model from the issue text.

Description

Split out from #370, where the symptom (a toroidal torque that should be zero) was traced to its
source.

The measurement

toroidal_torque = -2·n·Im(py) where py = ⟨Φ_resp, Λ⁻¹·Φ_resp⟩/4 — the imaginary part of a
quadratic form, identically zero for Hermitian Λ. Julia gives ~0.05 on the DIII-D-like ideal deck.

Fortran GPEC does not reproduce it. Rebuilt from the exact checkout (v1.5.7-611-ge68d7ac2,
gfortran 15.2), running docs/examples/run_ideal_example with inputs unmodified. The decks match:
byte-identical TkMkr_D3Dlike_Hmode.geqdsk, same d3d C-coil currents
[1000, 500, −500, −1000, −500, 500], hamada, nn=1, singfac_min=1e-4, and resp_index=0 —
the energy-identity method, which is the one Julia implements.

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 chain feeding it
(equilibrium, DCON solve, coil forcing, inductances, response flux) is faithfully ported. The
imaginary part is off by ~900×. An older Fortran build (v1.5.5-323, ~290 commits back) gives
8.1e-5, so the Fortran value is stable and ~1e-5 relative — zero to numerical precision, as the
Hermitian argument requires.

What is already ruled out

Not a discretization error. |T| is refinement-insensitive across five axes: mpsi (128→2048),
mtheta (256→512), mthvac (512→2048), eulerlagrange_tolerance (1e-10→1e-12), and
singfac_min (1e-3→1e-6). The last of these converges — first-order, each decade cutting the
change ~11× — onto a finite nonzero 5.0758e-2 rather than vanishing. Details in #370.

Not a response-method mismatch: Fortran runs resp_index=0, the method Julia ports.

Not the plasma-inductance assembly formula. Fortran gpresp.f:201-213 and Julia
ResponseMatrices.jl:256-264 are algebraically identical — Fortran's
t1 = -1/(chi1·s_i·2π·i) equals Julia's t1 = im/(chi1·s[i]·2π), likewise t2, and both form
2·t1·wt0[i,j]·t2 and invert it.

Lead: the surface inductance is not symmetrized

Fortran symmetrizes it explicitly before use (gpec/gpresp.f:284):

temp1=0.5*(temp1+CONJG(TRANSPOSE(temp1)))
surf_indmats=temp1

Julia's calc_surface_inductance has no equivalent step — it returns
inv(I_v) .* (μ₀ * (2π)^2) (ResponseMatrices.jl:306).

py is built from Λ (plasma), not the surface inductance directly, but response_flux is obtained
through the permeability and reluctance, both of which are built from surface_inductance
(Response.jl:42,47). An anti-Hermitian component surviving there propagates into Φ_resp and so
into Im(py), while leaving Re(py) essentially untouched — which is the signature observed.

Worth checking alongside: Response.jl:45 asserts "Λ (plasma inductance) is not Hermitian — its
anti-Hermitian part is the dissipative/torque response". For an ideal run with no dissipation that
part should vanish, so it is worth confirming that comment describes the intended physics here
rather than licensing the residual.

I have not verified that symmetrizing fixes the number — flagging it as the most specific divergence
found, not a diagnosed cause.

Reproducing

Julia: examples/DIIID-like_ideal_example as shipped; read
PerturbedEquilibrium/Energies/{toroidal_torque, plasma_energy} from gpec.h5.

Fortran: docs/examples/run_ideal_example with dcon then gpec; the values are printed as
"Required energy/torque to perturb plasma" and written as toroidal_torque (gpout.f:1353).

Consequence for regression testing

The torque is pinned in the golden-value work (#397) as a bound, |T| < atol with target 0, rather
than as a value — pinning a number that should be zero would encode this as truth. Fortran's ~1e-5
relative is the scale that bound should eventually reflect.

Dominant language
Julia
Stars
11
Forks
7
Avg merge
4d 1h
Merged PRs (30d)
68

Getting set up

  • No Dockerfile or Docker Compose file
  • Has a pull request template
  • No contributing guide

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

More from OpenFUSIONToolkit/GPEC

All issues in OpenFUSIONToolkit/GPEC

Similar issues

More Julia issues

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.