From cb5e7c66533598bb375c8b6647283ccd3557e2dc Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sat, 25 Jul 2026 14:37:48 +0200 Subject: [PATCH 1/2] 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/2] 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