Summary
The self-consistent kinetic matrix build uses adaptive Gauss–Kronrod (QuadGK) for the
pitch/energy integrals, but the default rtol_xlmda is set to match PENTRC's nominal
value. QuadGK bounds the overall-integral error while PENTRC/LSODE bounds per-step
error, so matching the nominal tolerance over-refines Julia by ~3× at no accuracy benefit.
Where
KineticForces.compute_calculated_kinetic_matrices ->
compute_kinetic_matrices_at_psi! -> integrate_pitch_gar_quadgk (QuadGK.quadgk! with
atol=atol_xlmda, rtol=rtol_xlmda). The control default is rtol_xlmda=1e-5,
atol_xlmda=1e-8.
Measurement (matched grid: mpsi=257, mpert=27, nl=6, 8 threads)
rtol_xlmda |
atol_xlmda |
matrix-build wall |
| 1e-5 (nominal-matched) |
1e-9 |
1668 s |
| 1e-3 (rtol-dominated) |
0.1 |
512 s (3.3× faster) |
The least-stable kinetic eigenvalue's real part is unchanged between the two (it matches
Fortran kinetic-DCON to 0.17% either way).
Why the nominal match is wrong
QuadGK's rtol is a bound on the converged integral; LSODE's rtol is a bound on each
step, so LSODE's effective integral accuracy is much looser than its nominal rtol.
rtol_xlmda=1e-3 (with a loose, rtol-dominated atol) is the matched-accuracy analogue of
PENTRC's step-level 1e-5.
Suggested change
Default rtol_xlmda to ~1e-3 with a loose (rtol-dominated) atol_xlmda, and/or document
that the QuadGK integral-level tolerance is not directly comparable to a PENTRC/LSODE
step-level tolerance. (Separately: the eigenvalue's imaginary/damping part is far more
tolerance/normalization sensitive — tracked separately.)
Summary
The self-consistent kinetic matrix build uses adaptive Gauss–Kronrod (QuadGK) for the
pitch/energy integrals, but the default
rtol_xlmdais set to match PENTRC's nominalvalue. QuadGK bounds the overall-integral error while PENTRC/LSODE bounds per-step
error, so matching the nominal tolerance over-refines Julia by ~3× at no accuracy benefit.
Where
KineticForces.compute_calculated_kinetic_matrices->compute_kinetic_matrices_at_psi!->integrate_pitch_gar_quadgk(QuadGK.quadgk!withatol=atol_xlmda, rtol=rtol_xlmda). The control default isrtol_xlmda=1e-5,atol_xlmda=1e-8.Measurement (matched grid: mpsi=257, mpert=27, nl=6, 8 threads)
rtol_xlmdaatol_xlmdaThe least-stable kinetic eigenvalue's real part is unchanged between the two (it matches
Fortran kinetic-DCON to 0.17% either way).
Why the nominal match is wrong
QuadGK's
rtolis a bound on the converged integral; LSODE'srtolis a bound on eachstep, so LSODE's effective integral accuracy is much looser than its nominal
rtol.rtol_xlmda=1e-3(with a loose, rtol-dominatedatol) is the matched-accuracy analogue ofPENTRC's step-level
1e-5.Suggested change
Default
rtol_xlmdato ~1e-3with a loose (rtol-dominated)atol_xlmda, and/or documentthat the QuadGK integral-level tolerance is not directly comparable to a PENTRC/LSODE
step-level tolerance. (Separately: the eigenvalue's imaginary/damping part is far more
tolerance/normalization sensitive — tracked separately.)