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

未关闭
  • 哪里不对？
#401 2 条评论 0 个 reaction 已指派 0 人 在 GitHub 查看

维护者通常 3 天内回复

还没有人认领这个 Issue。

评估

难度
4/5
预计耗时
3-5 天
新手友好度
42/100
Issue 类型
缺陷
描述清晰度
描述清楚
活跃度
冷清
技术栈
julia
领域
backend

调研方向

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).

由索引模型根据 Issue 内容生成。

描述

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.

主要语言
Julia
星标
11
派生
7
平均合并
4 天 1 小时
30 天内合并 PR
68

环境准备

  • 没有 Dockerfile 或 Docker Compose 文件
  • 有 Pull Request 模板
  • 没有贡献指南

从这里开始

  1. 先读完整个 Issue，再读项目的贡献指南。
  2. 在 Issue 下留言说明你要接手 —— 这能避免两个人做同样的事。
  3. Fork 仓库，在一个分支上完成修改。
  4. 提交 Pull Request，并在描述里引用这个 Issue 编号。

OpenFUSIONToolkit/GPEC 的其他 Issue

查看 OpenFUSIONToolkit/GPEC 的全部 Issue

相似的 Issue

更多 Julia Issue

把新 issue 发到你的邮箱

精选适合新手参与的 GitHub issue 摘要。