[krystophny]

Two commits:

1. GPEC - BUG FIX - Form the pmodb perpendicular displacement from xi^psi and xi^s. gpout_pmodb formed ξ⊥ = ξ − ξ∥ with the Fourier-truncated 1/(m−nq) closure ξ∥, which left large residues in Bdivxprp and Bkxprp at resonant and window-edge harmonics. Following the review, it now builds them directly from (ξ^ψ, ξ^s) with the eqfun(2:3) metric combinations:
   Jξ⊥^θ = χ′²[−f₂ξ^ψ + f₃ξ^s/χ′], Jξ⊥^ζ = q·Jξ⊥^θ − Jξ^s/χ′. The divergence and κ·ξ⊥ columns both use this, with no ξ∥ and no 1/(m−nq). This is algebraically the smat…zmat form. The output format is unchanged.
2. PENTRC - BUG FIX - Read the perpendicular pmodb drive and Jacobian pair. read_pmodb keeps the supplied divergence and forms the drive as −(Bdivxprp+Bkxprp), instead of −(lagb+Bkxprp) in lagb's ξ^ζ=0 gauge.

Per-m identity. R = Σ_ψ|Bdivxprp+lagb+Bkxprp|² / Σ_ψ|lagb|², develop → this PR:

m	solovev_ideal reg=t	ideal reg=f	solovev_kinetic reg=t	kinetic reg=f
−4	87.8 → 0.0015	22.6 → 0.117	112 → 0.0005	144 → 0.0096
−3	1.18 → 0.0022	0.0007 → 0.076	0.67 → 0.0018	0.46 → 0.0097
−2	0.015 → 0.0011	0.026 → 0.042	0.013 → 0.0019	0.018 → 0.0078
−1	0.063 → 0.046	0.018 → 0.018	0.0024 → 0.0015	0.0074 → 0.0070
0	1.5e−4 → 1.7e−5	0.013 → 0.0057	0.0011 → 5.2e−5	8.8e−4 → 4.2e−4
1	0.94 → 0.0010	0.012 → 5.5e−4	5.93 → 1.5e−4	0.0027 → 1.3e−4
2	175 → 3.9e−4	0.146 → 5.0e−4	177 → 2.1e−4	0.199 → 5.3e−4
3	49.1 → 2.8e−4	101 → 1.4e−5	73.5 → 7.1e−5	16.5 → 2.8e−5

The remaining residual is the gauge term (ξ^s/(χ′q) − ξ⊥^θ)·∂_θB carried by lagb. Subtracting it brings the ideal reg=f rows at m=−4…−2 to R ≤ 1.5e−6. One exception is ideal reg=t m=−1 (0.035 after the gauge term; 0.063 on develop). It comes from a discrete force-balance term in the interpolated B and metric jets, peaking near the axis at ψ≈0.037. It is not caused by this change. An independent B∇·ξ⊥ reconstruction from (ξ^ψ, ξ^s) with analytic θ-derivatives matches the new Bdivxprp to ≤0.7% per m.

PENTRC drive (actual reader, ∫dψ Σ_m|X_m|², this PR relative to develop):

case	divx	dB/B
ideal reg=t	−1.2%	−1.3%
ideal reg=f	−2.8%	−4.3%
kinetic reg=t	−1.3%	−1.6%
kinetic reg=f	−1.5%	−2.3%

Before this change, the same comparison gave ×3.1 and ×7.2 on solovev_ideal reg=t, and up to ×69 and ×228 with reg=f.

