KF/PE - BUGFIX! - Self-consistent xi_s for the NTV kernel input; static scheduling for the kinetic kernel loop - #425
Conversation
…iguous case); PR #425 Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk
2e168ab to
d5eec3c
Compare
7c47c07 to
3b37787
Compare
530ade4 to
b3d7ab9
Compare
3b37787 to
861c998
Compare
…eq.f kin_flag analog) The Clebsch xi_alpha fed to the NTV kernel was re-solved from the IDEAL A,B,C even in kinetic runs, inconsistent with the self-consistent kinetic tangential response exactly where resonant layers live. Mirror Fortran: scale the stored xi_s by the singfac regularization factor and invert nothing (analog of the #407 fix, adapted to the FourFitVars API; kinetic-ness detected via populated kwmats pending proper flag threading). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk
The per-surface kernel loop indexed per-thread scratch (thread_intrs and the four matrix buffers) by threadid() under DEFAULT dynamic scheduling -- the only threaded loop in the codebase missing :static. Task migration at yield points lets two tasks share a buffer and corrupt matrix rows, making every calculated-kinetic run nondeterministic at rare-migration frequency (observed: three runs of an identical DIII-D deck gave et[1] = 1.0053, 1.1355, 0.9895). All other threadid-indexed loops already use :static per the GPEC convention. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk
861c998 to
225ac39
Compare
|
Closing without merging — both halves are superseded, and one of my claims here needs correcting. 1. The Correction to this PR's description: it claimed "every historical calculated-kinetic run carries this latent race." That is wrong as written. The race is real and was real on this stack, but develop has been fixed since 2. The ξ_s self-consistency fix is owned by #407. Same file, same function, same physics (kinetic runs must regularize the stored self-consistent ξ_s rather than re-solve it from the ideal matrices, mirroring Fortran Worth flagging for #407 and #414 reviewers: #414 forces The measurement this PR produced still stands and is recorded in #414: with regularization off, the NTV torque and the Euler–Lagrange solution's own dissipation agree to 0.15% (0.1324 vs 0.1322 N·m) on the DIII-D kinetic-calculated case. Branch |
Stacked on #422 (
feature/kinetic-resonance-knots). The overnight resolution of the PE-vs-NTV torque discrepancy (issue #423's measured motivation; RESULTS.md §34–§35).The bug
In kinetic runs, the Clebsch ξ^α fed to the NTV kernel (via
dbob_m/divx_m) was re-solved from the ideal A, B, C matrices — inconsistent with the self-consistent kinetic tangential response exactly where resonant layers live. Fortran (gpeq.f,kin_flag) instead scales the stored self-consistent ξ_s by the singfac regularization factor and inverts nothing. This is the same physics flaw #407 fixes on the MatrixSplines branch (there: a silently broken Cholesky of the non-Hermitian kinetic A; here: a valid Cholesky of the wrong — ideal — A).The consequence (measured, DIII-D kinetic-calculated full chain)
The perturbative NTV evaluation was driven by a tangential displacement that ignores the kinetic layer response, producing a phantom 4× torque concentrated in a ~2e-3-wide resonant core at the pedestal ω_E crossing (fgar: 0.7241 N·m, of which +0.622 in ψ∈[0.90, 0.95]; the self-consistent EL solution collects 0.046 there — and the EL side was verified faithful at every level: boundary ≡ volume torque to 7 digits, FKG assembly and kinetic RHS sign-by-sign Fortran-conformant per fortran-physics-reviewer, matrix-spline resolution and integrator step size both exonerated by direct null experiments).
With the fix, the phantom collapses (+0.622 → +0.033) and the two independent torque calculations agree band-by-band (resonant layer: 0.0461 vs 0.0455). Domain-matched totals (ΔIm over [ψ_c, psilim] — fgar's exact domain; the sub-ψ_c Im(u1†u2) offset is non-dissipative initial-condition content, ±0.05 with opposite signs across cases):
Where the comparison is regularization-unambiguous (rot=0.2: resonances away from rationals), the two independent torque calculations now agree to 0.3%. The rot=1.0 residual sits inside fgar's own reg_spot sensitivity band (its answer moves 18% between reg 0.05 and 0.01 — the acknowledged near-rational input ambiguity, present in Fortran too); PE is reg-independent by construction.
Second fix in this PR: determinism
CalculatedKineticMatrices.jl's per-surface kernel loop was the only threaded loop in the codebase usingThreads.@threadswithout:staticwhile indexing per-thread scratch bythreadid()— dynamic-scheduling task migration corrupts kinetic matrix rows at rare-migration frequency. Observed live: three runs of one identical deck gave et[1] = 1.0053 / 1.1355 / 0.9895. With:static, identical decks reproduce bit-identically (verified), and et returns exactly to the pre-fix value — every historical calculated-kinetic result carries this latent race.Effects on outputs: kinetic-run regularized quantities (
xi_clebsch_alpha,*_regfields) and every NTV torque driven by them change; ideal runs andreg_spot = 0are bit-unaffected. Kinetic harness cases will move — new baselines with this physics justification.Review items
kwmatssplines (sentinel = 5 knots); reviewers may prefer threading an explicit flag fromkinetic_factor > 0— noted as the clean alternative.gpout_dwψ-profile port, issue Port gpout_dw / gpout_dw_matrix (psi-resolved torque response profiles and matrices) as the grid-logic diagnostic #423) is where the "which torque is which" delineation (PerturbedEquilibrium - DOCS - Delineate the boundary-response torque from the KineticForces NTV torque #424) completes.🤖 Generated with Claude Code
https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk