Equilibrium - BUGFIX - Converge PEST Δ′ in the radial grid

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

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

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

評価

難易度
5/5
見積もり時間
1週間以上
初心者へのやさしさ
6/100
issue の種類
バグ
明瞭さ
おおむね明確
活発さ
活発
技術スタック
julia
領域
compilers

調査の方向性

The defect sits in the Jacobian formula in DirectEquilibrium.jl (equilibrium_solver), with the same check needed in arclength_fieldline_int in DirectEquilibriumArcLength.jl. Start by reproducing the mpsi 256 vs 1024 table for PEST and hamada on the Solovev case. Done means PEST and hamada Δ′ agree to 1e-3 relative and are flat in mpsi, with the harness re-pinned.

索引モデルが issue の本文から書いたものです。

説明

Problem

On a direct (traced) equilibrium, Δ′ in PEST coordinates does not converge as mpsi grows. The same equilibrium in hamada coordinates does. Δ′ at Mercier index α is sensitive to coefficient error at the surface, amplified by roughly x^(−2α), so a small geometry error that is invisible in the energies shows up in Δ′.

Solovev, q0 = 1.2 (one rational at q = 2, α = 0.886), wall a = 0.05, etol 1e-10, eulerlagrange_tolerance 1e-12, Riccati, ldp grid. Run on #398's head (81006c3a6):

mpsi hamada PEST PEST + direct J
256 −148.86 −151.76 −151.90
1024 −149.81 −129.21 −135.20
  • Hamada moves 0.6% from 256 to 1024. PEST moves 15%, in the wrong direction.
  • et₁ agrees across all six runs to 2e-6, so this is a Δ′-only defect.
  • The exact-coefficient reference for this case is −150.73 (local analytic patch, settings differ slightly).

Likely sources

  1. J from a spline θ-derivative. PEST forms the Jacobian as (1 + ff_deriv_val[4]) * y_out[end, 2] * 2π * psio (DirectEquilibrium.jl, equilibrium_solver), i.e. by differentiating a periodic cubic spline of ∫dl/B_p. Hamada does not need that derivative. Before #398, ⟨∂J/∂ψ⟩ from this formula was off by 88% rms against the exact Solovev value.
  2. θ resolution of the traced surface. After #398 each surface is sampled only at the mtheta + 1 straight-fieldline nodes. A prototype that fixed this case to ≤1e-4 across mpsi 128–2048 in both PEST and hamada used 2048 field-line points per surface plus direct J. On #398, direct J alone recovers about a third of the gap at mpsi 1024.

Suggested approach

  • Compute J directly at each θ node from the field there: J = 2π ψ₀ ∮(w dl/B_p) / w(θ), where w is the Jacobian weight (power_bp, power_b, power_r, power_rc). Snippet below; it is a drop-in replacement for the one line above.
  • Then find the remaining error: rerun the table with mtheta doubled and quadrupled, and with a tighter dense-output tolerance on the field-line solve. Whichever closes the PEST/hamada gap at mpsi 1024 is the second half of the fix.
  • Check arclength_fieldline_int (DirectEquilibriumArcLength.jl) gets the same J.
rfac_j = sqrt(max(ff_val[1], 0.0)); eta_j = 2π * (theta + ff_val[2])
r_j = ro + rfac_j * cos(eta_j); z_j = zo + rfac_j * sin(eta_j)
direct_get_bfield!(bf_j, r_j, z_j, raw_profile.psi_in, raw_profile.sq_in, sqd_j, psio; derivs=1)
bp_j = sqrt(bf_j.br^2 + bf_j.bz^2); b_j = sqrt(bp_j^2 + (bf_j.f / r_j)^2)
w_j = bp_j^equil_params.power_bp * b_j^equil_params.power_b / (r_j^equil_params.power_r * rfac_j^equil_params.power_rc)
rzphi_nodes[ipsi, itheta, 4] = y_out[end, 5] / w_j * 2π * psio

(bf_j = DirectBField() and sqd_j = deriv1(raw_profile.sq_in) are allocated once before the ψ loop.)

Acceptance

  • Solovev α 0.89: PEST and hamada Δ′ agree to ≤1e-3 relative at mpsi 256 and 1024, and each is flat in mpsi to the same level.
  • DIII-D-like n = 1 Riccati: PEST and hamada Δ′ agree on all five surfaces once #506's profiles are in.
  • Every direct-equilibrium PEST result moves, so the PR takes ! and re-pins the harness.

Dependencies

  • Stacks on #398 (common straight-fieldline sampling; approved, needs its conflicts resolved). #398 removes the sparse-sampling noise that dominated before; this issue is what is left.
  • Independent of #506 (#504), #517 and the crossing work (#481).
主要言語
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 を短くまとめたダイジェスト。