Skip to content

KineticForces - BUG - nutype="zero" crashes with BoundsError when the resonant integrand is non-finite and no pole was subtracted #428

Description

@logan-nc

Symptom

A DIII-D-like kinetic-calculated run with [KineticForces] nutype = "zero" (collisionless) at nominal rotation dies during kinetic matrix formation:

ERROR: TaskFailedException
  nested task error: BoundsError: attempt to access 0-element Vector{ComplexF64} at index [1]
   [3] ... src/KineticForces/EnergyIntegration.jl:227 [inlined]
   [4] ... src/KineticForces/EnergyIntegration.jl:268
   [5] QuadGK evalrule ...

Reproduced on both a plain run and a variant; the same deck with nutype = "harmonic" runs fine, and nutype = "zero" plus toroidal_rotation_factor = 0.2 also runs fine — so it needs the collisionless operator at nominal rotation.

Mechanism

_resonant_integrand (EnergyIntegration.jl:~215-234) has a non-finite-value fallback that substitutes the nearest pole's analytic regular part:

if !isfinite(val)
    k = 1
    for j in 2:npole
        abs(x - real(x_poles[j])) < abs(x - real(x_poles[k])) && (k = j)
    end
    val = _real_pole_regular_part(real(x_poles[k]), ...)   # x_poles[1] on an empty vector

The fallback assumes at least one pole was found. But the caller builds the same closure for the npole == 0 branch (:~268) and integrates it:

integrand = x -> _resonant_integrand(x, p, residues, x_poles, npole, leff, wb, n, wd)
if npole == 0
    val, _ = quadgk(integrand, 0.0, X_ENERGY_MAX; ...)

With nutype = "zero" there is no collisional broadening, so Ω can round to exactly 0 at a quadrature node even when no pole was registered in [0, X_ENERGY_MAX]val becomes non-finite, the fallback runs with npole = 0, and x_poles[1] throws.

Suggested fix

Guard the fallback on npole > 0, and decide what the physical value is when the integrand is non-finite with no subtracted pole (returning 0 for that node, or excluding it, is a modelling choice worth stating in the annotation rather than implying). A regression case exercising nutype = "zero" on a DIII-D-like equilibrium would cover this — the existing solovev_kinetic_nuzero case passes, so Solovev alone does not reach the branch.

Not introduced by the issue #376 kinetic stack: EnergyIntegration.jl has no commits on that branch relative to develop.

Activity

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions