Skip to content

Equilibrium - FEATURE - Accept supplied Coulomb logarithm and collision rates in kinetic profiles - #627

Merged
matt-pharr merged 19 commits into
bugfix/fortran-parity-kinfrom
feature/kinetic-supplied-rates
Nov 2, 2026
Merged

matt-pharr merged 19 commits into
bugfix/fortran-parity-kinfrom
feature/kinetic-supplied-rates

Conversation

@logan-nc

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

Copy link
Copy Markdown
Collaborator

Release note

Kinetic profile files can now supply the Coulomb logarithm and the ion and electron collision rates, as HDF5 datasets lnLambda, nu_i and nu_e, or as a 9-column ASCII table read with kinetic_format = "rates-v1" in [Equilibrium] (API keyword format). They replace the NTV calculation's internal collision model; the tearing layers' resistivity keeps its own lnLambda_form closure. This is the Julia counterpart of Fortran GPEC 1.6.0's kinetic_format = "rates-v1".

Closes #624: it covers every item, in code or with the evidence below.

Closes #591: this is the last PR of the Fortran 1.6.0 parity series (#625, #626, #627). The feature-parity items #591 listed have their own issues (#622, #629, #630, #633, #429, #623; #628, #631 and #632 were closed as deliberate or covered), and the four unfinished parity checks moved to #640.

What changed

Commit #624 item Change
a12d936b4 Supplied collision rates (Fortran PR 345) See the bullets below.
3883dba76 Docs pass (Fortran PRs 282 and 266) conventions.md now states the laboratory direction of positive ω_E from the toroidal map φ_CCW = −h(2πζ + ν) used by sample_boundary_grid: positive ω_E is co-current exactly when B_t is clockwise, whatever the helicity. Fortran made the same correction. Julia's ℓ_eff = ℓ + nq already had Fortran's corrected sign.
a84e9731d review A ragged ASCII row is an error instead of being silently dropped.
25033d47e Inner quadrature accounting check_inner_quadrature warns once per integral kind and Julia session when a pitch or energy integral stops short of its tolerance or is not finite; see below.
d42c17bb4 Inner quadrature accounting The NTV rotation scan keeps one quiet logger, so a maxlog warning prints once per scan rather than once per scan point.
0787c5ebb, da929129f Overlap check (A13) The perturbed-equilibrium guide shows how to decompose any of the coupling matrices over any window of surfaces with the analysis helper, which replaces Fortran's per-model and local-window SVD outputs (#631, closed as deliberate).
83396e096 review Docstrings and docs say which rates bypass the model and are not rescaled.
f5b64a933 review A supplied lnΛ no longer feeds the resistivity closure. A rates-v1 lnΛ is the ion-type value PENTRC's collision model uses, not the electron lnΛ Spitzer η wants, and Fortran feeds it to nothing else.
3ede86983 review kinetic_format selects the ASCII layout, so a 9-column file is no longer read as rates by its column count alone. The column count must match the format.
0ed35062c review Harness case tj_kinetic_ntv_rates; see below.

What the supplied rates do

  • Inputs: three optional KineticProfileData fields and HDF5 datasets, lnLambda, nu_i and nu_e, each one row in _KINETIC_H5_PROFILES like chi_e; the ASCII reader takes Fortran's 9-column table under kinetic_format = "rates-v1".
  • Validation: values must be positive; an all-zero dataset counts as absent.
  • Interpolation: pchip in log space, as Fortran, so values stay positive between table points.
  • NTV collision model: the supplied values replace the internal ones in kinetic_profile_splines and resolve_ntv_species.
    • A supplied ν_i is the first main ion's rate. Other species keep the internal model's (z_s/z₁)²·√(m₁/m_s) ratio to it, which is exact: n_main, lnΛ, T_i and zpitch are common to every species.
    • A supplied lnΛ without ν feeds the internal ν.
    • The density and temperature scaling factors do not rescale supplied rates; collisionality_factor still multiplies them.

