diff --git a/src/KineticForces/CalculatedKineticMatrices.jl b/src/KineticForces/CalculatedKineticMatrices.jl index ec0c0997a..c3fd3faaa 100644 --- a/src/KineticForces/CalculatedKineticMatrices.jl +++ b/src/KineticForces/CalculatedKineticMatrices.jl @@ -138,7 +138,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 diff --git a/src/KineticForces/Compute.jl b/src/KineticForces/Compute.jl index d0fdd591b..62c55d7f1 100644 --- a/src/KineticForces/Compute.jl +++ b/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 diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl index e4eca2581..0b58da3ec 100644 --- a/src/KineticForces/KineticForcesStructs.jl +++ b/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). diff --git a/src/KineticForces/Torque.jl b/src/KineticForces/Torque.jl index 1cd40285e..772310abf 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 @@ -231,7 +238,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 @@ -894,7 +901,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 @@ -903,13 +912,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