Skip to content

Equilibrium/KineticForces - BUGFIX! - Port the Fortran 1.6.0 kinetic profile, bounce-integral and torque-method corrections - #625

Merged
matt-pharr merged 17 commits into
developfrom
bugfix/fortran-parity-kin
Nov 2, 2026
Merged

matt-pharr merged 17 commits into
developfrom
bugfix/fortran-parity-kin

Conversation

@logan-nc

@logan-nc logan-nc commented Nov 1, 2026 •

Copy link
Copy Markdown
Collaborator

Release note

  • Audience: users
  • Numerical impact: kinetic results move (harness @ e0bb4f7): up to 9 % in the TJ calculated-kinetic eigenfunction norms and about 1 % in Solovev/TJ kinetic energies, from the kinetic corrections (the pchip resample of the kinetic tables, the ω_D hold and the collision floors). SLAYER growth rates move by up to 0.7 % (coupled γ), and the decks that build the two-pass ψ grid from kinetic profiles move by at most 1e-6. A per-commit run shows the later fixes (inverse newq0, bounce-point order, scaling-factor and grid checks) do not move the two largest movers. RLAR/CLAR and newq0 are not used by any harness deck.
  • Migration: a deck that sets rlar_flag or clar_flag stops with a pointer to fgar_flag / tgar_flag; newq0 < 0 (including the old -1 flip) is rejected.

These are the kinetic-profile and KineticForces corrections from Fortran GPEC 1.6.0, the Equilibrium and KineticForces half of #591:

  • the precession integrand past the bounce midpoint is held from itself, not copied from the bounce-frequency integrand;
  • a cold edge node gets a finite collisionality;
  • kinetic tables are fit with monotone pchip;
  • q must be positive;
  • the placeholder RLAR and CLAR torque methods are removed.

What changed, one commit per item

Commit #591 item Fortran 1.6.0 Change
e1954db00 A4 PR 288 BounceAveraging.jl: for v∥ ≤ 0 past the midpoint, g_wd is held from g_wd[i-1], not g_wb[i-1]. The comment that kept the old copy "for parity" is gone. Julia's index runs one ahead of Fortran's, so Fortran fs(i-2,2) is Julia g_wd[i-1].
d45567603 survey §C analogue of PRs 312/314/334 check_psi_quadrature_convergence stops on a non-finite torque. NaN > tol is false, so a NaN torque used to return without a warning.
f00d20ff8 A11 PRs 289/349 (not ported) Removes calculate_rlar (a placeholder κ average) and calculate_clar (returned 0), with their registry entries, flags and the surface-state fields only they used. _MOVED_KEYS gives the deck message. The faithful port is #622.
4f2ad22ae A14 PR 269 newq0 < 0 errors at configuration. The -1 flip and the sign(newq0) factor are gone from Direct.jl and Inverse.jl.
e812b0e6c A5 PR 336 The collision model floors T at 1e-6 keV and lnΛ at 1, and errors on a negative input temperature. lnΛ, ν_i and ν_e are pchip interpolants: KineticProfileSplines gains a third type parameter. At a T = 0 node the old code returned ν = 0.
0c038188b A5 PR 300 _cubic_resample becomes _pchip_resample, so the kinetic tables are fit with pchip and cannot overshoot between table points.
5abf20f1b review PR 336 Densities must be positive and temperatures non-negative on the resampled grid, where Fortran also checks them. Before, a zero density gave an infinite collision rate. A one-time warning appears when the T or lnΛ floor engages, as in Fortran. The newq0 comment in every deck now says negative values are rejected.
09e25641d — — Formats the branch with the pinned hooks, on top of #356's sweep.

FastInterpolations' pchip uses the same weighted harmonic-mean interior slopes and three-point end slope as Fortran's spline_fit_pchip.

Regression report

Harness on SDCC (regress-fanout run pr625_e0bb4f7c1): develop at 2d12b7962 vs this branch at e0bb4f7c1.

