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.