Skip to content

PENTRC: omega_D in the GAR branch carries a spurious major radius - #281

Open
krystophny wants to merge 4 commits into
PrincetonUniversity:developfrom
krystophny:agent/pentrc-gar-precession-length
Open

PENTRC: omega_D in the GAR branch carries a spurious major radius#281
krystophny wants to merge 4 commits into
PrincetonUniversity:developfrom
krystophny:agent/pentrc-gar-precession-length

Conversation

@krystophny

@krystophny krystophny commented Jul 25, 2026 •
edited
Loading

Copy link
Copy Markdown
Collaborator

Stacked on #280 (agent/pentrc-nested-quadrature-tolerances); please read that one first, since the ITER case cannot run at all without it.

What this changes

One factor on one line in pentrc/torque.F90, in the general-aspect-ratio bounce average:

-                    wdbar = ro*ro*bo*wdfac*wbbar*2*(2-sigma)*bspl%fsi(bspl%mx,2)
+                    wdbar = ro*bo*wdfac*wbbar*2*(2-sigma)*bspl%fsi(bspl%mx,2)

Why

The operator is correct and is not touched. With vpar = 1 - lmda*b/bo, the second bspl integrand is exactly d/dPsi of the first's action companion at fixed lmda:

d/dPsi [ sqrt(vpar) J b ]
  = -lmda b (db/dPsi) J/(2 bo sqrt(vpar)) + sqrt(vpar)[ (dJ/dPsi) b + J db/dPsi ]
  = J (db/dPsi) (1 - 1.5 lmda b/bo)/sqrt(vpar) + (dJ/dPsi) b sqrt(vpar)

which is the line at torque.F90:708. So with dl = J b dtheta and J_par = oint m v_par dl, the bounce average is the canonical precession

|omega_phi| = |(1/chrg) (dJ_par/dpsi_p) / tau_b| = 4 pi (T/chrg) |I2/I1|

the 4 pi because chi1 = twopi*psio makes the tabulated eqfun_fx(1)/chi1 and rzphi_fx(4)/chi1**2 derivatives with respect to the full poloidal flux, while the canonical momentum wants psi_p = Psi_pol/(2 pi).

The overall sign is a convention question and this PR does not settle it. The familiar form omega_phi = -(1/chrg)(dJ_par/dpsi_p)/tau_b needs P_phi = +chrg psi_p; with P_phi = -chrg psi_p the sign reverses. GPEC cannot be read off as one or the other from here: EFIT input signs are normalised away in equil/read_eq.f:666-680, and the native toroidal coordinate is mapped to the laboratory angle with phi = -helicity*(2 pi zeta + dphi) in coil/field.F:169-190, with helicity declared in coil.in rather than derived from the equilibrium. So what is established below is the magnitude, and with it the removal of exactly one positive factor ro. The patch multiplies by a positive number and changes no sign; the sign question is tracked separately.

The prefactor is not. Substituting the code's own definitions,

wbbar = ro*twopi/((2-sigma)*I1)             bhat = sqrt(2T/mass)/ro
wdbar = ro*ro*bo*wdfac*wbbar*2*(2-sigma)*I2 dhat = (T/chrg)/(bo*ro*ro)

omega_b = wbbar*bhat = 2 pi/tau_b                      <- ro cancels, correct
omega_D = wdbar*dhat = 4 pi ro wdfac (T/chrg) I2/I1    <- one ro survives

wbbar is omega_b/bhat, so it carries a factor ro that bhat removes again. Reusing it inside wdbar imports that ro a third time, and dhat cancels only the two written explicitly.

The surviving factor is a length, not a 2 pi. This matters because ro = 6.411 m and 2 pi = 6.283 are only 2% apart. I1 is in metres since J b dtheta = dl; in I2 every term carries the same J*b product divided by Psi, so I2/I1 has units 1/Wb whatever internal normalisation b and J carry — any constant rescaling of J*b cancels in the ratio. With T/chrg in volts and V/Wb = 1/s, the correct expression is already a frequency and omega_D as computed has units of m/s. No dimensionless constant can repair that.

Verification

Against an independently validated canonical precession operator (checked against an integrated guiding-centre orbit) on the same equilibrium, read through a Boozer description rather than DCON's chart. Metric is the least-squares slope of omega_D against the reference over all pitches on a surface — a per-point ratio is ill conditioned where omega_D crosses zero.

