PENTRC: restore pitch-consistent clar normalizations (omega_b /sqrt(bo), omega_D /bo)

[logan-nc]

Fixes #287. Companion to #281, same root cause found via the git archaeology done for that review.

What this changes

Two factors in the clar bounce-average branch of pentrc/torque.F90:

-                    wbbar = pi*SQRT(2*epsr*lmda*bo)/(4*q*ro*ellipk(kappa**2))
+                    wbbar = pi*SQRT(2*epsr*lmda)/(4*q*ro*ellipk(kappa**2))
                     wdbar = (2*q*lmda*(ellipe(kappa**2)/ellipk(kappa**2)-0.5)&
-                        /(ro**2*epsr))*wdfac
+                        /(bo*ro**2*epsr))*wdfac

plus a comment recording why.

Why

In PENT (pent/profile.f, 2013-2015) the pitch variable was dimensional (vpar = 1 - lmda*b, lmda in 1/Tesla) and these formulas were consistent with it: lmda*bo under the square root was dimensionless, and wdbar*dhat = (T/chrg)*2*q*lmda*G/(ro**2*eps) was a frequency because lmda carried the 1/bo.

The PENT->PENTRC rewrite (da0e6067, June 2015) nondimensionalized the pitch (Lambda = lmda_pent*bo, vpar = 1 - (lmda/bo)*b) but kept the clar text verbatim. Substituting lmda_pent = Lambda/bo:

• kappa maps exactly (PENT's kk = (1/bo - lmda(1-eps))/(2*eps*lmda) is identical to PENTRC's kk = (1 - lmda(1-eps))/(2*eps*lmda)), so trapped-space geometry needs no change;
• wbbar retained a spurious sqrt(bo) and wdbar a spurious bo.

So clar has reported omega_b inflated by sqrt(bo) (~1.4-2.3x) and omega_D inflated by bo (~2-5.3x) since 2015. This PR restores exactly the forms PENT evaluated. This is the same failure mode as the GAR ro in #281 — normalization conventions changed under formulas that were not updated — and, like #281, it multiplies by positive numbers and changes no sign.

Scope

• Affects only the clar method (fbnce channels 1-3: omega_b, omega_D, and the wbbar-scaled action).
• The rlar method's flux-surface estimators (wbhat/wdhat, built from wtran/wgyro) are dimensionally sound and untouched.
• GAR methods are handled by PENTRC: omega_D in the GAR branch carries a spurious major radius #281; independent lines, applies cleanly in either order.

Verification status

Analytic + historical (the PENT mapping above); compiles cleanly with gfortran. It has not been benchmarked against an independent reference operator the way #281 was — a re-run of #281's slope-vs-reference measurement on the clar branch should read sqrt(bo) and bo before this patch and ~1 after, and would be a welcome check before merging.

🤖 Generated with Claude Code

https://claude.ai/code/session_01Rewp8poxoQyorK48SwxEUv

[logan-nc]

Ran the verification exercise recommended in the PR body: clar omega_b and omega_D vs pitch compared against the GAR bounce average on the same surfaces, before vs after this patch. Reference is tgar from a build with #281 applied (GAR omega_b was always correct; GAR omega_D is trusted after #281). Testbed: the in-repo DIIID_kinetic_example deck (shaped H-mode, DCON bo = 1.983 T), DCON -> GPEC -> PENTRC with xlmda_out recording on six surfaces, fitting least-squares slopes of clar against tgar over the recorded pitch grid (same-sign, away-from-zero mask for omega_D, which crosses zero).

Mechanical check — pointwise ratio of clar before/after this patch, all pitches, all surfaces:

quantity	measured ratio	predicted
omega_b	1.40816 +/- 0.00003	sqrt(bo) = 1.40824
omega_D	1.98321 +/- 0.00003	bo = 1.98321

The patch multiplies by exactly the claimed factors and nothing else (the omega_b ratio squared reproduces the omega_D ratio to 5 digits, which is how bo above was extracted from the code itself).

Physics check — slope of clar against tgar reference:

psi_n	eps	omega_b before	omega_b after	omega_D before	omega_D after
0.15	0.11	1.433	1.018	2.22	1.12
0.25	0.15	1.218	0.910	3.80	1.61
0.35	0.18	1.157	0.818	1.81	0.89
0.45	0.20	1.184	0.851	1.72	0.88
0.55	0.23	1.148	0.812	1.86	0.91
0.65	0.26	1.052	0.748	1.69	0.83

Predicted: before = 1.408 / 1.983, after = 1.0. On the most LAR-like surface (eps = 0.11) the post-fix slopes are 1.018 and 1.12 where the pre-fix code sat at 1.43 (= sqrt(bo)) and 2.22 (= bo x 1.12). The drift away from 1 at larger eps is the circular-large-aspect-ratio approximation error on a shaped, up-down-asymmetric equilibrium — it is identical before and after (the before/after ratio is sqrt(bo), bo at every surface), so it does not confound the measurement. The psi=0.25 omega_D fit is the least well conditioned (zero-crossing location differs most from GAR there).

This is the analogue of #281's slope test with the same outcome pattern: slope magnitude moves from the predicted spurious factor to ~1, with the residual explained by the approximation the method is defined by.

🤖 Generated with Claude Code

