PENTRC: fix OpenMP data race on the energy contour flag (and a nested-quadrature tolerance defect)

[krystophny]

PENTRC aborts with lambdaintgrl_lsode - too many steps in lambda required on an ITER equilibrium, independently of compiler, tolerances, xmax, ximag, nl and kinetic profiles. The root cause is a data race.

The race

energy_imaxis selects the integration contour inside xintgrnd:

if(energy_imaxis)then
   cx = xj*x            ! imaginary axis
else
   cx = x + xj*ximag    ! real axis, Landau offset
endif

xintgrl_lsode flips it between its imaginary-axis and real-axis legs (energy.f90:190, :223). It is per-integration state exactly like energy_wn, energy_wd, energy_leff and the rest — and it is the only one of them missing from the threadprivate list.

Shared, one thread's imaginary leg silently switches another thread's integrand from x + i*ximag to i*x. exp(-cx) stops decaying and starts oscillating, so values from deep in the Maxwellian tail return order 1e4 instead of 1e-2. The pitch integrand is then not a function of lambda at all: in pentrc_lambdaintrgl_lsode.err, neighbouring lambda separated by 1e-8 differ by 130% on average, spanning -9.9e3 to +1.4e4 around a median of 9e-3. LSODE shrinks its step without bound and gives up.

This explains why the failure looked impossible to pin down: it does not depend on the numerics, it depends on the thread count. Decks that set the deprecated openmp_threads silently run on every available core.

Scope of the threadprivate claim

For the ordinary parallel path the fix is complete: energy_wn, energy_wt, energy_we, energy_wd, energy_wb, energy_nuk, energy_leff, energy_n and now energy_imaxis are all thread-private (energy.f90:62-81), and LSODE's /DLS002/ common block is thread-private too (lsode2.f:29-30).

It is not a claim that all of energy.f90 is thread-safe. xintgrl_lsode(..., op_record=.true.) calls record_method (energy.f90:293), which mutates the shared energy_record object (energy.f90:478-513) without synchronisation. Production recording is serial today, so this does not affect the demonstrated fix, but the recorded path should not be parallelised as it stands.

Evidence

Same deck, same binary, unmodified ITER case:

	result
before	1006 energy-integral failures, pitch integral aborts, no output
after, 1 thread	clean, writes pentrc_tgar_n3.out and pentrc_pgar_n3.out
after, 16 threads	clean, and bit-identical to the 1-thread result

Bit-identical output across thread counts is the correctness gate for a race fix, and it holds on all 29 columns of both output files.

Second commit: the nested quadrature shared one tolerance

pentrc_interface.f90 aliased xatol/xrtol to the same atol_xlmda/rtol_xlmda pair used for the pitch integral. The pitch integrand is the energy integral, so the outer integrator was asked to resolve its integrand to rtol_xlmda while that integrand was itself computed only to rtol_xlmda, and it then chases the inner integrator's own quadrature error. The energy integration now takes its own atol_x/rtol_x.

This changes default numerical behaviour, deliberately. With atol_x/rtol_x left at their derive-me sentinel, the energy tolerances become nested_tolerance_margin = 1e-2 times the pitch tolerances, so the shipped defaults move

	energy atol	energy rtol
before (aliased to pitch)	1e-6	1e-3
after (derived, 1e-2 margin)	1e-8	1e-5
energy.f90 module defaults, for reference	1e-12	1e-9

Every deck that does not set the new controls therefore integrates the energy variable more tightly than before, and pays for it in runtime. That is the intent — the previous nesting was inconsistent — but it is a behaviour change, not a no-op, and reviewers should decide whether 1e-2 is the margin they want. Note the derived values do not restore energy.f90's own 1e-12/1e-9; a deck can ask for those explicitly with atol_x/rtol_x.

Withdrawn from this PR

An earlier revision also removed the degenerate endpoints from the trapped-only pitch interval. That hunk is gone. The reasoning was wrong in general: lmdamin = max(1/(1+epsr), bo/bmax) and lmdamax = min(1/(1-epsr), bo/bmin) are degenerate only when the bo/bmax or bo/bmin term wins the extremum; when the cylindrical epsr bound wins, that endpoint can have two regular bounce points and deleting it discards nonsingular contributions. A correct fix has to test which term won. Filed separately rather than bundled here.

Built and tested with gfortran 15.

[krystophny]

Updated after an independent adversarial review and a rebase onto current develop (this branch was 12 commits behind). Three changes, all narrowing the PR:

1. Dropped the degenerate-endpoint change. The review was right and my justification was wrong. lmdamin = max(1.0/(1+epsr), bo/bmax) and lmdamax = min(1.0/(1-epsr), bo/bmin) are degenerate only when the bo/bmax or bo/bmin term wins. When the cylindrical epsr bound wins, that endpoint has two regular bounce points, and excluding it unconditionally would shrink the physical integration interval and discard nonsingular contributions. My comment claimed "both ends of the trapped pitch space are degenerate", which is not true in general. pentrc/torque.F90 is no longer touched by this PR.

