PENTRC: preserve periodicity in parallel-velocity spline

[krystophny]

What changed

Fit both PENTRC parallel-velocity root splines with periodic rather than
extrapolating endpoint conditions.

Root cause

The sampled quantity,

1.0 - (lambda / B0) * B(theta)

is periodic on the closed normalized poloidal-angle interval. The equilibrium
B(theta) spline immediately upstream is periodic, but the derived root spline
was fitted with nonperiodic extrapolation conditions. In the TC24 phiI010
L5 case the resulting seam corrupted bounce roots, produced non-finite
omega_b and omega_D, and reached LSODE as ISTATE=-5.

Validation

• Clean GNU Fortran build from develop:
  make -C install -j16 FC=gfortran CC=gcc OPENBLASHOME=/usr NETCDFHOME=/usr pentrc
• Exact previously failing TC24 phiI010 L5 input now terminates normally
  without an LSODE or non-finite-frequency failure.
• The repaired single-thread and 16-thread results are bit-for-bit identical
  for pentrc_tgar_ell_n3.out, pentrc_pgar_ell_n3.out,
  pentrc_tgar_n3.out, pentrc_pgar_n3.out, and pentrc_output_n3.nc.
• The same executable completed the phiI000 and phiI010 L0/L5 production
  matrix.

This changes only the spline boundary condition. It does not change a
frequency, torque, helicity, or profile sign.

[logan-nc]

@claude review this pr, paying special attention to if the spline endpoints are correctly inclusive or exclusive for what the periodic boundary condition expects

[claude]

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

---

Review: PENTRC periodic v‖ spline

• Gather context (diff, surrounding code)
• Understand spline endpoint semantics (inclusive vs exclusive) for periodic BC
• Check the two changed lines against how the root spline is sampled/evaluated
• Post review feedback

Verdict: the fix is correct, and the endpoints are handled correctly.

Why the endpoints line up

vspl is allocated with vspl%xs(:) = tspl%xs(:) (pentrc/torque.F90:427 and :537), and tspl%xs(0:) = i/mthsurf for i=0,...,mthsurf (pentrc/torque.F90:239). That's mthsurf+1 points spanning a closed, inclusive interval [0,1], where index 0 (θ=0) and index mthsurf (θ=1) are two different array slots holding the same physical point (since normalized θ is 2π-periodic). This is exactly the duplicate-endpoint convention spline_fit_ahg's "periodic" case expects — see equil/spline.f:245-250: it computes fs1(0,:) from fs(1,:), fs(0,:), and fs(spl%mx-1,:), solves the Sherman–Morrison periodic system on 0:mx-1, then forces fs1(mx,:) = fs1(0,:). That only produces a consistent (non-kinked) fit if fs(0,:) == fs(mx,:) going in, which is the same requirement the pre-existing tspl (torque.F90:266), dbdtspl (:291), and cglspl (:385) periodic fits already rely on.

vspl%fs(:,1) = 1.0-(lmda/bo)*tspl%fs(:,1) (torque.F90:462, :599) is a pointwise scalar transform of tspl%fs(:,1) — lmda and bo are both scalars in this context, not θ-dependent — so vspl%fs(0,1) == vspl%fs(mx,1) holds exactly whenever tspl%fs(0,1) == tspl%fs(mx,1) does, which is guaranteed because tspl is itself already fit with "periodic" just above. So the periodicity assumption isn't approximate — it's algebraically inherited from the already-periodic upstream B(theta) spline, matching the PR description.

Previously, "extrap" on this same duplicate-endpoint grid was the actual bug: it fits polynomial extrapolation boundary conditions across a grid where the first and last points are supposed to coincide, which is exactly the kind of endpoint mismatch that produces a spurious seam/kink at θ=0↔1 in the root-finding spline — consistent with the corrupted bounce roots and non-finite omega_b/omega_D described in the PR.

Other observations (non-blocking)

• spline_roots (equil/spline.f:1634, DO ix=0, lx with lx=spl%mx-1) walks every knot interval from xs(0) to xs(mx) once — it doesn't need special-casing for the wraparound because the duplicate endpoint is already baked into the grid. If a bounce point lands exactly on θ=0/θ=1 it could in principle be found once in the [xs(0),xs(1)) interval and once in [xs(mx-1),xs(mx)], but that's a pre-existing property of this duplicate-endpoint grid shared with tspl/dbdtspl/cglspl, not something this PR introduces or worsens.
• The change is minimal and scoped exactly as described (2 lines, boundary condition only); it doesn't touch bspl/bjspl/bwspl/tphi_spl/trans/kelmm bounce-time splines, which correctly remain "extrap" since those aren't periodic quantities in θ.
• .F90 free-form file, so the repo's 72-column fixed-format rule for .f files doesn't apply here.

No changes requested.

[logan-nc]

@claude fix the bug that xs(0) and xs(1) can't both be acceptable root solutions for periodic splines

[claude]

Claude finished @logan-nc's task in 4m 29s —— View job

---

Fixed: duplicate root at periodic spline wrap-around point

• Gather context (diff, spline_roots, spline_fit_ahg periodic case)
• Understand the bug: xs(0) and xs(mx) (the periodic wrap-around point, θ=0 and θ=1) could both be reported as separate roots
• Implement fix in equil/spline.f
• Manually verify syntax/logic (build tools unavailable in this sandbox — see note below)
• Push commit

Root cause

