PerturbedEquilibrium - FEATURE - Write the pmodb quantities (Eulerian and Lagrangian δ|B|, B·∇·ξ⊥, B·κ·ξ⊥) to gpec.h5

Aperta
  • Qualcosa non va?
#429 0 commenti 0 reazioni 0 assegnatari Vedi su GitHub

I maintainer di solito rispondono entro 3 giorni

Nessuno ha ancora preso questa issue.

Valutazione

Difficoltà
5/5
Tempo stimato
Più di una settimana
Idoneità per principianti
25/100
Tipo di issue
Funzionalità
Chiarezza
Abbastanza chiara
Stato di attività
Tranquilla
Stack tecnologico
julia
Ambito
data, data-visualization

Direzione di ricerca

Start by reading src/Analysis/PerturbedEquilibrium.jl and src/Analysis/PerturbedEquilibriumModes.jl, especially plot_mode_spectrogram and modes_to_theta; inspect the listed response and equilibrium datasets in gpec.h5. Decide which δB/B route to use, then implement the on-demand helper and companion plot without adding datasets. Done means the acceptance criteria are met and a matched Fortran benchmark reports amplitude and phase agreement, tolerance, and figures.

Scritto dal modello di indicizzazione a partire dal testo della issue.

Descrizione

feature

Summary

Julia GPEC has no equivalent of the Fortran pmodb_flag outputs, and it should not
get one by adding datasets to gpec.h5. Instead, add an Analysis helper that
reconstructs these quantities on demand from the response fields the file already
carries. This keeps gpec.h5 small while still letting users plot δ|B| and
∇·ξ⊥ profiles whenever they want them.

What Fortran provides

gpout_pmodb (gpec/gpout.f:2745) writes GPEC_pmod_n#.out plus these netCDF
variables, on both the (ψ, m) mode grid and the (ψ, θ) function grid:

Fortran variable Quantity
b_eul / b_eul_fun Eulerian δ|B| = χ₁(b_θ + q·b_ζ)/(J·B) [T]
b_lag / b_lag_fun Lagrangian δ|B| = b_eul + ξ^ψ ∂B/∂ψ + (ξ^α/(χ₁q)) ∂B/∂θ [T]
Bdivxi_perp / _fun B·∇·ξ⊥
Bkappaxi_perp / _fun B·κ·ξ⊥

The Fortran docs describe this file as "necessary for NTV post processing" — its
original job was to hand δ|B| from GPEC to PENTRC across a file boundary.

Current state in Julia

The physics exists, but only as in-memory interpolants inside KineticForces:

  • set_perturbation_data! (src/KineticForces/KineticForcesStructs.jl:293) builds
    kf_intr.dbob_m (δB/B) and kf_intr.divx_m (∇·ξ⊥) as CubicSeriesInterpolants
    over ψ, [npsi, mpert].
  • This is a port of PENTRC's set_peq (pentrc/inputs.f90:826-877) — the JBB
    deweighting route (apply S,T,X,Y,Z in m-space → IDFT → divide by J·B² → DFT),
    see _jbb_deweight! at line 412 — not a port of gpout_pmodb, which builds
    δ|B| directly from the covariant perturbed field.
  • Neither is written to HDF5. src/KineticForces/Output.jl writes only
    kinetic_forces/<method>/ (torques, energies, dTdψ, matrices, optional
    per-record bounce diagnostics).
  • Because KineticForces runs in-process, the file hand-off Fortran needed is
    obsolete. The Eulerian/Lagrangian split is never formed anywhere in Julia.

So a user who wants to look at δ|B|(ψ, m) today has no supported path.

Proposal

Add an Analysis helper — e.g. compute_pmodb(h5path; ...) in
src/Analysis/PerturbedEquilibrium.jl, with a companion plot function following the
existing plot_mode_spectrogram (src/Analysis/PerturbedEquilibrium.jl:353)
conventions — that reconstructs the gpout_pmodb quantities from datasets already
present in gpec.h5.