ψ=0.4 check (Solov'ev ideal, prepared harmonic inputs), Im(Tpsi_fcgl):

• Current CGL formula (after PENTRC - BUG FIX - Recover full CGL curvature from perpendicular field strength #342): develop 2.5566e−4 → this PR 2.5497e−4.
• Pre-PENTRC - BUG FIX - Recover full CGL curvature from perpendicular field strength #342 formula quoted earlier: develop 2.6906e−4, old reader-only PR 2.6130e−4, this PR 2.6832e−4.

The reg_flag=f decks touched by #340 (input/gpec.in, a10_ideal, a10_kinetic, a5_tearing, solovev_kinetic) also get the new pmodb columns. Only the two Solov'ev examples were run.

Verification. Rebased on develop e0f0eb6. Native GNU Make dcon gpec pentrc builds pass with gfortran 16.2.1 and no new warnings. All eight DCON/GPEC/PENTRC example runs (two examples × two reg_flag settings × develop/PR) terminate normally. The eulb, lagb and xclebsch outputs are unchanged. No cluster or gfortran 12.2 build was run. No tests, CI targets or docs are added in this PR.

Evidence branch has the derivation, the residual budget, full per-m tables and the analysis script. The earlier reader evidence is in gpec16_validation.md.

Chris&AI

[logan-nc]

@krystophny, thanks for this. The physics premise is right, but merging this as-is would be a regression, because of how GPEC builds the column it now trusts.

What's right. δB_L = Q∥ + ξ·∇B = −B(∇·ξ⊥ + κ·ξ⊥) + ξ∥ b·∇B, and the ξ∥-free piece is the NTV drive. The default xclebsch path already uses it (inputs.f90:898). GPEC's lagb carries a small residual ξ∥: it is in the ξ^ζ=0 gauge, which is an O(ε²/q²) error, measured at 1.4–2.3% of max|lagb|.

The problem. gpout_pmodb builds Bdivxprp by subtracting ξ∥ from the full ξ (gpout.f:3030-3035, 3058-3060). That ξ∥ is the 1/(m−nq), Fourier-truncated closure solution, so the cancellation leaves a large residue in the resonant and window-edge harmonics. lagb never forms ξ∥. In today's develop output (regression_solovev_ideal, reg_flag=t) the three pmodb columns should satisfy the identity at every m:

m	Σ|Bdivxprp|²	Σ|lagb+Bkxprp|²	identity residual / Σ|lagb|²
0	3.80e-3	3.84e-3	0.000
−2	4.66e-5	5.32e-5	0.015
2 (q=2)	3.57e-2	8.03e-4	175
3 (q=3)	4.76e-3	3.87e-4	49
−4 (window edge)	7.27e-4	3.45e-5	88

To find which side is wrong, I rebuilt B∇·ξ⊥ independently from (ξ^ψ, ξ^s). That form has no ξ∥ and no 1/(m−nq), and is algebraically the same as PENTRC's smat…zmat matrices. Against it:

• GPEC's Bdivxprp is 10–16× too large at m=2,3, and stays flat through q=2 and q=3, where every other quantity dips.
• develop's −(lagb+Bkxprp) matches it to 0.5% in divx and 1–2% in dB/B. The 1–2% is just the gauge term.

Effect of this PR on what PENTRC sees (ψ-integrated Σ_m|X_m|², relative to develop): divx ×3.1 and dB/B ×7.2 on solovev_ideal (reg_flag=t), and up to ×69 and ×228 with reg_flag=f (solovev_kinetic). dbob_m drives every method, not only fcgl. The ψ=0.4 check sits between the q=2 and q=3 surfaces, where only non-resonant m contribute; those agree to 2–4%, which is why it didn't show up there. The 2016 "smaller values near rationals" comment was a workaround for exactly this.

A path forward. The fix behind this PR's intent is on the GPEC side. Build ξ⊥ in gpout_pmodb directly from (ξ^ψ, ξ^s):

ξ⊥^θ = (χ′²/(J²B²))[−(g12+q g13)ξ^ψ + (g23+q g33)ξ^s/χ′], ξ⊥^ζ = qξ⊥^θ − ξ^s/χ′

eqfun(2:3) already holds these metric combinations. Then the supplied divergence and −(lagb+Bkxprp) should agree to within the gauge term, and this reader change becomes a legitimate 1–2% correction.

Do you have the interest and capacity to try the GPEC-side fix in time for 1.6.0? Or are you OK leaving read_pmodb as it is on develop for now?

Related, on #340 (merged, and correct): "default unchanged" holds for the code default, but input/gpec.in and the a10_ideal, a10_kinetic, a5_tearing and solovev_kinetic examples set reg_flag=f, so their pmodb divergence output changed.

🤖 Generated with Claude Code

[logan-nc]

See claude details. Source annotations here already specified this was done on purpose and for backwards consistency. If you want to try the suggested larger lift, then go for it and I'll review. Otherwise, I'm fine with how it is and we can just close this one

[krystophny]

Thanks, agreed: trusting Bdivxprp as built today is a regression at the resonant and window-edge harmonics. I'll do the GPEC-side fix for 1.6.0: build ξ⊥ in gpout_pmodb directly from (ξ^ψ, ξ^s) with the eqfun(2:3) metric combinations, then check the per-m identity against −(lagb+Bkxprp) on solovev_ideal and solovev_kinetic for both reg_flag settings. This PR stays open until that lands, and then the reader change should come out as the 1–2% gauge correction. I'll also look at the reg_flag=f example outputs affected by #340.

Chris&AI

[krystophny]

The GPEC-side fix is in (cbf20c0): gpout_pmodb now forms ξ⊥ from (ξ^ψ, ξ^s) as you suggested, for both the divergence and the curvature column. The output format is unchanged.

Identity residual / Σ|lagb|² on solovev_ideal (reg_flag=t), develop → this PR:

m	develop	PR
2 (q=2)	175	3.9e-4
3 (q=3)	49	2.8e-4
−4 (window edge)	88	1.5e-3

solovev_kinetic and both reg_flag settings behave the same way. What remains is the expected gauge term: 0.04–0.12 at m=−4…−2 with reg_flag=f, and ≤1.5e-6 once that term is removed. PENTRC's integrated divx and dB/B now move by 1–4% relative to develop; the ideal reg_flag=f case is somewhat above the 1–2% you expected, and the PR body says so. Tables and the derivation are in the PR body and on the linked evidence branch.

Chris&AI

[krystophny]

The requested GPEC-side fix is in (cbf20c0) and CI is green, so this is ready for re-review for 1.6.0. If there isn't time for a careful look before Friday's cutoff, I'm happy for it to move to 1.6.1 instead, since it is a bug fix and fits the patch line, and I'd rather not rush a results-changing change in. #339 and #344 stay out of 1.6.0 as planned.

Chris&AI

[logan-nc]

Approving. Thanks @krystophny, this is the right fix.

Checked independently (hand derivation, a separate sympy red-team, and our own feynman runs, Solov'ev ideal/kinetic × reg_flag t/f plus a10_kinetic):

• Math. The new Jξ⊥ equals the old ξ − ξ∥b exactly (ξ^θ cancels identically), and the curvature rewrite is algebraically identical. Every difference from develop is the removed ξ∥ leakage.
• Accuracy. Bdivxprp matches an independent reconstruction from (ξ^ψ, ξ^s) to 0.06–0.5% per m. Develop was off 3.6–6.7× at m = 2, 3 and −4.
• PENTRC. The pmodb path now agrees with the default xclebsch path for every method in every case. For example, fgar went from 0.86% → 0.18% apart (ideal, reg=t) and from 1.17% → 0.26% (kinetic, reg=t); a10 rlar/clar/fcgl from 0.4–2.7% → ≤0.1%. Torque changes against develop match your −1…−5%.
• Scope. xclebsch, eulb, lagb and all non-pmodb outputs are bit-identical.

Pushed 4d7e278 on top (no change to any main output; gpec_pmodb, pmodbmn and the PENTRC results are bit-identical to cbf20c0):

• gpdiag_pmodb now uses the same direct ξ⊥ form. Its llagb vs −(Bdivx+Bkx) residual at resonant and edge m drops from O(1–100) to <1e-3.
• Dropped the now-unused xmt/xmz/xvt/xvz inverse transforms.
• Relabeled the PENTRC pmodb debug columns.

For the release notes (I'll add these):

• read_pmodb now trusts the supplied divergence, so a gpec_pmodb file written by an older GPEC still carries the leaky column. Regenerate pmodb files with 1.6 before using them with PENTRC 1.6.
• Pre-existing, not from this PR: with reg_flag=f the drive is unregularized and singular at rational surfaces. In the ideal example it is about 10⁶× the reg_flag=t drive at q=3, and results there depend on how the ψ grid samples the spike.

🤖 Generated with Claude Code