2. The tolerance change is no longer described as behaviour-preserving. It is not. Every deck that does not set the new atol_x/rtol_x gets effective energy tolerances of 1e-8/1e-5 instead of the formerly aliased 1e-6/1e-3. A deck can restore the old values explicitly. Note also that energy.f90's own module defaults are tighter still (1e-12/1e-9); they 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 and is not settled here.

3. Split out the documentation changes. docs/sign_conventions.rst did not belong in a numerical fix, and bundling it was actively harmful: that text predated #268, so merging it would have reverted #268's correction that inverse_run recomputes q by flux-surface integration. The sign-convention corrections are now in a separate PR against current develop, carrying only the additive parts.

The race fix itself is unchanged and re-verified on the new base. Output is byte-identical to the pre-rebase build, and bit-identical at 1 and 16 threads.

One caveat added to the commit message, also from the review: recorded integrations (op_record=.true.) additionally mutate the shared energy_record object without synchronization, so xintgrl_lsode with recording is not generally thread-safe. Production recording is serial, so this does not affect the fix. I checked that the pitch-resolved netCDF, including omega_D_tgar, is bit-identical across thread counts, so no diagnostic output we rely on is affected either.

[logan-nc]

Approve. Reviewed the diff line-by-line and independently reproduced both the bug and the fix's correctness gates on a different equilibrium (the in-repo a10 case) than the ITER deck in the PR.

Race fix (commit 1) — verified correct and complete within its stated scope.

• energy_imaxis is assigned before use on every call (= .true. at energy.f90:198 only when ximag /= 0, unconditionally = .false. at :231 before the real leg), with first read only inside the subsequent LSODE legs — so thread-local initialization is guaranteed per call and no copyin is needed.
• Completeness sweep: every sibling per-integration module variable is already threadprivate (pitch_ell, pitch_n, all pitch_* in pitch.f90:70-73; all energy_*), both LSODE instances' common blocks are threadprivate (/DLS002/ lsode2.f:30, /DLS011/ lsode1.f:31), and the remaining shared energy-module state (xatol, xrtol, ximag, xmax, xnufac, ...) is set serially and read-only in the parallel region. energy_imaxis was genuinely the last mutable hole in this class. The caveat about the serial-only energy_record path is accurate.
• The parallel region (over bounce harmonics, torque.F90:1413) is reached by both the dynamic and equil/input grid paths (tintgrl_grid also calls tintgrnd), so the fix covers all of them.
• Independent reproduction: a10 kinetic deck, tgar, ximag = 0.05, nl = 2, psilims = 0.3 0.6. Develop at 8 threads dies with the exact reported signature — cascading LSODE- Energy- ... error test failed repeatedly, pentrc_lambdaintrgl_lsode.err written, STOP ERROR: lambdaintgrl_lsode - too many steps in lambda required, no output — while 1 thread runs clean. This branch at 1 and 8 threads: clean and bit-identical outputs. The race requires ximag /= 0 (line 198 is the only .true. write) plus nl >= 1, which explains both why it masqueraded as a numerics problem and why decks with ximag = 0 never see it.

Tolerance split (commit 2) — verified a faithful re-plumbing, and the behavior change is empirically justified.

• Both namelist-reader sites are converted; xatol/xrtol have no other writers; the -1 sentinel semantics are safe (explicit 0 and the old aliased values remain expressible).
• Exact no-op when pinned: this branch with atol_x/rtol_x set to the old aliased values is bit-identical to a build carrying only the race-fix commit. The sentinel logic reproduces the old behavior exactly on request.
• The default change is doing real work: on the a10 case, tightening the energy tolerances by the 1e-2 margin moves the integrated torque by 0.64% — ~60x the requested rtol_psi = 1e-4 — i.e. the old aliased nesting was producing errors well beyond the accuracy the outer integrals claimed. The new default is not just tidier, it is necessary for the outer tolerances to mean what they say. The 1e-2 margin looks right.
• One measurement note for other reviewers: cross-build outputs differ at the ~1e-4 level even at 1 thread with identical tolerances (recompiling the modified translation unit under -O3 changes instruction scheduling/FMA fusion). That is codegen noise at the integration-tolerance level, not a semantic change — the meaningful gate is same-binary thread invariance, which holds.

Minor asks (non-blocking):

1. nested_tolerance_margin = 1e-2 — make it 1e-2_r8 per the project's precision convention.
2. The new atol_x/rtol_x controls are not documented in input/pentrc.in or any example deck, and every shipped pentrc.in still comments atol_xlmda as "tolerance in energy and pitch integration", which this PR makes untrue. Worth updating the template comment and adding the new controls in this PR or a quick follow-up.
3. The 0.64%-scale default shift will surface in any regression comparison against old PENTRC outputs — flagging for the eventual release notes.

🤖 Generated with Claude Code

https://claude.ai/code/session_01Rewp8poxoQyorK48SwxEUv