median slope range, 7 surfaces predicted
before -6.494 -6.795 to -6.087 -ro = -6.411
after -1.021 -1.063 to -0.965 -1.000

The minus signs in the predicted column are the observed relative orientation of the two operators carried across, not a derived one — see the convention caveat above. What the measurement establishes is that the ratio's magnitude goes from ro to 1.

omega_b, built from the same bspl%fsi(:,1), agrees to within [0.9964, 1.0011] both before and after. That is the control: it pins I1 to SI metres, and with it the units of I2/I1, and it shows the change is confined to the precession prefactor.

Residual 1–2% is the difference between the two equilibrium descriptions for this operator, independently measured at 0.6–5.3%.

Supporting checks:

  • chi1 is SI, not assumed to be. Integrating dpsi_pol/ds = |iota| Phi_edge/(2 pi) over the Boozer file gives psi_p(edge) = 11.8809 Wb against DCON's chi1/(2 pi) = 11.8828 Wb, agreeing to 1.6e-4. Had the 2 pi sat elsewhere in the chain the prediction would have been 2 pi ro = 40.3, which the measurement excludes.
  • Trapped domains match independently. PENTRC reports Lambda in [0.9064, 1.0773] at rho_tor = 0.251; bo times the Boozer chart's own [1/Bmax, 1/Bmin] is [0.9063, 1.0777].
  • PENTRC's own second estimator agrees with the fix. omega_d_rlar is wdhat = q T/(2 eps ro**2 chrg bo), which is dimensionally a frequency and reproduces to 8 digits (591.1770263 computed, 5.91177020E+002 written). It sits at 0.78x the canonical operator, while the GAR bounce average sat at 8.26x it. The two are not the same quantity, so only the order of magnitude is meant — but the dimensionally correct estimator inside PENTRC lands where the reference does.

Physical consequence on the ITER case

The superbanana resonance sits at x = -omega_E/omega_D. With omega_E ~ -6.4e3 rad/s, the corrected omega_D ~ 2.3e2 puts it at x ~ 28, out in the Maxwellian tail, where an inflated omega_D had put it at x ~ 4.3.

Same deck, byte-identical pentrc.in, only the binary differing, over 0.2 <= rho_tor <= 0.8:

channel before after
TGAR trapped -11.762 N m -0.493 N m
PGAR passing -0.544 N m -0.917 N m

The trapped channel drops by 24x and the balance becomes passing-dominated, which is what the drift-kinetic codes we compare against report for this case.

Gates

Adjacent, reported and not changed

The reduced-large-aspect-ratio branch has its own dimensional oddities — wbbar there carries sqrt(bo), and unlike the GAR branch bhat is not divided by ro (torque.F90:457-461). Our deck exercises TGAR/PGAR, so I have not touched RLAR; flagging it in case it is of interest.

Charge-verified TC24 production revalidation (2026-07-26)

The corrected single-D campaign was rerun from Xingting's archive (SHA-256 655ffbf16e35f6218d1b40f58b2e2c4c16fbbcd8b1fce58c3adad515a7f17585) at this exact head, 2e031ae675205ec9d69bb0298855155737f72579. It uses ni=ne=supplied density, physical Ti/Te, mi=2, zi=1, and the paired idcon_equil.out psi_pol -> s_tor map.

Over the exact 0.2 <= rho_tor <= 0.8 window, the executable-native endpoints are:

deck trapped [N m] passing [N m] total [N m]
L0 +0.0508510 -0.00330114 +0.0475499
L5 -0.469088 -0.948124 -1.417212

The channel sum reconstructs the native cumulative torque. These numbers supersede the earlier old-profile physical-consequence table as TC24 comparison evidence. No sign was transformed: the audit retains native signs, while publication remains magnitude-only until the exact runtime helicity and manufactured nonzero complex-mode torque/power gate close.

@krystophny
PENTRC: fix OpenMP data race on the energy contour flag
cb5e7c6 
energy_imaxis selects the integration contour inside xintgrnd -- cx = i*x on the
imaginary-axis leg against cx = x + i*ximag on the real one -- and xintgrl_lsode
flips it between its two legs (energy.f90:190, :223). It is per-integration state
exactly like energy_wn, energy_wd, energy_leff and the rest, and it was the only
one of them missing from the threadprivate list.

