Uh oh!
There was an error while loading. Please reload this page.
KineticForces - Validate calculated kinetic-DCON vs Fortran (#227); thread nutype, guard torque diagnostics - #280
Conversation
…guard NTV torque diagnostics The self-consistent kinetic_source="calculated" path hardcoded nutype="harmonic": compute_kinetic_matrices_at_psi! never accepted nutype/f0type, so a user setting them in [KineticForces] was silently ignored. Thread nutype/f0type/nufac from kf_ctrl through to kinetic_energy_matrices_for_euler_lagrange! (defaults match the prior hardcoded values, so all existing cases are byte-identical). Also guard the standalone NTV torque diagnostics on @isdefined(pe_state): they contract kinetic operators against perturbed-equilibrium displacements, which the calculated path does not produce, so running them there crashed. Skip with an info message instead. Regression harness (develop vs local): diiid_n1, solovev_n1, solovev_multi_n, solovev_kinetic_calculated all 0.0 diff; 120 kinetic unit tests pass. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…] vs Fortran kinetic DCON (#227) Add benchmarks/benchmark_solovev_kinetic_stability.jl, which generates a Fortran kinetic-DCON deck matched to the Solovev calculated regression fixture (same Solovev params, kinetic profile, mode band, nl, collision operator, grid, edge truncation), runs Fortran dcon, and compares its W_t[1] against Julia's XiNorm eigenmode_energies[1]. Result: Julia 15.888-0.711i vs Fortran 15.619-0.660i (Re 1.7%, Im 7.6%) — passing the DIIID kinetic acceptance thresholds. This physics-validates the PV+residue energy integral (issue #227): the pinned calculated-path et[1]=15.888 is correct; the pre-rewrite ximag value (34.176) was not. Update the ex5 fullruns comment to cite this Solovev reference. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
The collisionless (nutype="zero", or any nu==0 config) energy integral produced NaN/DomainError for resonance poles in the Maxwellian tail (x_res >~ 35): under the u = 1-exp(-x) substitution u_res rounds to 1.0, so the residue weight R proportional to (1-u_res) collapses to 0 while log1p(-u_res) -> -Inf, tripping 0*(-Inf) = NaN; in the [35,36] band a Gauss-Kronrod node coincident with the breakpoint gives 0/0 = NaN (QuadGK DomainError). Root cause is the u-substitution compressing the real-axis pole toward u=1. Fix: dispatch the structurally-collisionless case to a dedicated real-x-space integral over [0, X_ENERGY_MAX=72] (matching Fortran PENTRC energy.f90), where the pole x_res in (0,72) stays a well-conditioned O(1) number. Each pole is removed by analytic principal-value + residue (causal nu->0+ branch -/+ i*pi). The collisional (nu>0) path keeps the u-space formulation unchanged (off-axis complex pole, no truncation) - byte-identical output. Verified: x_res sweep 0.5..800 all finite; matches the nu->0+ harmonic limit to ~1e-6; 120 kinetic unit tests pass. Addresses #281. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…nce poles Follow-up to the x-space collisionless energy integral: the real-axis resonance denominator Ω(x) = leff·wb·√x + n·(we+wd·x) is a difference of O(1e4)-magnitude terms that cancel to ~0 at the pole, so it rounds to exactly 0.0 over a ~ULP-wide window (granularity ulp(1e4) ≈ 1.8e-12). Across that window the physical term N·exp(-x)/(iΩ) is 0/0 → NaN (QuadGK DomainError), which neither a Newton pole refinement (Ω already rounds to 0) nor a one-ULP nudge (whole window is flat-zero) can escape. Fix: in the flat-Ω window QuadGK is handed the analytic Laurent regular-part limit of the pole-subtracted integrand, [h′ - h·Ω″/(2Ω′)]/(iΩ′) with h = N·exp(-x) (new _energy_numerator_deriv + _collisionless_regular_part). Other poles' finite subtractions are retained. Collisional (ν>0) path is untouched. Also adds the nutype="zero" regression sentinel: - test fixture regression_solovev_kinetic_nuzero (calculated path, [KineticForces] nutype="zero"), - regression-harness case solovev_kinetic_nuzero, - fullruns ex6 and a tail-pole case in the runtests_kinetic ν→0⁺ limit test. Validation: Solovev calculated et[1] = 15.885 - 0.482i (finite). 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. Runtime 23s. 124 kinetic unit tests pass; x_res sweep 0.5..800 NaN-free. Closes#281. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
logan-nc
commented
Jun 15, 2026
@claude review this pr |
Code ReviewBug fixes (correct)
NTV torque diagnostic guard ( Issues1. 2. Comment block length ( 3. Test coverageNo test explicitly exercises the NTV skip path (KineticForces block present, no PE state). The Benchmark scriptThe deck construction logic (namelist patching, fixture-to-Fortran conversion) is well-structured. The env-var overrides ( |
…residue path Fold the two energy-integral paths (u-space for ν>0, real x-space for ν=0) into a single real-x-space principal-value + residue integrator, `_integrate_energy_resonant`. The pole is `x_pole = x_res − i·ν/Ω′` — off the real axis for ν>0, on it for ν=0 — so the collisionless case is the exact ν→0 limit of one formula rather than a separate branch. The causal ∓iπ·sign(Ω′) retarded branch now falls out of the signed zero of `pole_offset = ν/Ω′` (verified to 1e-14 for both signs of Ω′), and the flat-Ω regular-part limit (`_real_pole_regular_part`, renamed) only ever fires for ν=0. Deletes the ~40-line u-space collisional block (`u_poles`/`u_breaks`, the `u=1−exp(−x)` substitution, the `xr>700` guard). Net −47 lines. The collisional harmonic path now routes through x-space instead of u-space, moving `solovev_kinetic_calculated` et[1] by ~1e-6 relative (per-integral reassociation ~5e-12 amplified through the ill-conditioned eigenproblem; harness 0.00%, inside the fullruns rtol brackets). The 3 ideal harness cases are 0.0-diff. Collisionless ν=0 output is unchanged from #282. Verified: fortran-physics-reviewer PASS (residue, add-back, causal branch, regular-part guard, truncation); 126 runtests_kinetic + 24 runtests_fullruns green; regression harness 3 ideal cases 0.0-diff. Adds an Ω′<0 collisionless causal-branch unit test. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…uns to smoke (closes#273) Split the two conflated jobs the nerfed Solovev-kinetic cases were doing: - REGRESSION (harness): repoint the two kinetic cases at new full-band example decks examples/Solovev_kinetic_{calculated,nuzero}_example (mpsi=16, delta_mlow=delta_mhigh=8 → mpert=32), so the pinned eigenvalue is the mode-converged value (calc et[1]≈1.84−1.41i, nuzero≈1.72−1.08i) rather than the mpert=8-truncated 15.9. Mirrors how the ideal cases already point the harness at examples/ decks while fullruns use the nerfed test/test_data fixtures. Runtime ~105-120s/case (within the ~2-min budget). NB: the 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). - TEST (fullruns ex3-ex6): relax to run-through + physical-sanity smoke tests on the nerfed fixtures (isfinite, total energy > 0, damping imag < 0). Drop the numeric value pins — including the platform-fragile imag(et[1])≈-0.711 rtol=0.08 pin that failed on macOS aarch64 (-0.856). Closes#273. The mode band (not mpsi) is the dominant cost+physics lever here; the auto-mpsi (log_asymptotic) path is a no-op for analytic Solovev (no separatrix → flat slope estimate → near-uniform ~16-surface grid). Verified: 126 runtests_kinetic + 17 runtests_fullruns green; harness runs both new decks within budget at the expected mpert=32 values. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…ix/kinetic-collisionless-nan
…/kinetic-pv-regression-validation
…ation' into bugfix/kinetic-collisionless-nan
…files; collapse nu=0 deck via harness override Make the Solovev kinetic example physically rigorous and remove the duplicate collisionless deck. Profiles: replace the invented parabolic kinetic.dat with profiles tied to the Solovev equilibrium pressure P(psi)=P0*(1-psi), P0~4.27e4 Pa (beta_axis~11%). T(psi)=100+2600*(1-psi) eV (core 2700, edge 100); n=P0*(1-psi)/(2*e*T) (core ~4.9e19) with a small edge density floor so T stays finite. 2*n*T matches P_eq to ~1e-4. Written to the example deck and both kept test fixtures with a provenance header. nu=0 collapse: delete examples/Solovev_kinetic_nuzero_example. Add a generic [overrides] mechanism to the regression harness (CaseSpec.overrides -> _materialize_rundir copies the deck to a temp dir, patches dotted gpec.toml keys, cleans up). solovev_kinetic_nuzero now reuses the calculated deck with KineticForces.nutype="zero". Re-baselined: calculated et[1]~1.894-1.525i, nuzero et[1]~2.333-2.174i. Reader robustness: KineticProfiles readdlm(...; comments=true) so a '#' provenance header cannot widen the parsed matrix and pad the data rows. Deck migration: port the kinetic decks off the removed bal_flag/mer_flag onto local_stability_flag and drop the deprecated delta_mband (the develop merge removed those keys, which was breaking every kinetic run on this branch). Verified: fullruns 17/17 green (ex5/ex6 with new profiles), harness override runs end-to-end with the same equilibrium and distinct collisionless damping. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…100 and harden boundary-pole handling Manual PR review asked whether the real-x-space energy integral's upper limit X_ENERGY_MAX could be chosen more robustly, and whether the approach is robust to a resonance pole sitting at or within ~1e-15 of that limit (the direct analog of the original u→1 NaN this branch fixed). Findings: the integrand AND every resonance pole decay as x^p·exp(-x) (p ≤ 3.5 with the qt heat-flux factor) — at a simple pole Ω′ is O(the drift frequencies), so a dropped pole's residue scales as xr^2.5·exp(-xr), the same law as the smooth tail. exp(-x) does not underflow until x≈746, so the Fortran default 72 was a tolerance edge, not a precision edge. Benchmarking confirmed a tighter, tolerance-derived limit buys no speed (QuadGK accepts the near-zero tail in a single panel; 72→100 is byte-identical eval count). So extend the fixed limit 72 → 100 (100^3.5·exp(-100) ≈ 4e-37, ~30 orders below any tolerance) and document the real justification. A kept pole (xr < X_ENERGY_MAX) is subtracted analytically/exactly and is robust at any proximity to the limit; the xr >= X_ENERGY_MAX guard prevents the only NaN risk, log(X_ENERGY_MAX - x_pole)=log(0), at the exact endpoint. Near-degenerate poles (Ω′→0) are handled by the SINGULAR_EPS guard, not by the limit. Added a regression testset forcing a resonance to exactly X_ENERGY_MAX and ±1e-13 around it, asserting finiteness, continuity across the boundary, and agreement with the ν→0⁺ limit. All 200 kinetic unit tests pass. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…onless-nan KineticForces - BUG FIX - Collisionless (nutype=zero) energy integral NaN + nuzero sentinel (#281)
…nt/perf cleanups, count nstep_total on kinetic paths - benchmark_solovev_kinetic_stability.jl: rename local run->rundir (shadowed Base.run); mirror nufac into the matched Fortran pentrc.in deck. - EnergyIntegration.jl: trim over-long comments (X_ENERGY_MAX, resonant-integral docstring, inline overflow/non-finite notes) to the project conciseness rule; sizehint! the pole/residue buffers and replace the type-unstable breaks... splat into quadgk with an explicit breakpoint-count branch (byte-identical). - GeneralizedPerturbedEquilibrium.jl: shorten the NTV-skip comment to two lines. - Riccati.jl: record integration/nstep_total on the kinetic (Riccati + parallel-FM) paths, which previously reported 0 - count every accepted Riccati callback step and sum each propagator solve's accepted steps from the thread-local proxies after the BVP barrier. Verified: 200/200 kinetic unit tests pass; regression harness et[1] byte-identical on solovev_kinetic_calculated/solovev_kinetic_nuzero; nstep_total now 1808/1815. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Uh oh!
There was an error while loading. Please reload this page.
…erances The energy (x) integration was handed the same atol_xlmda/rtol_xlmda pair as the pitch (lambda) integration that encloses it. Since the pitch integrand IS the energy integral, the outer integrator was asked to resolve its integrand to the same tolerance to which that integrand was itself computed, so it chases the inner integrator's quadrature noise instead of converging. Adds atol_x/rtol_x, defaulting (NaN sentinel) to nested_tolerance_margin = 1e-2 times the pitch tolerances, which extends the nesting rule the struct already documents one level up for rtol_psi vs rtol_xlmda. Wired through both the psi-quadrature path (tpsi!) and the kinetic-matrix path. This deliberately changes default numerical behaviour: shipped decks now integrate the energy variable to atol 1e-10 / rtol 1e-7 rather than 1e-8 / 1e-5, and pay for it in runtime. A deck can set atol_x/rtol_x explicitly to override the derived values, or widen nested_tolerance_margin to recover the old cost. Julia analog of the second commit of Fortran GPEC PR #280. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Summary
Closes the regression-validation gaps in #227 (PV+residue energy-integral pole handling)
and folds in the collisionless energy-integral fix (#282, merged into this branch).
Investigation finding: the PV+residue rewrite #227 proposed is already on develop
(exact analytic Sokhotski-Plemelj;
ximagdeprecated; the proposedx_residue_radiusknob is unnecessary). The open gap was #227's acceptance criterion 2 — the calculated-path
Solovev
et[1]was pinned to a code-captured value rather than a physics baseline.Criterion 2 — physics validation against Fortran kinetic DCON
benchmarks/benchmark_solovev_kinetic_stability.jlgenerates a Fortran kinetic-DCON deckmatched to the Solovev calculated fixture (same Solovev params, kinetic profile, mode band,
nl, collision operator, grid, edge truncation), runs Fortrandcon, and compares itsW_t[1]to Julia'sXiNorm/eigenmode_energies[1]:Passes the DIIID kinetic acceptance thresholds (Re ≤5 %, Im ≤20 %): the PV-era value is
correct; the pre-rewrite
ximagvalue was not.Criterion 3 — collisionless energy integral (now included via #282)
The energy integral (
src/KineticForces/EnergyIntegration.jl) is a single real-x-spacepath over
[0, X_ENERGY_MAX]for all collisionalities (matching Fortran PENTRCenergy.f90). Each resonance polex_pole = x_res − i·ν/Ω′is removed analytically byprincipal-value + residue (Sokhotski-Plemelj); the collisionless case (ν ≡ 0) is the exact
ν→0⁺ limit, with its on-axis pole resolved by the analytic Laurent regular-part limit
[h′ − h·Ω″/(2Ω′)]/(iΩ′). This resolves the collisionless NaN that earlier blocked thenuzerosentinel, which is now added as both aruntests_fullruns.jlend-to-end smoke testand a
solovev_kinetic_nuzeroregression case (the calculated deck run collisionless via anew harness
[overrides]mechanism).Bug fixes
nutype/f0type/nufac—compute_kinetic_matrices_at_psi!never accepted them, so a user setting them in
[KineticForces]was silently ignored(hardcoded
"harmonic"/"maxwellian"). Now threaded fromkf_ctrl; defaults match theprior values, so existing cases are byte-identical.
[KineticForces]was present withoutPE data (the calculated path). Now skipped with an info message.
nstep_totalwas 0 on the kinetic paths — the Riccati callback and the parallel-FMpropagator solves never counted solver steps. Now counted (Riccati callback per accepted
step; parallel-FM chunks via per-proxy accumulation summed after the BVP barrier).
Verification
solovev_kinetic_calculatedet[1] ≈ 1.894−1.525i,solovev_kinetic_nuzeroet[1] ≈ 2.333−2.174i, both matching the pressure-consistentkinetic.dat. The calculated case'sexample_dirnow points at a newexamples/Solovev_kinetic_calculated_example(absent on develop), so the first baselineis captured with
--refs local.Review response (this PR)
run→rundir(shadowedBase.run); mirrorednufacinto theFortran deck.
X_ENERGY_MAX, the resonant-integraldocstring, the NTV-skip comment).
sizehint!the pole buffers and replaced the type-unstablebreaks...splat intoquadgkwith an explicit breakpoint-count branch.🤖 Generated with Claude Code