From b69e40ce5e17a023285372fb19572e3e3c907c52 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 15 Aug 2026 14:07:57 -0400 Subject: [PATCH] KineticForces - IMPROVEMENT - give the energy integration its own tolerances The energy (x) integration was handed the same atol_xlmda/rtol_xlmda pair as the pitch (lambda) integration that encloses it. Since the pitch integrand IS the energy integral, the outer integrator was asked to resolve its integrand to the same tolerance to which that integrand was itself computed, so it chases the inner integrator's quadrature noise instead of converging. Adds atol_x/rtol_x, defaulting (NaN sentinel) to nested_tolerance_margin = 1e-2 times the pitch tolerances, which extends the nesting rule the struct already documents one level up for rtol_psi vs rtol_xlmda. Wired through both the psi-quadrature path (tpsi!) and the kinetic-matrix path. This deliberately changes default numerical behaviour: shipped decks now integrate the energy variable to atol 1e-10 / rtol 1e-7 rather than 1e-8 / 1e-5, and pay for it in runtime. A deck can set atol_x/rtol_x explicitly to override the derived values, or widen nested_tolerance_margin to recover the old cost. Julia analog of the second commit of Fortran GPEC PR #280. Co-Authored-By: Claude Fable 5 --- .../CalculatedKineticMatrices.jl | 4 +++- src/KineticForces/Compute.jl | 4 +++- src/KineticForces/KineticForcesStructs.jl | 19 ++++++++++++----- src/KineticForces/Torque.jl | 21 +++++++++++++++---- 4 files changed, 37 insertions(+), 11 deletions(-) diff --git a/src/KineticForces/CalculatedKineticMatrices.jl b/src/KineticForces/CalculatedKineticMatrices.jl index d36fead99..66e969287 100644 --- a/src/KineticForces/CalculatedKineticMatrices.jl +++ b/src/KineticForces/CalculatedKineticMatrices.jl @@ -121,7 +121,9 @@ function compute_calculated_kinetic_matrices( kf_ctrl.zi, kf_ctrl.mi, kf_ctrl.wdfac, kf_ctrl.divxfac, kf_ctrl.electron, equil, intr_t, kinetic_profiles; nutype=kf_ctrl.nutype, f0type=kf_ctrl.f0type, nufac=kf_ctrl.nufac, - atol_xlmda=kf_ctrl.atol_xlmda, rtol_xlmda=kf_ctrl.rtol_xlmda + atol_xlmda=kf_ctrl.atol_xlmda, rtol_xlmda=kf_ctrl.rtol_xlmda, + atol_x=kf_ctrl.atol_x, rtol_x=kf_ctrl.rtol_x, + nested_tolerance_margin=kf_ctrl.nested_tolerance_margin ) full_w .+= block_w full_t .+= block_t diff --git a/src/KineticForces/Compute.jl b/src/KineticForces/Compute.jl index 4175a62b4..7a316a5cd 100644 --- a/src/KineticForces/Compute.jl +++ b/src/KineticForces/Compute.jl @@ -135,7 +135,9 @@ function integrate_psi_quadgk( tpsi!(thread_tpsi[tid], psi, n, l, zi, mi, wdfac, divxfac, electron, method, equil, thread_intrs[tid], kinetic_profiles; op_wmats=w, - atol_xlmda=ctrl.atol_xlmda, rtol_xlmda=ctrl.rtol_xlmda) + atol_xlmda=ctrl.atol_xlmda, rtol_xlmda=ctrl.rtol_xlmda, + atol_x=ctrl.atol_x, rtol_x=ctrl.rtol_x, + nested_tolerance_margin=ctrl.nested_tolerance_margin) harm_vals[ell_idx] = thread_tpsi[tid][] is_matrix_method && (harm_elems[ell_idx] .= w) end diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl index 6df40e6a8..b153e39ae 100644 --- a/src/KineticForces/KineticForcesStructs.jl +++ b/src/KineticForces/KineticForcesStructs.jl @@ -95,11 +95,20 @@ ctrl = KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForce nn::Int = 1 # Toroidal mode number nl::Int = 1 # Bounce harmonic number - # Tolerances. - # *_xlmda: shared tolerances for inner λ (pitch) and x (energy) integrations - # *_psi: tolerances for outer ψ quadrature - atol_xlmda::Float64 = 1e-8 # Absolute tolerance for inner pitch + energy integrations - rtol_xlmda::Float64 = 1e-5 # Relative tolerance for inner pitch + energy integrations + # Tolerances, outermost to innermost: ψ quadrature ⊃ λ (pitch) ⊃ x (energy). + # Each level must be resolved more tightly than the one enclosing it, or the outer + # integrator chases its integrand's own quadrature noise instead of converging. + # *_xlmda: tolerances for the λ (pitch) integration + # *_x: tolerances for the x (energy) integration nested inside it; NaN ⇒ derive as + # nested_tolerance_margin × the pitch tolerances + # *_psi: tolerances for the outer ψ quadrature + atol_xlmda::Float64 = 1e-8 # Absolute tolerance for the inner pitch integration + rtol_xlmda::Float64 = 1e-5 # Relative tolerance for the inner pitch integration + atol_x::Float64 = NaN # Absolute tolerance for the energy integration (NaN ⇒ derived) + rtol_x::Float64 = NaN # Relative tolerance for the energy integration (NaN ⇒ derived) + # The pitch integrand IS the energy integral, so the energy level is resolved this much + # tighter than the pitch level by default. + nested_tolerance_margin::Float64 = 1e-2 # Factor relating derived energy tolerances to the pitch ones # rtol_psi is the primary convergence knob: ~2 significant figures matches the validity # of the NTV model approximations. Do not set it tighter than the noise floor of the # inner integrals (keep rtol_psi ≳ 10 × rtol_xlmda). diff --git a/src/KineticForces/Torque.jl b/src/KineticForces/Torque.jl index f9af86485..5fd0c9475 100644 --- a/src/KineticForces/Torque.jl +++ b/src/KineticForces/Torque.jl @@ -35,7 +35,14 @@ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int, op_wmats::Union{Nothing,Array{ComplexF64,3}}=nothing, rex_override::Union{Nothing,Float64}=nothing, imx_override::Union{Nothing,Float64}=nothing, - atol_xlmda::Float64=1e-9, rtol_xlmda::Float64=1e-6) + atol_xlmda::Float64=1e-9, rtol_xlmda::Float64=1e-6, + atol_x::Float64=NaN, rtol_x::Float64=NaN, + nested_tolerance_margin::Float64=1e-2) + + # The pitch integrand is itself the energy integral, so the energy level is resolved + # tighter than the pitch level (explicit atol_x/rtol_x override the derived values). + atol_energy = isnan(atol_x) ? nested_tolerance_margin * atol_xlmda : atol_x + rtol_energy = isnan(rtol_x) ? nested_tolerance_margin * rtol_xlmda : rtol_x # Enforce bounds if psi > 1 @@ -229,7 +236,7 @@ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int, B_extrap=B_extrap, smat=smat_f, tmat=tmat_f, xmat=xmat_f, ymat=ymat_f, zmat=zmat_f, - energy_atol=atol_xlmda, energy_rtol=rtol_xlmda, + energy_atol=atol_energy, energy_rtol=rtol_energy, pitch_atol=atol_xlmda, pitch_rtol=rtol_xlmda, rex_override=rex_override, imx_override=imx_override) end @@ -892,7 +899,9 @@ function compute_kinetic_matrices_at_psi!( electron::Bool, equil, intr::KineticForcesInternal, kinetic_profiles::Equilibrium.KineticProfileSplines; nutype::String="harmonic", f0type::String="maxwellian", nufac::Float64=1.0, - atol_xlmda::Float64=1e-9, rtol_xlmda::Float64=1e-6) + atol_xlmda::Float64=1e-9, rtol_xlmda::Float64=1e-6, + atol_x::Float64=NaN, rtol_x::Float64=NaN, + nested_tolerance_margin::Float64=1e-2) # Bypass ψ > 1 (no kinetic contribution outside plasma) if psi > 1 @@ -901,13 +910,17 @@ function compute_kinetic_matrices_at_psi!( return nothing end + # See tpsi!: the energy integral is the pitch integrand, so it is resolved tighter. + atol_energy = isnan(atol_x) ? nested_tolerance_margin * atol_xlmda : atol_x + rtol_energy = isnan(rtol_x) ? nested_tolerance_margin * rtol_xlmda : rtol_x + state = _setup_surface_state(psi, zi, mi, electron, equil, intr, kinetic_profiles) kinetic_energy_matrices_for_euler_lagrange!( kwmat, ktmat, state, psi, n, l, wdfac, intr; nutype, f0type, nufac, - energy_atol=atol_xlmda, energy_rtol=rtol_xlmda, + energy_atol=atol_energy, energy_rtol=rtol_energy, pitch_atol=atol_xlmda, pitch_rtol=rtol_xlmda) return nothing