https://claude.ai/code/session_01Rewp8poxoQyorK48SwxEUv

[logan-nc]

Correction / downgrade of the previous comment's "physics check".

Red-teaming the verification exposed a flaw in it: the clar-vs-tgar slope comparison uses one in-repo method as the reference for the other, and on the DIII-D-like case bo/ro = 1.98/1.74 = 1.14, so the unpatched clar and tgar also mutually agree within the fit scatter (measured: original-vs-original wd slope median ~1.04). That table therefore cannot distinguish "this patch is right" from "the original pair was mutually consistent" — the two hypotheses only separate on a machine where bo/ro is far from 1.

I ran that discriminating case: the in-repo a10_kinetic_example (circular, bo = 0.9997 T, ro = 2.0046 m, so bo/ro = 0.50). Result: it does not discriminate either, for an instructive reason. The within-build clar/tgar omega_D slope lands between the two predictions (dev build 0.67, fixed build 1.35 at eps = 0.034, drifting to 0.21/0.43 by eps = 0.074, identically in both builds), and the omega_b control reads 0.85 flat where both hypotheses predict 1.0. Conclusion: the clar model's intrinsic error is O(1) for omega_D even on a circular equilibrium — its G(kappa) = E/K - 1/2 carries no magnetic-shear term, and the radial drift of the ratio is consistent with the missing shear contribution — plus a ~15% systematic in omega_b (epsilon-definition level). Cross-method comparison inside PENTRC can never resolve a factor-bo question against model error of that size.

What still stands, and what the case for this patch actually rests on:

1. The mechanical check from the previous comment is unaffected: the patch multiplies clar omega_b by exactly 1/sqrt(bo) and omega_D by exactly 1/bo (verified pointwise to 5-6 digits on both equilibria) and does nothing else.
2. The units chain, anchored only by quantities no one disputes: omega_b (unchanged by every PR in this family, and cross-method consistent to 15%) pins I1 = oint J b dtheta / sqrt(vpar) in metres; the diamagnetic frequencies (inputs.f90:267-268, validated physics) pin chi1 in webers. With the pitch dimensionless (vpar = 1 - (lmda/bo)*b, torque.F90:468), the original clar omega_D has units of V/m^2, not 1/s. No dimensionless constant can repair that; 1/bo is the unique one-factor fix consistent with the PENT-era formulas.
3. The archaeology: PENT (2013-2015) evaluated these exact formulas with a dimensional pitch (vpar = 1 - lmda*b), where they are consistent. The alternative reading — that the 2015 rewrite deliberately and silently applied two dimensionally distinct corrections (ro to GAR, bo to clar) while simultaneously redefining the pitch convention, with no documentation — requires the corrections to be wrong-by-construction on unit grounds (point 2).

So the honest status is: this patch is established by construction + units + provenance, not by an independent numerical benchmark, and the empirical gate proposed in the PR body — pointing the external reference operator of #281 at the clar branch (predicted slopes sqrt(bo) and bo before, ~1 after, with O(1) model scatter now expected from the shear-free G(kappa)) — remains the right merge criterion, with the caveat that the comparison should be restricted to low-shear surfaces or use a shear-corrected reference to be sharp.

🤖 Generated with Claude Code

https://claude.ai/code/session_01Rewp8poxoQyorK48SwxEUv

[logan-nc]

Documentary closure: checked against the PENT paper itself (Logan, Park, Kim, Wang & Berkery, Phys. Plasmas 20, 122507 (2013) — the paper cited at torque.F90:454 as the source of the bar normalization).

• Above Eq. (19): the pitch coordinate is Lambda = mu*B0/E — dimensionless. The clar lines predate this convention (PENT's internal lmda = mu/E was dimensional); this PR's diagnosis is exactly the paper-vs-code mismatch.
• Text after Eq. (20): "both wb_bar = wb*R0/sqrt(2xT/M) and dJ_bar^2 = dJ^2/2xTMR0^2 are unit-less quantities" — the code's bhat scheme, with bar quantities dimensionless by construction. Figs. 3 and 5 plot precession as omega_D/(E/e) * (R0^2*B0) — i.e. dhat = (T/Ze)/(bo*ro^2) per unit x, verbatim.
• Eq. (22) gives the RLAR precession as (q^3*Lambda/eps)*x*[E(k)/K(k) - 1/2] * omega_t^2/omega_g; expanding omega_t^2/omega_g yields x*(T/e)*2*q*Lambda*(E/K - 1/2)/(eps*B0*R0^2) — exactly the corrected line in this PR, with the B0 entering via the gyrofrequency. (The Lambda-averaged limit is the code's untouched wdhat estimator, which was always dimensionally sound.)
• Historical note: the paper's verification equilibria (Fig. 1) are Solov'ev with B0 = 1 T, R0 = 1 m — the unique point where both this factor (bo) and PENTRC: omega_D in the GAR branch carries a spurious major radius #281's (ro) equal 1 and are numerically invisible. Re-validation against the paper's own benchmarks could never have caught either.

With this, the paper, the 2013 PENT code, the 2017 companion paper (Eq. (83)), the units chain, and the exact mechanical verification all agree; the 2015-onward lines are isolated as the sole outliers.

🤖 Generated with Claude Code

https://claude.ai/code/session_01Rewp8poxoQyorK48SwxEUv