From cb5e7c66533598bb375c8b6647283ccd3557e2dc Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sat, 25 Jul 2026 14:37:48 +0200 Subject: [PATCH 1/4] PENTRC: fix OpenMP data race on the energy contour flag energy_imaxis selects the integration contour inside xintgrnd -- cx = i*x on the imaginary-axis leg against cx = x + i*ximag on the real one -- and xintgrl_lsode flips it between its two legs (energy.f90:190, :223). It is per-integration state exactly like energy_wn, energy_wd, energy_leff and the rest, and it was the only one of them missing from the threadprivate list. Shared, one thread's imaginary leg silently switches another thread's integrand from a decaying exponential to an oscillation, so values from deep in the Maxwellian tail come back four orders of magnitude too large and the pitch integrand stops being a function of lambda: neighbouring lambda separated by 1e-8 differ by 130% on average. lambdaintgrl_lsode then shrinks its step without bound and aborts with 'too many steps in lambda required'. That is why the failure resisted every numerical explanation -- tolerances, xmax, ximag, nl, compiler, kinetic profiles: it does not depend on the numerics, it depends on the thread count. Decks that set the deprecated and ignored openmp_threads = 1 silently ran on every available core. Gate: the ITER case that previously aborted now runs unmodified and gives bit-identical pentrc_tgar_n3.out and pentrc_pgar_n3.out at 1 and 16 threads, all 29 columns, including the pitch-resolved netCDF. Note for reviewers: recorded integrations (op_record=.true.) additionally mutate the shared energy_record object without synchronization. Production recording is serial so this does not affect the fix above, and it is left alone here. --- pentrc/energy.f90 | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/pentrc/energy.f90 b/pentrc/energy.f90 index ba1f9f1a7..e0a70b14e 100644 --- a/pentrc/energy.f90 +++ b/pentrc/energy.f90 @@ -69,7 +69,15 @@ module energy_integration energy_wb,& energy_nuk,& energy_leff - !$omp threadprivate(energy_wn,energy_wt,energy_we,energy_wd,energy_wb,& + ! energy_imaxis selects the integration contour inside xintgrnd, and + ! xintgrl_lsode flips it between its imaginary-axis and real-axis legs. It + ! is per-integration state exactly like the frequencies below it, so it must + ! be threadprivate too: shared, one thread's imaginary leg silently switches + ! another thread's integrand from x + i*ximag to i*x, turning exp(-cx) from + ! a decaying exponential into an oscillation and injecting values orders of + ! magnitude too large into the pitch integrand. + !$omp threadprivate(energy_imaxis,& + !$omp& energy_wn,energy_wt,energy_we,energy_wd,energy_wb,& !$omp& energy_nuk,energy_leff,energy_n) type record From f2cd94071310b7602406dde5b37fc689b223c31d Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sat, 25 Jul 2026 14:38:13 +0200 Subject: [PATCH 2/4] PENTRC: give the energy integration its own tolerances pentrc_interface aliased both xatol/xrtol and lambdaatol/lambdartol to the single atol_xlmda/rtol_xlmda pair. The pitch integral's integrand IS the energy integral, so this is a nested adaptive quadrature, and with one tolerance for both the outer integrator chases the inner one's quadrature error -- which is not a smooth function of lambda. LSODE then shrinks its step without bound. New pent_control entries atol_x/rtol_x carry the energy tolerances. Their default of -1 means 'derive from the pitch tolerances', as nested_tolerance_margin = 1e-2 times them. This DOES change numerical behaviour for every deck that does not set the new controls: the effective energy tolerances go from the formerly aliased 1e-6/1e-3 to 1e-8/1e-5, i.e. tighter. Decks can restore the old values explicitly by setting atol_x = 1e-6, rtol_x = 1e-3. Note that energy.f90's own module defaults are tighter still, 1e-12/1e-9; those 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. This is NOT what unblocked the previously-aborting ITER case; the OpenMP race on energy_imaxis was. It is a real defect in its own right and is kept separate. Dropped from an earlier version of this change: an exclusion of both endpoints of the RLAR trapped pitch space. lmdamin = max(1/(1+epsr), bo/bmax) and lmdamax = min(1/(1-epsr), bo/bmin) are degenerate only when the bo/b term wins; when the cylindrical epsr bound wins the endpoint has two regular bounce points, and excluding it unconditionally would shrink the physical interval and discard nonsingular contributions. --- pentrc/pentrc_interface.f90 | 46 +++++++++++++++++++++++++++++++++---- 1 file changed, 41 insertions(+), 5 deletions(-) diff --git a/pentrc/pentrc_interface.f90 b/pentrc/pentrc_interface.f90 index 8b890b4d9..35571f139 100644 --- a/pentrc/pentrc_interface.f90 +++ b/pentrc/pentrc_interface.f90 @@ -111,9 +111,15 @@ module pentrc_interface pentrc_threads = 0,& openmp_threads = 0 + !> How much tighter the energy integration is than the pitch integration + !> when atol_x/rtol_x are left at their derive-me default. + real(r8), parameter :: nested_tolerance_margin = 1e-2 + real(r8) :: & atol_xlmda=1e-6, & rtol_xlmda=1e-3, & + atol_x=-1, & + rtol_x=-1, & nfac=1.0, & tfac=1.0, & wefac=1.0, & @@ -148,7 +154,8 @@ module pentrc_interface jac_in, jsurf_in, tmag_in, power_bin, power_bpin, power_rin, power_rcin namelist/pent_control/nfac, tfac, wefac, wdfac, wpfac, nufac, divxfac, & - atol_xlmda, rtol_xlmda, atol_psi, rtol_psi, nlmda, ntheta, ximag, xmax, psilims, & + atol_xlmda, rtol_xlmda, atol_x, rtol_x, atol_psi, rtol_psi, & + nlmda, ntheta, ximag, xmax, psilims, & use_classic_splines,pentrc_threads,openmp_threads,force_xialpha namelist/pent_output/moment, output_ascii, output_netcdf, & @@ -221,8 +228,7 @@ subroutine initialize_pentrc(op_kin,op_deq,op_peq) !if(any(psilim/=0)) print *, "!! WARNING: psilim has been deprecated. Use psilims." ! distribute some simplified inputs to module circles - xatol = atol_xlmda - xrtol = rtol_xlmda + call set_nested_tolerances xnufac= nufac xnutype= nutype xf0type= f0type @@ -279,8 +285,7 @@ subroutine get_pentrc(get_nl,get_zi,get_mi,get_wdfac,get_divxfac,& close(i) ! distribute inputs to PENTRC module circles - xatol = atol_xlmda - xrtol = rtol_xlmda + call set_nested_tolerances xnufac= nufac xnutype= nutype xf0type= f0type @@ -307,6 +312,37 @@ subroutine get_pentrc(get_nl,get_zi,get_mi,get_wdfac,get_divxfac,& end subroutine get_pentrc + + !======================================================================= + subroutine set_nested_tolerances + !----------------------------------------------------------------------- + !*DESCRIPTION: + ! Give the energy integration its own, tighter tolerances. + ! + ! The pitch integral's integrand IS the energy integral, so this is a + ! nested adaptive quadrature. With one tolerance for both, the outer + ! integrator chases the inner one's quadrature error, which is not a + ! smooth function of lambda; LSODE then shrinks its step without bound + ! and stops with "too many steps in lambda required". + ! + ! Negative atol_x/rtol_x, the default, means "derive from the pitch + ! tolerances". A deck may set them explicitly, including back to the + ! old aliased values. + !----------------------------------------------------------------------- + implicit none + if(atol_x < 0) then + xatol = atol_xlmda * nested_tolerance_margin + else + xatol = atol_x + endif + if(rtol_x < 0) then + xrtol = rtol_xlmda * nested_tolerance_margin + else + xrtol = rtol_x + endif + end subroutine set_nested_tolerances + + end module pentrc_interface From 1ad89ff386e2e3c9a7f9ea19a831254b28893e00 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sat, 25 Jul 2026 14:10:01 +0200 Subject: [PATCH 3/4] PENTRC: omega_D in the GAR branch carried a spurious major radius wdbar reused wbbar, which is omega_b/bhat and therefore carries one factor of ro that bhat takes back out again. dhat cancels only the ro**2 written explicitly, so the third ro survived and omega_D came out with units of m/s rather than 1/s. The operator itself is correct and is not touched: bspl%fsi(:,2) is exactly d/dPsi of the parallel action integrand at fixed lmda -- the 1 - 1.5*lmda*b/bo coefficient follows from sqrt(vpar) = vpar/sqrt(vpar) in the product rule -- so the bounce average is the canonical precession -(1/chrg)(dJ_par/dpsi_p)/tau_b = -4 pi (T/chrg) I2/I1, the 4 pi because chi1 = twopi*psio makes the tabulated derivatives d/dPsi_pol while the canonical momentum wants d/dpsi_p. Only the prefactor changes: ro, not ro*ro. Verified against an independently validated canonical operator on the same equilibrium, as the least-squares slope over all pitches on a surface (a per-point ratio is ill conditioned at the omega_D zero crossing): -6.494 before against the predicted -ro = -6.411, and -1.021 after against -1, over 7 surfaces. omega_b, built from the same bspl(:,1), agrees to better than 0.4% before and after and is untouched. The dimensional argument is what makes this a length rather than a missing 2 pi, which is only 2% away numerically: I2/I1 has units 1/Wb whatever normalisation b and J carry, since any constant rescaling of J*b cancels in the ratio, so (T/chrg) I2/I1 is already a frequency. Confirmed that no hidden normalization supplies the extra radius: rzphi stores physical R**2 and a physical-coordinate Jacobian, eqfun stores B in tesla, ro is in metres. A second witness is internal: omega_d_rlar in the same output file is wdhat = q T/(2 eps ro**2 chrg bo), dimensionally a frequency, reproducing to 8 digits, and it sits at 0.78x the canonical operator while the GAR bounce average sat at 8.26x it. SCOPE OF THE SIGN CLAIM. Only the magnitude is established here. The overall sign of the canonical form depends on whether the thin-orbit canonical momentum is +chrg*psi_p or -chrg*psi_p, and GPEC cannot settle that from its stored radial coordinate: EFIT input signs are normalized away in equil/read_eq.f and the native toroidal coordinate reaches the laboratory frame only through phi_CCW = -helicity*(2 pi zeta + dphi) in coil/field.F. This change removes one positive factor of ro and nothing else; no sign is altered, and the residual orientation question is independent of it. On the ITER case this moves the trapped superbanana resonance from x ~ 4.3 to x ~ 28, so TGAR torque falls from -11.8 to -0.5 N m and the balance becomes passing-dominated. Output stays bit-identical at 1 and 16 threads. --- pentrc/torque.F90 | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/pentrc/torque.F90 b/pentrc/torque.F90 index 9d7426673..0522d30f1 100644 --- a/pentrc/torque.F90 +++ b/pentrc/torque.F90 @@ -738,7 +738,14 @@ subroutine tpsi(tpsi_var,psi,n,l,zi,mi,wdfac,divxfac,electron, & ! Bounce averaged Lambda functions fbnce%xs(ilmda-1) = lmda wbbar = ro*twopi/((2-sigma)*bspl%fsi(bspl%mx,1)) - wdbar = ro*ro*bo*wdfac*wbbar*2*(2-sigma)*bspl%fsi(bspl%mx,2) + ! wbbar is omega_b/bhat, so it carries one factor of ro that + ! bhat takes back out. dhat removes only the explicit ro**2, + ! so a third ro from wbbar would survive and leave omega_D a + ! velocity rather than a frequency. Hence ro, not ro*ro: + ! wdbar*dhat = 4 pi (T/chrg) I2/I1 = -(1/chrg)(dJ/dpsi)/tau_b, + ! the canonical precession, with I1 in bspl%fsi(:,1) and I2 + ! in bspl%fsi(:,2) and d/dpsi taken against chi1. + wdbar = ro*bo*wdfac*wbbar*2*(2-sigma)*bspl%fsi(bspl%mx,2) bhat = sqrt(2*kin_f(s+2)/mass)/ro dhat = (kin_f(s+2)/chrg)/(bo*ro*ro) fbnce%fs(ilmda-1,1) = wbbar*bhat From 2e031ae675205ec9d69bb0298855155737f72579 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sat, 25 Jul 2026 16:14:54 +0200 Subject: [PATCH 4/4] PENTRC: state the precession identity as a magnitude in the comment too The PR body was corrected to claim only a magnitude, since the familiar minus sign needs P_phi = +chrg*psi_p and GPEC does not establish that convention. The source comment still asserted the signed identity, which is the exact claim the body disclaims. No executable change. --- pentrc/torque.F90 | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/pentrc/torque.F90 b/pentrc/torque.F90 index 0522d30f1..b3ead346d 100644 --- a/pentrc/torque.F90 +++ b/pentrc/torque.F90 @@ -742,9 +742,14 @@ subroutine tpsi(tpsi_var,psi,n,l,zi,mi,wdfac,divxfac,electron, & ! bhat takes back out. dhat removes only the explicit ro**2, ! so a third ro from wbbar would survive and leave omega_D a ! velocity rather than a frequency. Hence ro, not ro*ro: - ! wdbar*dhat = 4 pi (T/chrg) I2/I1 = -(1/chrg)(dJ/dpsi)/tau_b, + ! |wdbar*dhat| = 4 pi (T/chrg) |I2/I1| = |(1/chrg)(dJ/dpsi)/tau_b|, ! the canonical precession, with I1 in bspl%fsi(:,1) and I2 - ! in bspl%fsi(:,2) and d/dpsi taken against chi1. + ! in bspl%fsi(:,2) and d/dpsi taken against chi1. Magnitudes + ! only: the familiar minus sign needs P_phi = +chrg*psi_p, + ! and GPEC does not fix that convention here -- EFIT signs + ! are normalised away in equil/read_eq.f and the toroidal + ! angle carries a helicity declared in coil.in. This change + ! multiplies by a positive number and alters no sign. wdbar = ro*bo*wdfac*wbbar*2*(2-sigma)*bspl%fsi(bspl%mx,2) bhat = sqrt(2*kin_f(s+2)/mass)/ro dhat = (kin_f(s+2)/chrg)/(bo*ro*ro)