Case Moved Largest relative change
tj_kinetic_calculated 13/20 9.3e-2 (‖dξ/dψ‖, ‖ξ_α‖); Im et₁ 1.0e-2, NTV torque 7e-4
solovev_kinetic_nuzero 7/14 1.1e-2 (total energy Re et₁)
diiid_slayer_n1_coupled 5/7 6.8e-3 (coupled γ)
diiid_slayer_n1 (+ chi_perp_e, winding) 14/22 each 1.1e-3 (γ, CRF Q_peak)
diiid_slayer_n1_toroidal 3/5 4e-4 (γ)
tj_kinetic_multiion 3/6 5e-5 (NTV tritium)
diiid_error_field 5/10 1e-4 (EFC residual torque vs rotation)
diiid_kinetic_multi_n 4/8 4e-5 (Im et₁); ODE step count 1.7e-2
diiid_n1 31/53 5e-7 (PE surface energy); ODE step count 5e-3
diiid_truncation_n1 2/12 2e-7 (energies)
tj_kinetic_ntv 1/7 2e-5 (NTV fgar)
the other 16 cases 0 — (gal_resistive_pe N/A by design)

Attribution: a per-commit run (pr625_attrib) of tj_kinetic_calculated and diiid_slayer_n1_coupled at 6e2796718 (develop merged, before the review fixes), e38b2c995 (bounce-point order) and 50a7c8efb (scaling-factor and grid checks) is bit-identical across all three, so those two cases' moves come from the kinetic corrections in the earlier commits, not from the review fixes. The other cases were not split by commit.

Why each moves:

  • Kinetic cases: the pchip resample of the kinetic tables (Fortran PR 300) changes the profiles between table points, and the ω_D hold (A4) and collision floors change the bounce and collision terms.
  • SLAYER and the two-pass-grid decks (diiid_n1, diiid_truncation_n1, diiid_error_field): the SLAYER layer builds its own cubic splines of the raw table, so it does not read the resampled profiles directly; these decks place two-pass ψ-grid knots at kinetic-resonance surfaces computed from the NTV view, so the grid, and everything on it, shifts slightly.

Notes for reviewers

🤖 Generated with Claude Code

https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF

logan-nc and others added 6 commits November 1, 2026 17:33
…ast the bounce midpoint

Fortran GPEC 1.6.0 fixed the v∥ ≤ 0 fill that copied the bounce-frequency
column into the precession column; Julia had kept the old copy for parity.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF
…urning it without a warning

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF
…thods

RLAR used a placeholder κ average in place of the Fortran κ integral and CLAR
returned zero, so neither matched Fortran GPEC. A deck setting rlar_flag or
clar_flag now stops with a pointer to fgar_flag or tgar_flag.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF
GPEC requires q > 0, matching Fortran GPEC 1.6.0: newq0 = -1 no longer flips
the sign of q, and any negative newq0 stops at configuration.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF
…garithm and take monotone collision-rate slopes

Match Fortran GPEC 1.6.0: the collision model floors T at 1 meV and lnΛ at 1,
so a cold edge node gets a large finite ν instead of zero, negative input
temperatures stop the run, and lnΛ, ν_i and ν_e are pchip interpolants.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF
…tone pchip fits

Match Fortran GPEC 1.6.0, which fits the kinetic table with pchip so a steep
edge cannot overshoot between table points.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF
@logan-nc logan-nc added this to the GPEC v2.0.0 milestone Nov 1, 2026
@logan-nc logan-nc self-assigned this Nov 1, 2026
@github-actions github-actions Bot added bugfix Something was wrong and now is not changed-results Results move or an interface breaks - read before upgrading labels Nov 1, 2026
@logan-nc
logan-nc removed the request for review from jaesun57 November 1, 2026 23:39
… BounceAveraging

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF
logan-nc and others added 4 commits November 1, 2026 19:22
… resampling and warn when a collision floor engages

Densities must be positive and temperatures non-negative on the resampled grid,
as Fortran PENTRC checks them; a zero density used to give an infinite collision
rate. The decks' newq0 comments now match the rejection of negative values.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF
# Conflicts:
#	examples/Solovev_kinetic_NTV_example/gpec.toml
#	examples/Solovev_kinetic_calculated_example/gpec.toml
#	examples/Solovev_kinetic_multiion_example/gpec.toml

@d-burg d-burg left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Adversarial review at b7837705f, by model reviewers working from the code and the Fortran v1.6.0 source. I checked the two pre-existing bugs and items 1–3 by hand. Items are marked ran or traced.