Shared, one thread's imaginary leg silently switches another thread's integrand
from a decaying exponential to an oscillation, so values from deep in the
Maxwellian tail come back four orders of magnitude too large and the pitch
integrand stops being a function of lambda: neighbouring lambda separated by 1e-8
differ by 130% on average. lambdaintgrl_lsode then shrinks its step without bound
and aborts with 'too many steps in lambda required'.

That is why the failure resisted every numerical explanation -- tolerances, xmax,
ximag, nl, compiler, kinetic profiles: it does not depend on the numerics, it
depends on the thread count. Decks that set the deprecated and ignored
openmp_threads = 1 silently ran on every available core.

Gate: the ITER case that previously aborted now runs unmodified and gives
bit-identical pentrc_tgar_n3.out and pentrc_pgar_n3.out at 1 and 16 threads, all
29 columns, including the pitch-resolved netCDF.

Note for reviewers: recorded integrations (op_record=.true.) additionally mutate
the shared energy_record object without synchronization. Production recording is
serial so this does not affect the fix above, and it is left alone here.
@krystophny
PENTRC: give the energy integration its own tolerances
f2cd940 
pentrc_interface aliased both xatol/xrtol and lambdaatol/lambdartol to the single
atol_xlmda/rtol_xlmda pair. The pitch integral's integrand IS the energy integral,
so this is a nested adaptive quadrature, and with one tolerance for both the outer
integrator chases the inner one's quadrature error -- which is not a smooth
function of lambda. LSODE then shrinks its step without bound.

New pent_control entries atol_x/rtol_x carry the energy tolerances. Their default
of -1 means 'derive from the pitch tolerances', as nested_tolerance_margin = 1e-2
times them.

This DOES change numerical behaviour for every deck that does not set the new
controls: the effective energy tolerances go from the formerly aliased 1e-6/1e-3
to 1e-8/1e-5, i.e. tighter. Decks can restore the old values explicitly by setting
atol_x = 1e-6, rtol_x = 1e-3. Note that energy.f90's own module defaults are
tighter still, 1e-12/1e-9; those were overridden by the aliasing before this change
and remain overridden by the derived values after it, which is a separate question
about what the right energy tolerance actually is.

This is NOT what unblocked the previously-aborting ITER case; the OpenMP race on
energy_imaxis was. It is a real defect in its own right and is kept separate.

Dropped from an earlier version of this change: an exclusion of both endpoints of
the RLAR trapped pitch space. lmdamin = max(1/(1+epsr), bo/bmax) and
lmdamax = min(1/(1-epsr), bo/bmin) are degenerate only when the bo/b term wins;
when the cylindrical epsr bound wins the endpoint has two regular bounce points,
and excluding it unconditionally would shrink the physical interval and discard
nonsingular contributions.
@krystophny
PENTRC: omega_D in the GAR branch carried a spurious major radius
1ad89ff 
wdbar reused wbbar, which is omega_b/bhat and therefore carries one factor of ro
that bhat takes back out again. dhat cancels only the ro**2 written explicitly, so
the third ro survived and omega_D came out with units of m/s rather than 1/s.

The operator itself is correct and is not touched: bspl%fsi(:,2) is exactly
d/dPsi of the parallel action integrand at fixed lmda -- the 1 - 1.5*lmda*b/bo
coefficient follows from sqrt(vpar) = vpar/sqrt(vpar) in the product rule -- so the
bounce average is the canonical precession
-(1/chrg)(dJ_par/dpsi_p)/tau_b = -4 pi (T/chrg) I2/I1, the 4 pi because
chi1 = twopi*psio makes the tabulated derivatives d/dPsi_pol while the canonical
momentum wants d/dpsi_p. Only the prefactor changes: ro, not ro*ro.

Verified against an independently validated canonical operator on the same
equilibrium, as the least-squares slope over all pitches on a surface (a per-point
ratio is ill conditioned at the omega_D zero crossing): -6.494 before against the
predicted -ro = -6.411, and -1.021 after against -1, over 7 surfaces. omega_b,
built from the same bspl(:,1), agrees to better than 0.4% before and after and is
untouched.

The dimensional argument is what makes this a length rather than a missing 2 pi,
which is only 2% away numerically: I2/I1 has units 1/Wb whatever normalisation b
and J carry, since any constant rescaling of J*b cancels in the ratio, so
(T/chrg) I2/I1 is already a frequency. Confirmed that no hidden normalization
supplies the extra radius: rzphi stores physical R**2 and a physical-coordinate
Jacobian, eqfun stores B in tesla, ro is in metres.

