From 3b8e873379cee10ef0cd77a92a27731246795447 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Fri, 6 Nov 2026 22:26:18 -0500 Subject: [PATCH] KineticForces - BUGFIX! - Start the torque quadrature where the perturbation data starts With qlow above q0, ForceFreeStates raises psilow and the plasma response starts off axis, but the post-PE torque quadrature still started at the equilibrium's psilow and integrated extrapolated perturbation data below the response. It now starts at the first psi of that data, as it already stopped at psilim. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01BgNxCFidkS5NEJM3EvSoBs --- src/KineticForces/Compute.jl | 6 +++--- src/KineticForces/KineticForcesStructs.jl | 3 +++ test/runtests_solve_api.jl | 7 +++++++ 3 files changed, 13 insertions(+), 3 deletions(-) diff --git a/src/KineticForces/Compute.jl b/src/KineticForces/Compute.jl index dd087dc94..506fac24c 100644 --- a/src/KineticForces/Compute.jl +++ b/src/KineticForces/Compute.jl @@ -84,9 +84,9 @@ function integrate_psi_quadgk( is_matrix_method = occursin("mm", method) mpert = intr.mpert - # Start at psilow (degenerate bounce/drift frequencies at the axis); stop at intr.psilim, where the - # perturbation data ends. - x0 = max(psi_min, equil.config.psilow) + # Start at psilow (degenerate bounce/drift frequencies at the axis) or where the perturbation data starts, if later; + # stop at intr.psilim, where it ends. + x0 = max(psi_min, equil.config.psilow, intr.psilow) xout = min(psi_max, intr.psilim, 1.0 - 1e-6) if x0 >= xout diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl index 33cadb66b..dc8d05494 100644 --- a/src/KineticForces/KineticForcesStructs.jl +++ b/src/KineticForces/KineticForcesStructs.jl @@ -199,6 +199,8 @@ perturbation and geometric-matrix interpolants. Equilibrium and kinetic profiles # Upper ψ bound (FFS psilim); the perturbation interpolants are invalid beyond it, # so the outer ψ quadrature clips here. psilim::Float64 = 1.0 + # Lower ψ bound: the first ψ of the perturbation data, above psilow when FFS starts at qlow. + psilow::Float64 = 0.0 # Rational-surface ψ locations for ψ-quadrature paneling (see docstring). sing_psis::Vector{Float64} = Float64[] @@ -281,6 +283,7 @@ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state::Pertur xi_modes = pe_state.xi_modes psi_grid = pe_state.psi_grid npsi = length(psi_grid) + kf_intr.psilow = psi_grid[1] mpert = ffs.mpert # Build xs_m: 3 CubicSeriesInterpolants from Clebsch displacement matrices diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl index 0b03b4c30..f4d574709 100644 --- a/test/runtests_solve_api.jl +++ b/test/runtests_solve_api.jl @@ -627,6 +627,13 @@ using HDF5 @test read(api_h5["KineticForces/Validity/psi_c"]) == read(deck_h5["KineticForces/Validity/psi_c"]) end end + + # With qlow above q0 the response starts off axis, and the torque quadrature starts with it. + ntv_deck["ForceFreeStates"]["qlow"] = 1.9 + open(io -> TOML.print(io, ntv_deck), joinpath(dir, "gpec.toml"), "w") + raised = GPEC.main([dir]) + @test raised.pe.psi_grid[1] > ntv_deck["Equilibrium"]["psilow"] + @test minimum(raised.kinetic.total.method_results["fgar"].psi_grid) >= raised.pe.psi_grid[1] end end end