KineticForces - BUGFIX! - Give the energy integration its own tolerances (#392)

Author: logan-nc

src/KineticForces/CalculatedKineticMatrices.jl

@@ -136,7 +136,9 @@ function compute_calculated_kinetic_matrices(
                         z_s, m_s, kf_ctrl.wdfac, kf_ctrl.divxfac,
                         el_s, equil, intr_t, prof_s;
                         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

src/KineticForces/Compute.jl

@@ -137,7 +137,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

src/KineticForces/KineticForcesStructs.jl

@@ -121,11 +121,20 @@ builds a second control for its differing tolerance).
     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).

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
@@ -233,7 +240,7 @@ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int,
                                    B_vpar=B_vpar,
                                    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
@@ -897,7 +904,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
@@ -906,13 +915,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