Everything required is already in the file:

  • perturbed_equilibrium/response/b_cova_theta, b_cova_zeta → b_θ, b_ζ
  • perturbed_equilibrium/response/clebsch_psi → ξ^ψ
  • perturbed_equilibrium/response/clebsch_alpha → ξ^α/χ₁ (already divided by χ₁)
  • integration/psi, integration/q → radial axis and q
  • splines/rzphi/* and the equilibrium field data → J, B, ∂B/∂ψ, ∂B/∂θ
  • info/mlow, info/mpert, info/nlow, info/npert → mode metadata

Analysis.PerturbedEquilibriumModes.modes_to_theta
(src/Analysis/PerturbedEquilibriumModes.jl:41) already handles the m → θ
transform, so the helper should reuse it rather than reimplementing the DFT.

Explicit non-goal: do not add b_eul/b_lag/divxi/kappaxi datasets to
gpec.h5. These are [npsi, mpert] (and potentially [npsi, mtheta]) complex
arrays per quantity per n, which is real disk cost for a diagnostic most runs never
look at. Everything needed to regenerate them is already stored; recomputing on
demand is cheap.

Open question for the implementer

Which route to use for δB/B:

  1. gpout_pmodb route — from the covariant field; gives both Eulerian and
    Lagrangian, and is exactly reproducible against the Fortran output.
  2. set_peq route — reuse the existing _jbb_deweight! machinery that
    KineticForces already runs, giving the single δB/B the action integral consumes.

Route 1 is the better fit for a user-facing plotting helper (it is what the Fortran
diagnostic actually emits, and it separates Eulerian from Lagrangian). Route 2 has
the advantage of exercising code already on the NTV hot path. Cross-checking the two
against each other would itself be a useful result and would surface any discrepancy
between the two formulations.

Benchmarking requirement

Any PR addressing this issue must include a benchmark against Fortran GPEC.
The helper is only useful if its numbers match, and there is a nontrivial risk of
convention mismatch: the Julia response fields live in the working jac_type
coordinate, whereas Fortran applies gpeq_bcoordsout with jac_out/tmag_out
before writing, and the θ-space *_fun outputs additionally carry the
exp(i·n·Δφ) factor to machine toroidal angle.

Concretely, a PR should:

  • Run a matched case in both codes (e.g. examples/DIIID-like_ideal_example against
    the corresponding Fortran gpec.in with pmodb_flag=t), with jac_out set so the
    coordinates genuinely correspond.
  • Compare b_eul, b_lag, Bdivxi_perp, and Bkappaxi_perp in both mode space
    and θ space, over the full ψ range — not just at the edge surface.
  • Report amplitude and phase agreement per m, and state the achieved tolerance.
  • Include the comparison figure(s) in the PR description, following the
    benchmarks/ conventions (inputs referenced from examples/, outputs written into
    benchmarks/, nothing committed).
  • Explicitly document any remaining discrepancy and its cause, rather than tuning a
    tolerance until it passes.

Note that the existing cylindrical-component helpers (xi_R, xi_Z, xi_phi,
b_R, b_Z, b_phi) are documented as beta with up to ~20% disagreement against
Fortran; this helper should not ship in that state.

Acceptance criteria

  • Analysis helper returns Eulerian and Lagrangian δ|B|, B·∇·ξ⊥, and B·κ·ξ⊥ in
    mode space, with an option for θ space.
  • Companion plot function following existing Analysis plotting conventions.
  • No new datasets written to gpec.h5.
  • Docstrings added and covered by an @autodocs block in docs/src/analysis.md.
  • Fortran benchmark performed, tolerance stated, figures in the PR description.
  • Regression harness run per project policy (this should be output-neutral —
    confirm it is).
Lingua principale
Julia
Stelle
11
Fork
7
Merge medio
4g 3h
PR unite (30g)
70

Preparare l'ambiente

  • Nessun Dockerfile né file Docker Compose
  • Ha un modello di pull request
  • Nessuna guida per i contributori

Come iniziare

  1. Leggi tutta la issue e poi la guida ai contributi del progetto.
  2. Commenta sulla issue per dire che te ne occupi tu — evita che due persone facciano lo stesso lavoro.
  3. Fai un fork del repository e lavora su un branch.
  4. Apri una pull request che faccia riferimento al numero della issue.

Altre issue di OpenFUSIONToolkit/GPEC

Tutte le issue di OpenFUSIONToolkit/GPEC

Issue simili

Altre issue su Julia

Ricevi le nuove issue nella tua casella

Un breve riepilogo di issue GitHub adatte ai principianti.