#624 items closed without code, and the inner-quadrature check

  • Overlap check (A13, Fortran PRs 330 and 185).
    • Fortran's SVD was in a weighted basis, so its exported vectors had to be W^H·V.
    • Julia takes the SVD of C_resonant_area_weighted_field, which acts on b̃ (the root-area-weighted field), and the overlap V^H·b̃ uses that same basis.
    • test/runtests_dominant_coupling.jl:97 already checks that U * (σ .* overlap) reproduces the solve's resonant field to 1e-8.
    • forcing_overlap carries units (T), which was PR 185's other point.
  • Inner quadrature accounting (commit 25033d47e).
    • The gap: the pitch and energy integrals discarded QuadGK's error estimates, so an integral that stopped short of its tolerance returned its estimate silently.
    • The fix: check_inner_quadrature warns once per integral kind (pitch, energy, CGL energy) when the error estimate exceeds max(atol, rtol·|I|), which happens when QuadGK stops at its evaluation cap or on endpoint roundoff, or when the result is not finite.
    • Cost: one norm per integral.
  • Spline end slopes (A15, Fortran PR 308). cubic_interp + deriv1 end slopes stay within a few times the round-off limit (machine epsilon × offset / grid spacing) up to an offset of 1e10. There is no cancellation like the one Fortran fixed.

Regression report

Harness on SDCC (regress-fanout run pr627_2667d1294): #625 at e0bb4f7c1 vs this branch at 2667d1294.

Notes for reviewers

🤖 Generated with Claude Code

@logan-nc logan-nc added this to the GPEC v2.0.0 milestone Nov 2, 2026
@logan-nc logan-nc self-assigned this Nov 2, 2026
@logan-nc
logan-nc requested review from d-burg and jaesun57 November 2, 2026 00:02
@github-actions github-actions Bot added the feature New capability label Nov 2, 2026
@logan-nc
logan-nc requested review from matt-pharr and removed request for jaesun57 November 2, 2026 00:13
logan-nc and others added 3 commits November 1, 2026 19:23
…on rates in kinetic profiles

Kinetic files may carry lnLambda, nu_i and nu_e (HDF5 datasets, or a 9-column
ASCII table as Fortran PENTRC rates-v1). Each replaces the internal collision
model, interpolated in log space; a supplied lnLambda also replaces the
resistivity closure's, so KineticForces and the tearing layers see one value.
A supplied nu_i is the first main ion's rate, and other species keep the
model's (z_s/z_1)^2 sqrt(m_1/m_s) ratio to it.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF
… the toroidal map and B_t

Mirrors the corrected Fortran GPEC sign-convention statement: positive omega_E
is co-current exactly when B_t is clockwise, not according to helicity alone.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF
… supplied rates meet the closure and scaling factors

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF
@logan-nc
logan-nc force-pushed the feature/kinetic-supplied-rates branch from d6878aa to a84e973 Compare November 2, 2026 00:24
logan-nc and others added 2 commits November 1, 2026 20:57
…hort of its tolerance or is not finite

The inner QuadGK calls discarded their error estimates, so an integral that
stopped at QuadGK's evaluation cap returned its estimate silently; the psi
quadrature check sees only the values. Each kind now warns once.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF
…n so maxlog warnings print once per scan

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF
… over any window of surfaces

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF
…eld the real one carries

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

@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 0404a9028, by model reviewers working from the code and the Fortran v1.6.0 source. I checked items 1, 2 and 5–7 by hand and reran the A15 numbers. Items are marked ran or traced.

What held up: column order and units against inputs.f90; 7- and 8-column files error; the species map and D/T ratios; n/T factors leave supplied ν unscaled, as Fortran; all four closure callers pass one lnΛ, with none left on lnLambda_form; the tolerance formula, zero allocations on the passing path, and thread safety (3 warnings in 20 of 20 repeats); the ω_E direction; the stamp, since src/ is unchanged from d42c17bb4. No tolerance loosens.