spline_roots (equil/spline.f:1592) walks each knot interval [xs(ix), xs(ix+1)] independently, ix=0,...,mx-1, and accepts a root whenever the normalized position z within that interval falls in [-eps, 1+eps]. For a periodic spline, xs(0) and xs(mx) are two different array slots holding the same physical point (e.g. normalized poloidal angle 0 and 1). If the sampled function has a root exactly at that shared point, the search finds it twice — once as z≈0 in the first interval (x=xs(0)), and again as z≈1 in the last interval (x=xs(mx)) — and nothing in the existing "reject repeated roots" logic (last1/last2/last3, equil/spline.f:1732-1774) catches this, because that dedup only compares against the immediately preceding accepted root, not across the whole domain.

This directly matches the case this PR is fixing: vspl (and also dbdtspl, pentrc/torque.F90:291, which finds extrema of dB/dθ) are now fit "periodic", so a bounce point or extremum sitting on the θ=0/1 seam is legitimately reachable and can be double-counted, corrupting the bounce-point count/pairing downstream.

Fix

In equil/spline.f:1776-1802, after the root-search loop completes, for periodic splines only: scan the collected roots for one within tolerance of xs(0) and one within tolerance of xs(mx); if both are present, drop the xs(mx) entry and decrement nroots. This is a minimal, generic fix in the shared spline_roots routine, so it also correctly benefits dbdtspl's periodic extrema search, not just the parallel-velocity spline in this PR. It has no effect on non-periodic ("extrap") splines, including the gpec/gpdiag.f self-test splines, which never set spl%periodic.

       INTEGER:: ix, jroot, lx, nzvalid, ir, ir0, irmx
       ...
       nroots = jroot
+      ! periodic splines duplicate the wrap-around point: xs(0) and
+      ! xs(mx) are the same physical location, so at most one may
+      ! survive as a root
+      IF(spl%periodic .AND. nroots>1)THEN
+         tol = eps*(spl%xs(1)-spl%xs(0))
+         ir0 = 0
+         DO ir=1,nroots
+            IF(ABS(roots(ir)-spl%xs(0))<tol)THEN
+               ir0 = ir
+               EXIT
+            ENDIF
+         ENDDO
+         tol = eps*(spl%xs(spl%mx)-spl%xs(spl%mx-1))
+         irmx = 0
+         DO ir=1,nroots
+            IF(ir/=ir0 .AND.
+     $         ABS(roots(ir)-spl%xs(spl%mx))<tol)THEN
+               irmx = ir
+               EXIT
+            ENDIF
+         ENDDO
+         IF(ir0>0 .AND. irmx>0)THEN
+            roots(irmx:nroots-1) = roots(irmx+1:nroots)
+            nroots = nroots - 1
+         ENDIF
+      ENDIF
       ! sort the roots lowest to highest

All added/modified lines stay within the 72-column fixed-format limit (verified with a regex check for lines ≥73 characters — only pre-existing long comments remain).

Note on verification

This sandbox's --allowedTools does not permit running gfortran/make (both required approval that wasn't grantable here), so I was unable to compile/run the TC24 phiI010 L5 case as done for the original PR. I manually re-verified the Fortran syntax (matching IF/ENDIF/ENDDO, correct 1-based array indexing consistent with the existing jroot/roots bookkeeping, and that array-section assignment on the assumed-size roots(*) dummy is legal with explicit bounds). If you'd like, please run the existing build/test workflow on this commit, or update my allowed tools to include gfortran/make so I can self-verify in a follow-up.

[logan-nc]

Verification of the periodic root dedup (e70a156f) + one hardening commit

The dedup commit was pushed without compile/run verification (sandbox limitation noted in the job log), so I built and tested it locally, and added one small hardening fix.

New commit: ccc31334

spline_fit set spl%periodic=.TRUE. for "periodic" but never cleared it on refit with another endmode. Previously this stale flag only affected interval wrapping in spline_eval, but the new dedup in spline_roots keys off it: refitting a previously-periodic spline object with "extrap" would silently drop a legitimate endpoint root. I reproduced this in a unit test before the fix. The commit resets the flag at the top of the spline_fit dispatcher (the only entry point — spline_fit_ha has no callers).

Verification done

• Unit test (sin(2πθ) on the closed [0,1] grid, mx=64): periodic fit returns 2 roots (0, 0.5) instead of the seam-duplicated 3; a fresh "extrap" fit still returns 3 roots including θ=1; an "extrap" refit of a previously periodic spline now also correctly returns 3.
• Clean gfortran build of dcon, gpec, pentrc at ccc31334.
• Full solovev_kinetic_example chain (dcon → gpec → pentrc): all exit 0, no warnings, all output variables finite.
• Baseline comparison: built pentrc at cb062a81 (periodic fit, no dedup) and ran it on identical inputs — pentrc_output_n1.nc is bit-for-bit identical across all 44 variables. The dedup is a pure safety net on this case; it only engages when a bounce point or dB/dθ extremum lands exactly on the θ=0/1 seam.

Remaining caveat

The TC24 phiI010 L5 validation in the PR description was run at cb062a81, i.e. before the dedup. Since dbdtspl (already periodic before this PR) can have its θ=0 extremum double-counted, the dedup can legitimately change results on cases where the seam is hit — the TC24 matrix and the 1-vs-16-thread bit-for-bit check should be rerun at ccc31334 before merge.

🤖 Generated with Claude Code

https://claude.ai/code/session_01WChssZ9PaArcmv9h583W9x