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

オープン
  • 何か問題がありますか？
#401 コメント 2 件 リアクション 0 件 担当者 0 名 GitHub で見る

メンテナーはふだん 3 日以内に返信

まだ誰も着手していません。

評価

難易度
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日 3時間
マージ済み PR（30日）
70

環境構築

  • Dockerfile・Docker Compose ファイルなし
  • プルリクエストのテンプレートあり
  • コントリビューションガイドなし

はじめの一歩

  1. issue を最後まで読み、次にプロジェクトのコントリビューションガイドを読みます。
  2. 着手することを issue にコメントします — 二人が同じ作業をするのを防げます。
  3. リポジトリをフォークし、ブランチを切って変更します。
  4. issue 番号を参照したプルリクエストを送ります。

OpenFUSIONToolkit/GPEC のほかの issue

OpenFUSIONToolkit/GPEC の issue をすべて見る

似ている issue

Julia の issue をもっと見る

新しい issue をメールで受け取る

初心者向けの GitHub issue を短くまとめたダイジェスト。