A second witness is internal: omega_d_rlar in the same output file is
wdhat = q T/(2 eps ro**2 chrg bo), dimensionally a frequency, reproducing to 8
digits, and it sits at 0.78x the canonical operator while the GAR bounce average
sat at 8.26x it.

SCOPE OF THE SIGN CLAIM. Only the magnitude is established here. The overall sign
of the canonical form depends on whether the thin-orbit canonical momentum is
+chrg*psi_p or -chrg*psi_p, and GPEC cannot settle that from its stored radial
coordinate: EFIT input signs are normalized away in equil/read_eq.f and the native
toroidal coordinate reaches the laboratory frame only through
phi_CCW = -helicity*(2 pi zeta + dphi) in coil/field.F. This change removes one
positive factor of ro and nothing else; no sign is altered, and the residual
orientation question is independent of it.

On the ITER case this moves the trapped superbanana resonance from x ~ 4.3 to
x ~ 28, so TGAR torque falls from -11.8 to -0.5 N m and the balance becomes
passing-dominated. Output stays bit-identical at 1 and 16 threads.
@krystophny
krystophny force-pushed the agent/pentrc-gar-precession-length branch from 71e57fc to 1ad89ff Compare July 25, 2026 12:46
@krystophny

krystophny commented Jul 25, 2026

Copy link
Copy Markdown
Collaborator Author

Rebased onto current develop (was 12 commits behind) and narrowed after an independent adversarial review. The rebased build reproduces the pre-rebase output byte for byte, so nothing numerical moved.

Scope of the sign claim, corrected. The review is right that the PR asserted more than GPEC can support. The relation

omega_phi = -(1/q) (dJ_par/dpsi_p) / (dJ_par/dE)

follows only if the thin-orbit canonical momentum contains +q*psi_p; with -q*psi_p the sign reverses. GPEC cannot settle which it is from its stored radial coordinate, because EFIT input signs are normalized away in equil/read_eq.f and the native toroidal coordinate reaches the laboratory frame only through phi_CCW = -helicity*(2*pi*zeta + dphi) in coil/field.F.

So this PR establishes the magnitude and nothing about orientation. It removes one positive factor of ro. No sign is changed anywhere, and the residual orientation question is independent of it. The commit message now says this explicitly.

What the review independently confirmed, having been asked specifically to refute it:

  • The 1 - 1.5*lmda*b/bo coefficient at torque.F90:707-710 is exact — d/dPsi(J*b*sqrt(v)) reproduces the source line.
  • rzphi_fx(4)/chi1**2 is consistent with rather than evidence against the full-poloidal-flux reading: if rzphi_f(4)/chi1 is the Jacobian w.r.t. full flux, its derivative necessarily carries chi1**2.
  • No hidden normalization supplies the extra radius: rzphi stores physical R**2 and a physical-coordinate Jacobian (equil/direct.f:157-175), eqfun stores B in tesla (equil/direct.f:228-249), ro is in metres, bo = abs(sq%f(1))/(twopi*ro). So I1 is in metres, I2/I1 in 1/Wb, and the old ro*ro leaves one metre after dhat.
  • No consumer applies a compensating radius: the corrected fbnce%fs(:,2) flows into pitch.f90:343-369 and the resonance denominator and into the pitch netCDF at pitch.f90:808-824. The independent RLAR estimate written from torque.F90:1914-1916 is intentionally unaffected.
  • The RLAR dimensional warning is valid, and leaving RLAR alone was appropriate for a one-line GAR fix.

Also dropped from the stack below this one: the degenerate-endpoint change in the RLAR branch, which the review showed is unsound when the cylindrical epsr bound sets the interval end. See #280.

@krystophny
PENTRC: state the precession identity as a magnitude in the comment too
2e031ae 
The PR body was corrected to claim only a magnitude, since the familiar minus
sign needs P_phi = +chrg*psi_p and GPEC does not establish that convention.
The source comment still asserted the signed identity, which is the exact claim
the body disclaims. No executable change.
@krystophny

krystophny commented Jul 26, 2026

Copy link
Copy Markdown
Collaborator Author

