KineticForces - BUG FIX - Collisionless (nutype=zero) energy integral NaN + nuzero sentinel (#281)

[logan-nc]

Summary

Fixes the collisionless (nutype="zero") energy-integral NaN tracked in #281 (deferred out of #280), and adds the end-to-end nuzero regression sentinel that was blocked on it.

Stacked on #280 (test/kinetic-pv-regression-validation) — it depends on that PR's nutype threading. Base will be retargeted to develop once #280 merges.

The bug

The collisionless energy integral (src/KineticForces/EnergyIntegration.jl) produced NaN / a QuadGK DomainError whenever a resonance pole fell in the Maxwellian tail (x_res ≳ 35). Two mechanisms, both rooted in the u = 1 − exp(−x) substitution compressing the real-axis pole toward u = 1:

• u_res = 1 − exp(−x_res) rounds to 1.0, so R ∝ (1−u_res) = 0 while log1p(−u_res) = −Inf → 0·(−Inf) = NaN;
• and, in real x-space, the resonance denominator Ω(x) = leff·wb·√x + n·(we+wd·x) — a difference of O(10⁴)-magnitude terms — rounds to exactly 0.0 over a ~ULP-wide window (granularity ulp(1e4) ≈ 1.8e-12), where the physical term N·exp(−x)/(iΩ) is 0/0. Neither a pole-location Newton refinement (Ω already rounds to 0) nor a one-ULP nudge (the whole window is flat-zero) escapes it.

The fix

Split the energy integral by collisionality (integrate_energy dispatches on _energy_collision_frequency(1.0,p) == 0.0):

• Collisional (ν > 0): unchanged — keeps the u-space [0,1) integral with the off-axis complex-pole add-back. The harness confirms byte-identical output (see below).
• Collisionless (ν ≡ 0): new real-x-space branch (_integrate_energy_collisionless) over [0, xmax = 72], matching Fortran PENTRC energy.f90. Each real pole is removed by analytic principal-value + residue (causal ν→0⁺ branch ∓iπ·sign(Ω′)), an improvement over the deprecated ximag contour offset. Inside the flat-Ω window the integrand is replaced by the analytic Laurent regular-part limit[h′ − h·Ω″/(2Ω′)]/(iΩ′) (h = N·e^{−x}), via new _collisionless_regular_part + _energy_numerator_deriv.

Implementation note: the x-space route was chosen to keep the real-axis pole well-conditioned (and faithful to Fortran), but it still required the analytic regular-part limit because Ω's catastrophic cancellation makes a multi-float flat-zero window — there is no way around handling the singularity analytically there.

Validation

• Physics:fortran-physics-reviewer independently re-derived the residue, the causal ∓iπ·sign(Ω′) branch, the regular-part Laurent limit (the crux), the numerator derivative, and the xmax=72 truncation — all match the code (PASS), against pentrc/energy.f90 and Logan & Park 2013 Eq. (8).
• Regression harness (regression-guardian):diiid_n1, solovev_n1, solovev_multi_n, solovev_kinetic_calculated are 0.0-diff vs develop (81 quantities, zero regressions — the collisional/ideal paths are byte-identical).
• nuzero sentinel: Solovev calculated collisionless et[1] = 15.885 − 0.482i. Re sits on the harmonic 15.888 to 0.02% (collisions are a small perturbation to the real energy); Im is the resonant ∓iπ damping, smaller than the harmonic −0.711 as expected. Runs in ~68 s (≤ 2 min budget).
• 124 runtests_kinetic.jl unit tests pass (incl. a new tail-pole case in the ν→0⁺-limit test); x_res sweep 0.5…800 is NaN-free; runtests_fullruns.jlex6 added.

Changes

• src/KineticForces/EnergyIntegration.jl — collisionless x-space branch + regular-part limit + dispatch.
• test/test_data/regression_solovev_kinetic_nuzero/ — new fixture ([KineticForces] nutype="zero").
• regression-harness/cases/solovev_kinetic_nuzero.toml — new harness case.
• test/runtests_fullruns.jl (ex6), test/runtests_kinetic.jl (tail-pole case).

Closes#281.

🤖 Generated with Claude Code

[logan-nc]

@claude review this pr

[claude]

Code Review

Overview

This PR fixes a real NaN bug (nutype=\"zero\" energy integral blowing up for resonances deep in the Maxwellian tail) by splitting the collisionless path into a dedicated real-x-space integral. The mathematical approach (PV + residue + Laurent regular-part limit) is sound and matches the Fortran PENTRC convention. Good test coverage: new unit test for the tail-pole case, a new end-to-end fixture, and a regression harness case.

---

Issues

1. Dispatch on function call instead of field value (EnergyIntegration.jl, new dispatch block in integrate_energy)

if_energy_collision_frequency(1.0, p) ==0.0return_integrate_energy_collisionless(p, leff, wb, n, we, wd, atol, rtol)
end

Since p.nutype is a String field directly accessible on EnergyParams, and _energy_collision_frequency is just a switch on that string, this is an indirect way to write p.nutype == "zero". The comment tries to justify the indirection ("ν vanishes identically iff it vanishes at any x"), but it actually introduces a subtle assumption: that any collision frequency returning 0.0 at x=1.0 is identically zero everywhere. That is currently true only because nutype=\"zero\" returns a literal 0.0, not by general contract. Use p.nutype == \"zero\" directly — it is explicit, no floating-point comparison, and reads as an intent check rather than a numerical probe.

---

2. Reactive NaN guard in the integrand closure (_integrate_energy_collisionless)

if!isfinite(val)
# find nearest pole k, replace val with _collisionless_regular_part...end

This works in practice (the only non-finite path is the flat-Ω window around the pole), but it is reactive: it intercepts NaN after the fact and silently replaces val with the regular part of the nearest pole. If a different, unexpected source of NaN were introduced later (e.g., in _energy_integrand_real under some edge-case parameter), the guard would silently return a wrong-but-finite value with no indication anything went wrong.

A proactive threshold check is more explicit:

k =argmin(abs(x - xp) for xp in x_poles) # or the existing loopifabs(x - x_poles[k]) < threshold
val =_collisionless_regular_part(x_poles[k], p, leff, wb, n, wd)
...else
val =_energy_integrand_real(x, p)
...end

The threshold would need to cover the flat-Ω ULP window — roughly eps(Float64) * abs(x_poles[k]). This approach also avoids evaluating _energy_integrand_real at all for the degenerate points, removing the NaN generation entirely.

---

3. Nested @test inside @test begin ... true end (runtests_fullruns.jl, ex6 block)

@testbegin...@testisfinite(real(et[1]))
@testreal(et[1]) ≈15.885 rtol =0.01...trueend

If any inner @test fails, the outer block still returns true and records an outer PASS alongside the inner failures. This is a pre-existing pattern in the file (ex1–ex5 do the same), so not a regression introduced here, but it means the outer @test is vacuous. The inner assertions do register independently in the test log, so failures are visible — but the signal is muddied. Worth fixing consistently across the file when convenient.

---

Minor

• unique(sort(x_poles)) — sort already returns a copy in Julia so x_poles is safe, but the comment explaining why ("Sort a copy — x_poles must stay aligned with residues") correctly documents the intent. No issue, just confirming it is safe.
• The new gpec.toml fixture omits [PerturbedEquilibrium] and [ForcingTerms] sections — consistent with the existing regression_solovev_kinetic_calculated/gpec.toml, which does the same. Fine.
• X_ENERGY_MAX = 72.0 as a module-level const with a comment explaining the Fortran provenance is the right approach.

---

Summary

The math and the fix are correct. Items 1 and 2 are the issues worth addressing before merge: item 1 is a one-line change, item 2 requires a bit more refactoring but eliminates a potential silent-wrong-answer path.

[logan-nc]

Follow-on: unified onto a single x-space integrator (0d7c02c)

Per reviewer feedback that the two parallel integration paths were a maintainability cost, the collisional (ν>0) and collisionless (ν=0) cases are now one real-x-space PV+residue integrator, _integrate_energy_resonant. The pole is x_pole = x_res − i·ν/Ω′ (off-axis for ν>0, on-axis for ν=0), so the collisionless case is the exact ν→0 limit of one formula. The causal ∓iπ·sign(Ω′) retarded branch now falls out of the signed zero of pole_offset = ν/Ω′ (verified to 1e-14 for both signs of Ω′); the flat-Ω regular-part limit only ever fires for ν=0. The ~40-line u-space collisional block is deleted (net −47 lines).

This supersedes two claims in the PR description above:

• The collisional ν>0 path is no longer "unchanged / byte-identical" — it now routes through x-space.
• Consequently solovev_kinetic_calculated (harmonic) moves by ~1e-6 relative on et[1] (per-integral reassociation ~5e-12 amplified through the ill-conditioned eigenproblem; harness reports 0.00%, well inside the fullruns rtol brackets). The 3 ideal harness cases remain 0.0-diff. Collisionless ν=0 output is unchanged from the prior commit.

Validation:fortran-physics-reviewer PASS on residue, add-back, the signed-zero causal branch, the ν=0-only regular-part guard, and truncation — confirmed faithful to Fortran energy.f90 and Logan & Park 2013 Eq. (8). 126 runtests_kinetic (+ a new Ω′<0 causal-branch test) and 24 runtests_fullruns green; regression harness shows the 3 ideal cases 0.0-diff.

Post-merge note: re-cache the solovev_kinetic_calculated harness baseline at the new x-space values (regress --cases solovev_kinetic_calculated --refs local --force) so future comparisons start from the unified path.

[logan-nc]

Follow-on: kinetic regression on a mode-converged deck; fullruns → smoke tests (closes #273, 8bf5d01)

Separated the two jobs the nerfed Solovev-kinetic cases were conflating:

Regression harness — the two kinetic cases now point at new full-band example decks examples/Solovev_kinetic_{calculated,nuzero}_example (mpsi=16, delta_mlow=delta_mhigh=8 → mpert=32). This pins the mode-converged eigenvalue (calc et[1]≈1.84−1.41i, nuzero ≈1.72−1.08i) instead of the mpert=8-truncated 15.9. Mirrors the existing pattern (ideal harness cases point at examples/ decks; fullruns use the nerfed test/test_data fixtures). Runtime ~105–120s/case (within the ~2-min budget).

The mode band — not mpsi — is the dominant lever: restoring delta_m=8 moves et[1] ~8× (mpert=8→32). The pinned mpert=32 value is a Julia self-consistent anchor; the Fortran cross-check in benchmark_solovev_kinetic_stability.jl is locked to mpert=8 (15.888 vs 15.619, 1.7%). (The auto-mpsi log_asymptotic path is a no-op for analytic Solovev — no separatrix → flat slope estimate → near-uniform ~16-surface grid.)

Tests (fullruns ex3–ex6) — relaxed to run-through + physical-sanity smoke tests on the nerfed fixtures (isfinite, total energy > 0, damping imag < 0). Dropped the numeric value pins, including the platform-fragile imag(et[1]) ≈ −0.711 rtol=0.08 that failed on macOS aarch64 (−0.856). Closes #273. The precise collisionless-integral correctness (#281) stays guarded by the deterministic runtests_kinetic.jl unit tests.

Validation: 126 runtests_kinetic + 17 runtests_fullruns green; harness runs both new decks within budget at the expected mpert=32 values.

So this PR now: fixes the collisionless NaN (#281), unifies the energy integral onto one x-space path, and re-architects the kinetic coverage into fast smoke tests + mode-converged harness pins (closing #273).