KineticForces: QuadGK adaptive quadrature for energy/pitch/psi integration#220
Closed
logan-nc wants to merge 1 commit into
Closed
KineticForces: QuadGK adaptive quadrature for energy/pitch/psi integration#220logan-nc wants to merge 1 commit into
logan-nc wants to merge 1 commit into
Conversation
…itch/psi integration
Replace ODE accumulation (Tsit5) with purpose-built adaptive Gauss-Kronrod
quadrature (QuadGK.jl) for all three KineticForces integration levels.
Energy integration uses a u-substitution u=1-exp(-x) mapping [0,∞)→[0,1)
that absorbs the Maxwellian weight exactly, eliminating the xmax=72
truncation. Resonance poles are handled analytically via Sokhotski-Plemelj
decomposition: the singular part is subtracted and its contribution computed
via the complex logarithm formula, which handles both collisional (ν>0) and
collisionless (ν→0) regimes without the ximag complex contour workaround.
Pitch integration passes the trapped-passing boundary (λ=B₀/Bmax) as a
QuadGK breakpoint to handle the leff/nueff discontinuity. Psi integration
uses QuadGK for scalar torque methods; matrix methods remain on the ODE
path.
Both paths are switchable via ctrl.integration_method ("quadgk" or "ode").
Default is "quadgk". All existing tests pass; 20 new unit tests added
including CGL analytical verification, ODE cross-validation, and resonance
root finding. Regression: 0 diff on solovev_kinetic_calculated and
solovev_n1.
Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
logan-nc
added a commit
that referenced
this pull request
May 23, 2026
…lytic resonance poles Port the QuadGK energy-integration robustness from PR #220 into the KineticForces energy integrator, with a corrected pole decomposition. - Map x = E/T to u = 1-exp(-x) on [0,1), absorbing the Maxwellian weight into du, so the integral covers [0,infinity) without an xmax cutoff. - Remove each resonance pole analytically by a Sokhotski-Plemelj decomposition: subtract R/(u-u_pole) and add back the elementary integral R*(log(1-u_pole)-log(-u_pole)). The same complex pole u_pole (x_pole = x_res - i*nu/Omega') is used in both the subtraction and the add-back, so the decomposition is exact and converges for nu > 0 as well as nu -> 0 (PR #220 subtracted at the real u_res but added at the complex u_pole, which diverges for collisional cases). - Collisionless limit uses the causal Sokhotski-Plemelj branch. - CGL has no resonance denominator: integrate the physical x-space integrand directly over [0,infinity) via QuadGK. - ximag is accepted for signature compatibility but is now unused. Adds find_resonance_energies and unit tests: resonance-root cases, a collisional cross-check against direct x-space quadrature, and a collisionless-limit continuity check. 110 kinetic unit tests pass. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
Collaborator
Author
|
Superseded by #224. #224 ported the three robustness techniques from this PR into
Two additional NaN guards were needed in the calculated-source pipeline that this PR did not exercise: Closing as superseded; the |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
u = 1 - exp(-x)mapping[0,∞) → [0,1)absorbing the Maxwellian weight, eliminatingxmax=72truncationν > 0) and collisionless (ν → 0) withoutximagworkaroundλ = B₀/Bmax) passed as QuadGK breakpoint for pitch integrationctrl.integration_method("quadgk"default,"ode"fallback); matrix methods stay on ODEKey changes
Project.tomlEnergyIntegration.jlfind_resonance_energies,collision_frequency,energy_numerator,integrate_energy_quadgkPitchIntegration.jlintegrate_pitch_gar_quadgk,integrate_pitch_quadgkCompute.jlintegrate_psi_quadgk, dispatch incompute_torque_all_methods!Torque.jlintegration_methodkwarg throughtpsi!→calculate_garKineticForcesStructs.jlintegration_method,xmax,ximagtoKineticForcesControlruntests_kinetic.jlTest plan
solovev_kinetic_calculated— 15 quantities, 0 diff vskineticsolovev_n1— 20 quantities, 0 diff vskinetic🤖 Generated with Claude Code