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.
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)wherepy = ⟨Φ_resp, Λ⁻¹·Φ_resp⟩/4— the imaginary part of aquadratic 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_examplewith inputs unmodified. The decks match:byte-identical
TkMkr_D3Dlike_Hmode.geqdsk, samed3dC-coil currents[1000, 500, −500, −1000, −500, 500],hamada,nn=1,singfac_min=1e-4, andresp_index=0—the energy-identity method, which is the one Julia implements.
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 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) gives8.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), andsingfac_min(1e-3→1e-6). The last of these converges — first-order, each decade cutting thechange ~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-213and JuliaResponseMatrices.jl:256-264are algebraically identical — Fortran'st1 = -1/(chi1·s_i·2π·i)equals Julia'st1 = im/(chi1·s[i]·2π), likewiset2, and both form2·t1·wt0[i,j]·t2and invert it.Lead: the surface inductance is not symmetrized
Fortran symmetrizes it explicitly before use (
gpec/gpresp.f:284):Julia's
calc_surface_inductancehas no equivalent step — it returnsinv(I_v) .* (μ₀ * (2π)^2)(ResponseMatrices.jl:306).pyis built from Λ (plasma), not the surface inductance directly, butresponse_fluxis obtainedthrough the permeability and reluctance, both of which are built from
surface_inductance(
Response.jl:42,47). An anti-Hermitian component surviving there propagates intoΦ_respand sointo
Im(py), while leavingRe(py)essentially untouched — which is the signature observed.Worth checking alongside:
Response.jl:45asserts "Λ (plasma inductance) is not Hermitian — itsanti-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_exampleas shipped; readPerturbedEquilibrium/Energies/{toroidal_torque, plasma_energy}fromgpec.h5.Fortran:
docs/examples/run_ideal_examplewithdconthengpec; 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| < atolwith target 0, ratherthan 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.