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

Abierto
  • ¿Algo no está bien?
#429 0 comentarios 0 reacciones 0 asignados Ver en GitHub

Los mantenedores suelen responder en 3 días

Nadie ha tomado este issue todavía.

Evaluación

Dificultad
5/5
Tiempo estimado
Más de una semana
Aptitud para principiantes
25/100
Tipo de issue
Nueva funcionalidad
Claridad
Bastante claro
Estado de actividad
Tranquilo
Stack tecnológico
julia
Área
data, data-visualization

Línea de trabajo

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.

Escrito por el modelo de indexación a partir del texto del issue.

Descripción

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).
Lenguaje dominante
Julia
Estrellas
11
Forks
7
Merge medio
4 d 3 h
PR fusionados (30 d)
70

Preparar el entorno

  • Sin Dockerfile ni archivo de Docker Compose
  • Tiene una plantilla de pull request
  • Sin guía de contribución

Primeros pasos

  1. Lee el issue completo y luego la guía de contribución del proyecto.
  2. Comenta en el issue que vas a ocuparte — evita que dos personas hagan lo mismo.
  3. Haz un fork del repositorio y trabaja en una rama.
  4. Abre un pull request que haga referencia al número del issue.

Más de OpenFUSIONToolkit/GPEC

Todos los issues de OpenFUSIONToolkit/GPEC

Issues similares

Más issues de Julia

Recibe los nuevos issues en tu correo

Un resumen breve de issues de GitHub para principiantes.