KineticForces - BUG FIX - Coulomb logarithm uses natural log, not log10 - #337
Conversation
|
@claude review this |
|
Reviewed the diff. Summary of findings: Physics correctness — verified. Re-deriving the NRL formulary electron-ion Coulomb log λ_ei = 23 − ln(nₑ^½Tₑ^−³ᐟ²) (nₑ in cm⁻³, Tₑ in eV) and renormalizing to GPEC units (nₑ/10²⁰ m⁻³, Tₑ/1 keV) does reproduce the 17.3/−0.5/+1.5 coefficients under natural log (17.24, matching the PR's math), and the old Scope/isolation — good. Grepped the repo and this is the only place this Coulomb-log formula appears, so there's no other Regression harness — used correctly per the repo policy; the reported shifts in Comment — three lines; could be tightened to one per the "concise, one line where possible" guidance in CLAUDE.md, but it's documenting a physics formula (encouraged) rather than session narration, so this is a minor nit, not a blocker. Test coverage — this is the real gap. There's no unit test that pins down Minor/optional: the NRL formulary's 23−ln(...) form (used here) is the low/mid-temperature regime; for Tₑ ≳ 10·Z² eV the formulary switches to λ = 24 − ln(nₑ^½/Tₑ). At the 20 keV core cited in the PR description that threshold is far exceeded, so the regime choice itself may be worth a follow-up look — but that's a pre-existing modeling choice this PR doesn't touch, not something to hold up this fix on. |
The lnLambda formula 17.3 - 0.5*log(ne/1e20) + 1.5*log(Te/1keV) is the NRL Plasma Formulary electron-ion Coulomb logarithm (lambda_ei = 23 - ln(ne^1/2 * Te^-3/2)) re-normalized to (ne/1e20 m^-3, Te/1keV); its 17.3/-0.5/+1.5 coefficients are calibrated for natural log. The Julia port used log10, so lnLambda -- and the collision frequencies nue/nui that scale linearly with it -- ran ~13% low at a 20 keV core (19.2 vs the correct 21.7). Restore natural log, matching both Fortran PENTRC inputs.f90:238 and the formulary. Regression (diiid_n1): NTV FGAR and dW FGAR shift within tolerance (~0.06-0.2%; small only because that case sits near the ne=1e20/Te=1keV reference point where log10 ~= ln). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01SsLUydP2gJbE1tGoaHDiS1
76370a6 to
cd0bba6
Compare
Two conflicts, both from develop moving under this branch. The Coulomb logarithm was corrected to a natural log on both sides independently: develop in #337, this branch in its own refactor. They are numerically identical, so the helper is kept and develop's annotation already sits on its definition. The kinetic-matrix block conflicted because develop hoisted it and this branch still carried the older copy, which git left duplicated. One copy is kept, in develop's position, with this branch's `species=kf_species` argument carried into it. The branch's guard on `ctrl.ode_flag` is dropped because develop deprecated and removed that flag; the surviving guard is `singfac_min > 0`. The module loads against the merged tree.
Summary
The Coulomb logarithm in the kinetic-profile collisionality (
src/Equilibrium/KineticProfiles.jl) was computed withlog10where it must use natural log:The
17.3 / -0.5 / +1.5coefficients are calibrated for natural log — they are the NRL Plasma Formulary electron-ion Coulomb logarithmλ_ei = 23 − ln(nₑ^{1/2} · Tₑ^{−3/2}) (nₑ in cm⁻³, Tₑ in eV)
re-normalized to GPEC units (nₑ / 10²⁰ m⁻³, Tₑ / 1 keV = 1.602e-16 J): the constant becomes 17.24 ≈ 17.3, and the +1.5 coefficient on ln(Tₑ) is the signature of the Tₑ^{−3/2} regime (impossible with log₁₀). This also matches the Fortran reference
pentrc/inputs.f90:238, which useslog()(natural log in Fortran).Impact
lnΛ, and the collision frequencies
nue/nuithat scale linearly with it, ran ~13% low at a 20 keV core (19.2 vs the correct 21.7). The error vanishes near the reference point (nₑ=10²⁰, Tₑ=1 keV) and grows with temperature.Verification
diiid_n1):NTV torque FGARandNTV kinetic energy dW FGARshift +0.06% / +0.22% — within tolerance, correct direction, small because that case sits near the log reference point. No other quantity is attributable to the change (the touched code is the collision frequency; the delta-prime/resonant-field scatter is pre-existing threaded-reduction non-determinism).🤖 Generated with Claude Code
https://claude.ai/code/session_01SsLUydP2gJbE1tGoaHDiS1