What held up: the A4 index map (hold, restart, backfill, first point past the midpoint); pchip slopes and extrapolation against a Fortran reimplementation (6e-16 on non-uniform grids); floor logic and constants against read_kin; RLAR/CLAR removal, with no residue; every shipped kinetic file passes the new checks. No tolerance or convergence criterion loosens. The added floors are T ≥ 1e-6 keV and lnΛ ≥ 1, and the density check tightens from < 0 to ≤ 0.

To fix here

  1. The new checks see unscaled profiles, but the model uses scaled ones. KineticProfiles.jl:427-431 validates before the scaling, and KineticForcesStructs.jl:117-118 validates nothing. Ran: temperature_factor = -1 on a 5 eV table stores T = −5 eV and gives ν = 2.7e10 with no error; density_factor = 0 gives ν = NaN; temperature_factor = 1e-4 is floored with no warning. Suggested fix: check after scaling, or reject non-positive factors.
  2. The checks can stop decks that don't use NTV. Traced: a file with n_e = 0 at ψ = 1 now aborts an ideal or SLAYER deck that reads it only for two-pass grid resonance surfaces, through CalculatedKineticMatrices.jl:180-182. This is the same class as #614. Suggested fix: run them only for a kinetic source.
  3. The SLAYER explanation doesn't hold. The layer's n, T and ω_E are cubic splines of the raw table (KineticView.jl:18, Utilities/KineticProfiles.jl:47-50), so they don't read the resampled profiles. The resample moves the profiles by at most 2.8e-5 on the DIII-D-like file, so γ moving 1.4e-3 is unexplained. The two-pass grid may be the cause, and a per-commit run would settle it. NTV and SLAYER also now interpolate one table differently.
  4. A4 has no test. No test references g_wd, so reverting e1954db00 should pass (traced). Ran: the hold fired 0 times on single-well model surfaces, so the attribution to "particles reflecting past the midpoint" needs the per-commit run too.

Pre-existing bugs

  • Axisymmetric/Inverse.jl:288, in the block this PR edits: f0 reads column 2 (μ₀p), where Fortran inverse.f:247 uses column 1 (F). With newq0 > 0 on the inverse path, F is barely rescaled, yet q0 = newq0 is still assigned. The new test covers only Direct.jl. Traced. Fix: use column 1.
  • BounceAveraging.jl:873,895: bounce points sort descending, while Fortran's spline_roots sorts ascending. _find_deepest_well pairs neighbors, so with 4 roots every pair spans a barrier. Ran: on a model surface with B ∝ (1 + 0.1cos(4πθ + 1.1))/(1 + 0.05cos2πθ) at λ = 0.9488, the roots are [0.918, 0.882, 0.536, 0.314] and the window is (0.882, 1.536). The hold fires on 45 of 128 pitch points, and 0 of 128 once sorted ascending. runtests_kinetic.jl:303-304 pins the wrong order. We didn't run a harness deck, so whether a shipped equilibrium reaches this is open.

Minor

  • maxlog = 1 is per process and site, so a second file's floor warning is silent (ran).
  • "Cannot overshoot" holds only for the resample. A linear ramp to T = 0 passes the checks, but Ti_spline reaches −12 to −93 eV, and Utils.jl:76 and BounceAveraging.jl:325 take its square root (ran to the negative T, traced after).
  • A NaN torque didn't return silently before this PR: QuadGK's evalrule throws DomainError first (ran). The new check is reachable only by overflow.
  • Untested: floor warnings and the rlar_flag/clar_flag entries, and rlar_flag = false also stops a run. The cold-edge test's (1 − ψ)² looks reproducible by a cubic, so 0c038188b may be unguarded.

This is not an approval. A third-party human reviewer's approval is still required before merge.

@matt-pharr
matt-pharr added this pull request to stack #637 November 2, 2026 09:37
matt-pharr and others added 3 commits November 2, 2026 12:15
…the inverse path and require a positive CHEASE q

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…-well search pairs neighbouring roots

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
matt-pharr and others added 2 commits November 2, 2026 12:31
…d keep unusable profiles out of the two-pass grid

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@matt-pharr
matt-pharr added this pull request to the merge queue Nov 2, 2026
Merged via the queue into develop with commit 2e54c1e Nov 2, 2026
19 of 20 checks passed
@matt-pharr
matt-pharr deleted the bugfix/fortran-parity-kin branch November 2, 2026 12:28
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

bugfix Something was wrong and now is not changed-results Results move or an interface breaks - read before upgrading

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants