KineticForces Module (formerly PENTRC)

[logan-nc]

KineticForces Module (formerly PENTRC)

Complete refactor of the Fortran PENTRC module into a Julia KineticForces module, fully integrated into the GPEC architecture.

Refactoring (Commits 1-9)

• Renamed PENTRC → KineticForces throughout
• Converted module-level globals to struct-based parameter passing (EnergyParams, PitchParams)
• Standardized IO: accumulate results in KineticForcesState, write to gpec.h5
• Renamed LSODE functions to Julia conventions
• Implemented pitch angle ODE integration
• Cleaned up docs, removed all PENTRC references

API Streamlining (Commits 10-11)

• Deleted Input.jl (file IO stubs no longer needed — data passed via structs)
• Deleted Grid.jl (only dynamic ODE integration supported)
• Inlined get_method_flags, added KineticForcesInternal(equil) constructor
• Fixed b0/bo naming bug in equilibrium initialization

GAR NTV Implementation (Commit 12)

• New BounceAveraging.jl: powspace() grid generation, compute_bounce_data() for bounce-averaged ωb(λ), ωd(λ), |δJ|²(λ), and optional W vector outer products for kinetic matrix assembly
• Full calculate_gar(): Replaces stub with complete implementation — bounce data → fbnce interpolant → median normalization → pitch ODE → torque normalization (Eq. 19, Logan et al. 2013)
• Dynamic ODE ψ integration: integrate_psi_ode() / psi_integrand!() replaces trapezoidal rule on fixed grid, matching Fortran tintgrl_lsode
• Multi-n support: Block-diagonal numpert_total × numpert_total × 6 kinetic matrix assembly
• All three integration levels (ψ, λ, energy) use adaptive ODE solvers
• Supports fgar/tgar/pgar scalar torque and wmm/tmm/kmm/rmm matrix paths

Remaining TODOs

• build_kinetic_metric_matrices! in Fourfit.jl (s/t/x/y/z geometric matrices for matrix path)
• compute_kinetic_contribution() body (contracts matrices with ξ displacements)
• Populating dbob_m/divx_m/xs_m interpolants from PerturbedEquilibriumState

🤖 Generated with Claude Code

[logan-nc]

@mfairborn23 Please add whatever notes you took in our last meeting to the above summary (you can edit my post)

It is nice to have the PR open so we can al track the files changed and exchange questions / suggestions using the in-line comments. It is clear it is a WIP - so no worries there.

Let's open a separate branch for the kinetic modifications in fourfit.jl that put the kinetic terms into the stability calculation. I don't want this one to depend on that one or that one to have to wait if this one takes some more time. That one can use dummy kinetic matrices for now until we get this one merged.

[logan-nc]

@claude review this PR

[claude]

Code Review: PR #112 — KineticForces Module

This is a substantial PR that introduces the new KineticForces module. The architecture is sound — the callback injection pattern to avoid inverting ForceFreeStates→KineticForces load order is clean, and the struct-based parameter passing replacing Fortran globals is a good improvement. The benchmark scripts are thorough and the regression case infrastructure is welcome. Below are issues worth addressing before merge.

---

Bugs

Hardcoded absolute paths in examples/a10_kinetic_example/gpec.toml

eq_filename = "/Users/nlogan/Code/gpec/docs/examples/a10_ideal_example/fix_a100_k10_q2_bn010_prof1"
kinetic_file = "/Users/nlogan/Code/gpec/docs/examples/a10_ideal_example/a10_prof1.txt"

These paths will fail for everyone else. The committed example file is non-functional. Either document it as a template that requires local editing, or make the paths relative / use a placeholder.

bmax != B_vals[ibmax] check is dead code (Torque.jl, tpsi!)

bmax = maximum(B_vals)
ibmax = argmax(B_vals)
if bmax != B_vals[ibmax]
    error(...)
end

B_vals[ibmax] == bmax by construction — argmax returns the index of maximum. This error can never trigger; the check should be removed.

build_kinetic_metric_matrices is single-n only without adequate warning

At Fourfit.jl:3508 (in the diff), the y/z matrix loop uses n = intr.nlow as a fixed constant. For multi-n runs (npert > 1) this silently gives wrong matrices for n ≠ nlow. The function docstring doesn't flag this as a known limitation. At minimum, add an assertion or @warn when intr.npert > 1.

---

Performance

Any-typed fields on hot-path structs prevent type inference

In KineticForcesInternal, the fields dbob_m, divx_m, smats, tmats, xmats, ymats, zmats, xs_m are all typed ::Any. Similarly MethodResult.torque_profile is ::Any. These fields are evaluated inside the torque integration loops. Julia cannot infer return types through Any field access, so every interpolant call incurs dynamic dispatch. These should be typed (e.g., as a union of Nothing and the concrete interpolant type, or by making the struct parametric).

tspl spline rebuilt on every tpsi! call

A new CubicSeriesInterpolant tspl is constructed from scratch at each call to tpsi! (Torque.jl), which is invoked O(hundreds) of times per ψ quadrature. The pre-allocated buffer arrays (tpsi_B, tpsi_jac, etc.) on KineticForcesInternal suggest the intent to amortize setup cost, but the cubic_interp(xs, Series(hcat(...))) call at the end still allocates a new spline object each time. Pre-computing and caching the poloidal spline per ψ surface (e.g., via a small ψ → tspl dict) or restructuring to build it once per integrate_psi_ode call would avoid repeated construction.

expm allocation inside calculate_fcgl inner loop

expm = exp(im * twopi * intr.mfac * theta)

This allocates a new vector every θ iteration. Pre-allocate once and use @. expm = exp(im * twopi * intr.mfac * theta).

---

Code Quality

calculate_rlar is an incomplete placeholder but lacks a warning

calculate_rlar has a full docstring and no warning, but uses:

kappaint_val = sqrt(mean(abs.(dbob_m_f).^2))  # Placeholder: simplified estimate

This is not a physical implementation. calculate_clar correctly emits @warn "CLAR method not yet fully implemented, returning zero". calculate_rlar should do the same, or its docstring should state explicitly that it is a stub.

Bare catch swallows unexpected exceptions in find_kinetic_singular_surfaces!

try
    cond_vals[i] = evaluate_fbar_condition(...)
catch
    cond_vals[i] = Inf
end

A bare catch will silently absorb any error (e.g., a MethodError or BoundsError in the interpolants) and substitute Inf, causing silent incorrect behavior. At minimum catch e and rethrow unless the exception is LinearAlgebra.SingularException or DomainError.

Typo in citations.md

"providing the theoretical foundation for the the KineticForces functionality."

---

Minor

• The FourceFreeStates rename dmats → dmats_prim and emats → emats_prim is a breaking change to the FourFitVars public API; the old names no longer exist. Downstream code (or any user calling ffit.dmats directly) will get a FieldError. If these are considered internal fields, no action needed — just noting it for release notes.

• The new regression case solovev_kinetic_calculated adds test data but the test itself is commented out (the kinetic validation is explicitly deferred). This is useful infrastructure, but the regression harness will report this case as having no baseline — worth a note in the PR that this case is seeded for future use.

[logan-nc]

@matt-pharr you don't necessarily have to review all this - I just had to put someone. The merge of #224 into here actually contains all the documentation, plots and such.

I have claude working through it's own review points now