KineticForces - BUGFIX! - Start the torque quadrature where the perturbation data starts

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

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

test/runtests_solve_api.jl

@@ -636,6 +636,13 @@ include("display_check.jl")
                         @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