To fix here

  1. Any 9-column ASCII file becomes rates. KineticProfiles.jl:107 goes by column count. Ran: a table with χ values near 1 in columns 7–9 is accepted. Traced: they then act as lnΛ ≈ 1 and ν ≈ 1/s in the torque and η, with no warning. Fortran's default legacy format ignores those columns, and Julia errored before this PR. #624 called for a kinetic_format key, so "covers every item" isn't quite true. Suggested fix: add the key next to kinetic_group.
  2. "Closing it drops nothing" doesn't hold. These appear only in #591: the open _rs_at / _da_dpsi_at normalization check (Fortran PR 249), qa = qlim not traced (Fortran PRs 225/234), the direct path not inspected (Fortran PRs 298/335), and the pmodb gauges docs (Fortran PR 339). #632 was closed as covered by #465, which is still open with changes requested. Suggested fix: move them to a v2.1 issue before #591 closes.
  3. One lnΛ feeds both the ion NTV model and the electron η closure, with no opt-out, and no Fortran layer reads it. At n_e = 5e19 and T_e = 1 keV, the NTV form gives 17.65 and the closure gives 15.14, so a PENTRC-style file raises Spitzer η about 17% (traced). The body intends this, but it's a physics decision for the maintainers.
  4. The harness covers none of this. No case, example or fixture carries rates, so "27 unchanged" says nothing about the feature. The only closure test goes through layer_eta, so dropping lnLambda= at TearingProblem.jl:64 or LayerOverlap.jl:415, or lnLamb= at LayerParameters.jl:290,297, would pass (traced). Also untested: electron=true (KineticProfiles.jl:333), ν_e-only supply, scale factors and an end-to-end torque. Suggested fix: one harness case with a rates file.
  5. Partial supply is silent. At PlasmaEquilibrium.jl:218, an all-zero column reverts that rate to the model alongside the supplied ones, where Fortran stops (inputs.f90:231-232). Nothing records which rates ran, in the log or in gpec.h5.
  6. maxlog is per process, not "per run". check_inner_quadrature is silent on a second failing run in the same session (ran), so a Julia sweep reports only its first failure.
  7. The quiet logger is wider than one scan. NTVLimits.jl:534 spans all coils and both perturbations (:558, :580-581). The warning names no coil, ψ or rotation, and EFCCoupling has no failure flag (traced).

Open questions

  • A15: "no cancellation" isn't shown, because the test data's noise hides it. With exact constant 1e10, end slopes on a uniform 101-point grid are 6.8e-5 and 5.6e-4, and with 3.3e6 they are −3.1e-8 and 1.2e-7 (ran). Fortran's PR 308 differences from fs(0), so it gives 0 (traced). This is within about 3× round-off, so low priority, but #624 says to file upstream in this case.
  • Exit crash: nothing shows the free(): invalid pointer came from a native library's cleanup, and there's no backtrace. An earlier out-of-bounds write fits equally. The harness did record it as failed, and the reported numbers are from clean reruns. A --check-bounds=yes run would help.

Minor

  • 131e17c06 and d6878aa26 are stale; they're now a12d936b4 and 3883dba76.
  • QuadGK can also stop before the cap: check_endpoint_roundoff stopped 1/(x − π/4) at 1395 evaluations (ran). The check catches this, so only the docstring is wrong. Its non-finite branch is unreachable from QuadGK, as on #625.
  • With lnΛ supplied and ν absent, a cold edge logs "flooring the Coulomb logarithm" though no floor applies (KineticProfiles.jl:35).
  • 1.0D+19 exponents, legal in Fortran, stop with a misleading "fewer columns" error (ran).
  • perturbed_equilibrium.md says five matrices are written, but C_penetrated_area_weighted_field is skipped when ffs.bpen is empty. No test covers a custom-window overlap, and runtests_dominant_coupling.jl:118 tests the core window in memory, not from the saved file.

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 matt-pharr 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.

lgtm. made some fixes and merged develop up the stack.

@matt-pharr
matt-pharr added this pull request to the merge queue Nov 2, 2026
Merged via the queue into develop with commit 329ee59 Nov 2, 2026
17 of 19 checks passed
@matt-pharr
matt-pharr deleted the feature/kinetic-supplied-rates 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

feature New capability

Projects

None yet

3 participants