Final stacked production revalidation (after #280, this PR, and the independent periodic-spline #284 fix): all four phiI000/phiI010, L0/L5 lanes complete. The previously absent lower-row curves were not zero—they were missing because phiI010 L5 hit a nonperiodic parallel-velocity spline seam and then non-finite omega_b/omega_D; #284 fixes that without changing any frequency or torque sign.

The clean production endpoints are:

deck L0 total [N m] L5 total [N m]
phiI000 +0.0475499 -1.417212
phiI010 +3.9190 -4.3130

These are executable-native cumulative GAR totals on the GPEC field. Publication remains magnitude-only until the manufactured complex-mode helicity gate closes. One-thread and 16-thread outputs for the repaired failing lane are bit-for-bit identical.

@krystophny krystophny mentioned this pull request Jul 30, 2026
Closed
@logan-nc

logan-nc commented Aug 1, 2026

Copy link
Copy Markdown
Contributor

@claude review this pr

@claude

claude Bot commented Aug 1, 2026 •
edited
Loading

Copy link
Copy Markdown

Claude finished @logan-nc's task in 2m 26s —— View job

Reviewing PR #281

Summary

This PR's own contribution (commits 1ad89ff/2e031ae) is a single-line, well-contained fix in pentrc/torque.F90:753:

-wdbar = ro*ro*bo*wdfac*wbbar*2*(2-sigma)*bspl%fsi(bspl%mx,2)
+wdbar = ro*bo*wdfac*wbbar*2*(2-sigma)*bspl%fsi(bspl%mx,2)

(pentrc/energy.f90 and pentrc/pentrc_interface.f90 only show up in the diff because this branch is stacked on unmerged #280 — no new changes to those files here.)

Verification performed

  • Algebra: re-derived omega_D = wdbar*dhat by hand from the code's own definitions of wbbar, bhat, dhat. Confirms the claim exactly: wbbar = ro*twopi/((2-sigma)*I1) carries one factor of ro; dhat's explicit ro**2 only cancels two of the three ro's that accumulate when wbbar is reused inside the old wdbar = ro*ro*bo*.... Before: omega_D = 4*pi*ro*wdfac*(T/chrg)*I2/I1 (units m/s). After: omega_D = 4*pi*wdfac*(T/chrg)*I2/I1 (units 1/s). The fix is dimensionally necessary, not just numerically plausible.
  • Blast radius: wdbar is a local variable used nowhere else in the codebase (only one hit for wdbar repo-wide), and fbnce%fs(:,2) flows into pitch.f90:344-369 (wd = real(flmabda_f(2))) with no compensating factor applied downstream — consistent with the PR's claim that no consumer was silently correcting for the extra ro.
  • RLAR branch untouched, as claimed: confirmed bhat at torque.F90:450 (RLAR branch) is not divided by ro, unlike the GAR branch's bhat at line 754 — so the two branches are genuinely independent and leaving RLAR alone is correct, not an oversight.
  • Could not compile/run the Fortran build in this environment (no NetCDF/toolchain configured here), so I relied on static/algebraic verification and cross-referencing consumers rather than reproducing the numerical tables in the PR description.

Assessment

  • The fix is minimal, dimensionally justified, and self-consistent with the surrounding code. I did not find any other latent ro factor or missed consumer.
  • The PR is explicit and appropriately conservative about not resolving the overall sign convention (P_phi = ±chrg*psi_p), correctly scoping this change to magnitude only. That's the right call given equil/read_eq.f normalizes away EFIT sign and coil/field.F's helicity is a declared, not derived, quantity — nothing in this diff needs to (or does) touch sign.
  • The inline comment at torque.F90:741-752 is long (12 lines) for a one-line change, but given how subtle and easy to reintroduce this bug is (a 2% numerical difference between ro and 2*pi in the ITER case masked it), I think the extra derivation is warranted here rather than a nitpick.
  • No .f-fixed-format files were touched (this is .F90, free-form), so the 72-column rule doesn't apply; line lengths in the diff are reasonable.
  • Nothing to request changed. The omega_b control check described in the PR body (agreement to [0.9964, 1.0011] before/after) is a good sanity anchor, since it pins I1's units and shows the change is confined to the precession prefactor as intended.

Branch: agent/pentrc-gar-precession-length

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

Reviewers

No reviews

Assignees

No one assigned

Labels

None yet

Projects

None yet

Milestone

No milestone

Development

Successfully merging this pull request may close these issues.

2 participants

@krystophny @logan-nc