Repository navigation
Repo - GPEC 2.0 Fortran parity plan #591
Description
Survey of the 80 PRs merged into PrincetonUniversity/GPEC between 2026-04-25 and 2026-10-24, extended back to the start of the Julia port (2025-07-01, one week after the Fortran v1.5.6 tag; only 18 further Fortran PRs merged in between — section E). Every "lacks" claim was re-checked by hand against the Julia tree. Scope: physics correctness, not Fortran-specific performance or memory hygiene. Fortran PR numbers refer to PrincetonUniversity/GPEC.
Assigned to myself; this is the 2.0 parity work I am taking on. @matt-pharr — tagging you because several items build directly on your open API stack (A4 and A6 on the KineticForces/ForceFreeStates paths #569/#574 and #520 touch, A9 on the [Equilibrium] kinetic_file home #581 creates, A13 on the SingularCoupling reorganisation in #574), so the landing order of that stack decides where these branch from. Nothing here changes your PRs; it only sequences after them.
A. Physics items Julia must act on (priority order)
| # | Fortran PR | What Fortran corrected | Julia location | Status | Action | Size |
|---|---|---|---|---|---|---|
| A1 | PrincetonUniversity/GPEC#357 (10-24) | Conjugate the bounce-averaged W vectors before forming W_X/W_Y/W_Z; wrong kinetic EL matrices for up-down-asymmetric equilibria | KineticForces/BounceAveraging.jl:_bounce_harmonic_tail (wmu_ba[mi] = conj(dot(scr.int_w, wsamp))) |
done — Julia commit 92e2c95e0 (matt-pharr), merged 2026-10-25 in #472; Fortran 357 is the port of that finding. An earlier survey draft marked this as missing because it checked a tree 9 commits behind develop |
none; when the KF stack (#579–584) merges develop, confirm both conj(dot(…)) lines survive in the restructured _bounce_harmonic_tail |
— |
| A2 | PrincetonUniversity/GPEC#342 (10-22) | CGL curvature drive κ·ξ⊥ = -(δB/B + ∇·ξ⊥) (was -½·…) and the method aborted on a vector exp |
KineticForces/Torque.jl:calculate_fcgl |
in flight: PR #588 | review + merge | done |
| A3 | PrincetonUniversity/GPEC#320 (10-19) | nutype="small" Krook width `1e-5· |
ω_E | ` (sign flipped where precession reverses) | KineticForces/EnergyIntegration.jl:66 (1e-5 * p.we) |
lacks |
| A4 | PrincetonUniversity/GPEC#288 (08-02) | v∥ ≤ 0 fill of the bounce integrals copied the ω_b column into the ω_D column (fs(i-2,1) → fs(i-2,2)) |
KineticForces/BounceAveraging.jl:1049-1050 (fill!(view(g_wd,…), g_wb[i-1]), comment says "reproducing the Fortran behavior for parity") |
lacks, deliberately | trace the index convention (Julia i-1 ↔ Fortran i-2), fix to g_wd[i-1], re-baseline the KF harness rows (ω_D near the trapped/passing boundary moves). KineticForces - BUGFIX! |
S |
| A5 | PrincetonUniversity/GPEC#336 + PrincetonUniversity/GPEC#300 (10-19/20) | Kinetic-profile edge: floor T at 1e-6 keV and lnΛ at 1; monotone (pchip) slopes for ν_i, ν_e, lnΛ and the Zeff table (2-point tables allowed) | Equilibrium/KineticProfiles.jl:330-353 (loglam = … : 0.0, _cubic_resample with ExtendExtrap) |
lacks | FastInterpolations ships pchip; switch the tabulated ν/lnΛ/Zeff resamples to it, add the floors, test on a profile with a cold edge node. Equilibrium - BUGFIX! |
S |
| A6 | PrincetonUniversity/GPEC#332 (10-20) | EL ODE: a norm-ratio reduction set new=.TRUE., and the next singular-surface call reset the baseline instead of doing the mandatory singular reduction; fix fix_required = sing_flag ∨ uratio>ucrit independent of the reset |
ForceFreeStates/EulerLagrange.jl:compute_solution_norms! (sing_flag tested only in the else branch of if odet.new) |
lacks (likely) | confirm the singular call can arrive with odet.new == true (instrument one run); fix; expect ODE step-pattern rows to move. ForceFreeStates - BUGFIX! |
S |
| A7 | PrincetonUniversity/GPEC#305 + PrincetonUniversity/GPEC#301 (10-15) | Trace-aware Frobenius exponents p = tr/2 ∓ √(tr²/4 − det) (was ∓√(−det)); Galerkin Cholesky failure falls back to LU |
ForceFreeStates/Surfaces/Asymptotics.jl:69 still sqrt(-det) on develop |
fixed on an unmerged branch: Julia #517 merged into feature/galerkin-cholesky = PR #513 (draft, CHANGES_REQUESTED, 15 commits ahead of develop) |
land #513 (it is the Julia origin of Fortran 305/301); nothing to port | review |
| A8 | PrincetonUniversity/GPEC#303 (10-15) | Galerkin F built from the Cholesky-factor spline F = (QL)(QL)^H so the Galerkin matrix and the asymptotics see the same F |
ForceFreeStates/Galerkin/GalerkinAssembly.jl:78 uses the un-factored F_spline_gal |
verify | check whether Julia's asymptotics use a factored F at all (Fourfit.jl builds both); if they differ, align as Fortran did; may be covered by #513 — check its diff first | S–M |
| A9 | PrincetonUniversity/GPEC#345 (10-20) | kinetic_format="rates-v1": supplied lnΛ, ν_i, ν_e columns replace the internal collision model (log-interpolated) |
Equilibrium/KineticProfiles.jl:326 (always recomputes) |
lacks (input capability) | add the 9-column format and a kinetic_format key; reuse the pchip path of A5. Equilibrium - FEATURE |
M |
| A10 | PrincetonUniversity/GPEC#297 (10-15) | profile_source="integrate": build F and p in the profile spline from FF′ and p′ when the file's F/p are inconsistent |
Equilibrium/ReadEquilibrium.jl:86-87 reads ffprime/pprime but never uses them |
lacks | port as an equilibrium input option with a consistency diagnostic. Equilibrium - FEATURE |
M |
| A11 | PrincetonUniversity/GPEC#289 + PrincetonUniversity/GPEC#349 (08-02, 10-21) | CLAR normalisations (ω_b without bo under the root, ω_D over an extra bo); RLAR κ-integral buffer |
Torque.jl:calculate_rlar is a placeholder (RMS of δB/B), calculate_clar returns 0 with a warning |
not ported | 2.0 decision: port RLAR/CLAR faithfully (use the post-289 forms) or drop the flags from the registry so a deck cannot silently get a stub. Recommend drop unless a user needs them | S (drop) / L (port) |
| A12 | PrincetonUniversity/GPEC#323 (10-20) | Guard the rational-surface list when it is empty or exhausted | ForceFreeStates/Surfaces/Finding.jl (no explicit isempty guard found) |
verify | run a q range with no rational surface; add a guard + test if it errors | XS |
| A13 | PrincetonUniversity/GPEC#330 + PrincetonUniversity/GPEC#185 (10-16/20) | Singular-coupling overlap export: W^H V instead of CONJG(W)·V; units and descriptions of the overlap outputs; dominant-mode overlap |
PerturbedEquilibrium/SingularCoupling.jl (no SVD overlap export found) |
verify / feature gap | decide whether 2.0 ships the overlap diagnostics; if yes, port with the corrected adjoint | M |
| A14 | PrincetonUniversity/GPEC#269 (06-12) | Positive-q enforcement; drop the newq0 == -1 → -q0 convention |
Equilibrium/InverseEquilibrium.jl:315-316 still flips |
lacks (low impact) | align with the documented sign convention; Equilibrium - MINOR |
XS |
| A15 | PrincetonUniversity/GPEC#308 / PrincetonUniversity/GPEC#317 (10-15/16) | Spline endpoint-derivative constant offsets; r8 Cardano exponent in spline roots | FastInterpolations (cubic/cubic_solver.jl boundary conditions); Julia's own cubic roots (_real_cubic_roots) already conditioned (#483) |
n/a (likely) | one check of FastInterpolations' PolyFit endpoint derivative against a linear + constant test function | XS |
B. Already in Julia (evidence in the fragments)
| Fortran PR | Item | Julia evidence |
|---|---|---|
| PrincetonUniversity/GPEC#281 | ω_D in the GAR branch carried a spurious major radius | BounceAveraging.jl:361 wdbar = ro * bo * … |
| PrincetonUniversity/GPEC#284 | periodic v∥ spline | Torque.jl B_vpar = cubic_interp(…; bc=PeriodicBC()) |
| PrincetonUniversity/GPEC#341 | pmodb drives: keep the supplied ∇·ξ⊥, δB/B = -(∇·ξ⊥ + κ·ξ⊥) |
KineticForcesStructs.jl:443 jbb_dbob = -(jbb_divx + jbb_kapx) |
| 280 (ladder half) | energy rtol = 1e-2 × pitch rtol | Julia #582 |
| PrincetonUniversity/GPEC#329 | sample direct equilibria at common straight-field-line angles | Julia #398 is the origin |
| PrincetonUniversity/GPEC#328 | qlim within qhigh | Finding.jl:125 qlim = min(qmax, qhigh) |
| PrincetonUniversity/GPEC#295 | vacuum matrix in the plasma Fourier frame | Julia removes the θ reversal (commits 11971072c, 29b4de02d) |
| PrincetonUniversity/GPEC#298 / PrincetonUniversity/GPEC#335 | ldp_i ψ-bounds record; grid_type="original" bounded by psilow/psihigh |
ReadEquilibrium.jl:413-456; InverseEquilibrium.jl:181-184 (direct path not inspected) |
| PrincetonUniversity/GPEC#304 | Taylor-slot typos in sing_mmat |
Asymptotics.jl:280, 373 indices look right (not compared line by line) |
| PrincetonUniversity/GPEC#326 / PrincetonUniversity/GPEC#309 / PrincetonUniversity/GPEC#310 | Galerkin element sizing n·q1 both sides; Jacobi-scaled LU; complex W_v at the edge |
GalerkinGrid.jl:101,133 (docstring stale); GalerkinSolve.jl:149; wv_edge::ComplexF64 |
| PrincetonUniversity/GPEC#173 | SLAYER growth-rate machinery | Julia is the origin (#463, #403, #434, #431, #433) |
| PrincetonUniversity/GPEC#296 | Zeff profile | zeff vector in kinetic profiles; SLAYER consumes it (no RDCON mufrac path) |
C. Not applicable to Julia, or feature parity rather than a fix
- LSODE-specific (312 stale DINTDY export, 314/334 ISTATE stops): Julia's ψ, λ and energy integrals are adaptive Gauss–Kronrod;
check_psi_quadrature_convergencewarns on non-convergence and the energy integrand has anisfinitecheck. Worth one audit that a failed quadrature cannot pass silently (XS). - Fortran memory / types / formats (299, 294, 256, 250, 325, 315, 352, 351, 347, 327, 302, 354, 356, 337, 346, 348, 285, 279, 264): no Julia analogue, or Julia's structure makes the bug impossible (356: R computed per θ locally; 354: no fixed-width buffers).
- pmodb diagnostic output (340, 343, 333): Julia has no pmodb file output; the drives it consumes are formed correctly (341). Feature gap only if 2.0 wants the diagnostic.
- SLAYER thresholds (265, 259, 274
Phi_res_crit_callen, 239 island-saturation pressure loss): no torque-balance threshold scan in Julia's InnerLayer; Julia has per-surface Prandtl numbers. Feature-parity list. - 213 Bernstein δW_p breakdown: no counterpart in this tree (the PR mentions a Julia port elsewhere — ask the author). 210 resistive singular coupling: partial (Tearing has the coupled dispersion; RMATCH-chain parity unverified).
- Docs (266, 268/269 sign conventions, 282 PENTRC rotation/harmonic sign statements, 290, 339 pmodb gauges):
docs/src/conventions.mdexists; a 2.0 docs pass should mirror 282 and 339. - Release/merge PRs (358, 359, 263, 260, 271, 245, 278).
E. Earlier window, 2025-07-01 → 2026-04-24 (18 PRs)
| Fortran PR | Item | Julia status |
|---|---|---|
| PrincetonUniversity/GPEC#249 (03-24) | SLAYER r_s, shear and b_l normalisation | has: InnerLayer/SLAYER/LayerParameters.jl uses the true shear rs·q1/(q·da_dpsi); one open check that _rs_at is ρ_tor and _da_dpsi_at is dρ_tor/dψ_n in LayerInputs.jl (the Fortran bug was passing the derivative as r_s) |
| PrincetonUniversity/GPEC#229 (08-27) | Vacuum Green's functions with Chance 2007 Legendre improvements | has: Vacuum/Kernel2D.jl green(…; uselegacygreenfunction=false) default, 32-point Gauss cache; nodes/weights not diffed against Fortran |
| PrincetonUniversity/GPEC#225 / PrincetonUniversity/GPEC#234 (09-19/23) | RDCON/STRIDE truncation at qlow/qhigh standardised | has (qlow handling mirrors RDCON sing_min); the qa = qlim change in free.f not traced |
| PrincetonUniversity/GPEC#248, PrincetonUniversity/GPEC#247, PrincetonUniversity/GPEC#227 | nn == 0 guard; LSODE step-cap control; eq_type names | has (Julia's gal_gnstep caps QuadGK evaluations rather than LSODE steps) |
| PrincetonUniversity/GPEC#237 (02-06) | Callen vacuum-island widths from the cubic A·w³ + B·w + C = 0 (hw_min, hw_sat, hw_v_crit) |
lacks — feature parity (goes with 239/274 in section C) |
| PrincetonUniversity/GPEC#224 + PrincetonUniversity/GPEC#242 (07-27, 03-10) | Modified Rutherford equation surface terms; fixed prefactors and the trapped-fraction form 1 − √((1−ε)/(1+ε))·(1−ε_eff)/(1+2√ε_eff), ε_eff = max(0, 0.67(1 − 1.4δ |
δ |
| PrincetonUniversity/GPEC#246, PrincetonUniversity/GPEC#219, PrincetonUniversity/GPEC#230, PrincetonUniversity/GPEC#254, PrincetonUniversity/GPEC#255, PrincetonUniversity/GPEC#238, PrincetonUniversity/GPEC#232, PrincetonUniversity/GPEC#231 | field-line step cap, uninitialised-variable guard, build, headers, docs | n/a |
Conclusion for the earlier window: no missing correctness fix; two feature-parity items (Callen island widths, MRE terms) join the section C list.
D. Proposed order of work for 2.0
- Land what exists: KineticForces - BUGFIX! - Make the fluid CGL energy run and use the full curvature drive #588 (A2), ForceFreeStates.Galerkin - FEATURE! - Solve the Galerkin system by banded Cholesky, with LU fallback #513 after its review items (A7, and check A8 inside it).
- KineticForces correctness: A3 is KineticForces - BUGFIX! - Keep the small collision width positive when precession reverses #561 (review it); A4 alone in one
KineticForces - BUGFIX!(tiny, moves KF rows near the trapped/passing boundary). A1 landed in PerturbedEquilibrium - FEATURE - Add the psi-resolved torque response matrices T_xe and T_source #472. - ForceFreeStates: A6 (instrument, fix, re-baseline the step-pattern rows — the harness noise thresholds flagged in KineticForces - PERF! - Integrate the energy resonance operator in u = √x #584 should be revisited at the same time); A12 guard.
- Equilibrium inputs: A5 (pchip + floors) then A9 (
rates-v1) and A10 (profile_source="integrate") as two FEATURE PRs sharing the pchip path. - Decisions for the maintainer: A11 (drop or port RLAR/CLAR), A13 (ship overlap diagnostics?), pmodb output, SLAYER thresholds, Bernstein breakdown — feature-parity scope for 2.0 rather than correctness.
- Docs pass mirroring 282/339/266.
Verification standard for every port: unit test pinning the corrected expression, harness run against develop with the moved rows explained, and where the Fortran PR changed physics (A4, A6) a direct comparison with Fortran 1.6.0 (Fortran v1.6.0, not the 1.5-era gpec/0.0 module) on the same deck.
F. Plan items already covered or re-based by open PRs
| item | open PR | verdict |
|---|---|---|
| A2 CGL | #588 (mine) duplicates krystophny's #541 + #542 (10-19, drafts, stacked; broadcast + the -(dbob+divx) curvature with a FortSym symbolic proof and a closure test). #588 is a superset (both fixes, closed-form unit test, example runs fcgl, harness row) |
consolidate: ask krystophny to review #588, credit 541/542 in its body, close them — or the reverse; maintainer's call |
| A3 Krook width | #561 is exactly the plan's fix (1e-5 * abs(p.we) + sign test); reviewer matt-pharr |
done once reviewed; link Fortran 320 in its body |
| A10 profile_source | #506 (d-burg, draft) adds profile_source = "derivatives" | "values" — but makes derivatives the default and moves every g-file/IMAS result (Δ′ diag 15 %, NTV 2.8 %, SLAYER γ 6 %). Fortran 297 added it as an opt-in |
review #506; decide the default consciously (the Fortran default is still the tabulated values) |
| A7 Frobenius | #513 | unchanged |
| A6 norm reset | #520 (base #489, approved) rewrites compute_solution_norms! and half-fixes it: an ideal crossing (resonant_rows given) now always reduces, but cross_kinetic_singular_surf! still calls with sing_flag=true and no rows, so odet.new && !at_ideal_crossing resets and skips the forced reduction — the Fortran 332 hole survives on the kinetic path. Fortran's fix is fix_required = sing_flag before the branch |
build A6 on #520: one-line if odet.new && !sing_flag, plus the Fortran IF(sing_flag) unorm(index(1:msol)) = … pivot-order line (1f8abc54) if Julia's apply_gaussian_reduction! has the same ordering dependence. Do not do it on develop; #520 also changes what triggers a reset |
| A12 empty surface list | #492 adds named errors and an msing = 0 fixture but no empty-list guard |
still verify |
| A1 (357 conj) | merged in #472 (92e2c95e0), 2026-10-25 |
nothing to do; the KF stack must carry the two conj lines through its merge of develop |
| A4 (288 ω_D fill) | BounceAveraging.jl is contended by #583 and #569 |
build on the KF stack tip once it is reconciled with #569; not on develop |
| A5 pchip + floors | KineticProfiles.jl resample internals are not restructured by #581 (it changes the entry points) |
can go on develop now; small |
A9 rates-v1 |
needs the [Equilibrium] kinetic_file home that #581 creates |
wait for #581 |
| A13 overlap export | #574 reorganises SingularCoupling/ResonantCoupling |
build on #574 |
| #543 (mine, draft) | #574 moves response data to pe.response |
rebase onto #574 before un-drafting |
Metadata
Metadata
Assignees
Labels
No labels
Activity