ForceFreeStates - BUGFIX! - Give every Euler-Lagrange solve an explicit absolute tolerance and default to Vern7

examples/DIIID-like_ideal_example/gpec.toml

@@ -45,6 +45,8 @@ mthvac = 512                   # Number of points used in splines over poloidal
 kinetic_source = "fixed"       # Kinetic matrix source: "fixed" test matrices, or "calculated" from the kinetic NTV model
 kinetic_factor = 0.0           # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode)
 eulerlagrange_tolerance = 1e-10 # Relative tolerance for ODE integration of Euler-Lagrange equations
+ode_abstol = 1e-8              # Absolute tolerance for the Euler-Lagrange sweeps and propagator chunks
+ode_solver = "Vern7"           # Explicit Runge-Kutta method for every Euler-Lagrange solve: Vern6, Vern7, Vern8, Vern9 or DP8
 save_interval = 3              # Save every Nth ODE step (1=all). Always saves near rational surfaces.
 singfac_min = 1e-4             # Fractional distance from rational q at which ideal jump enforced
 ucrit = 1e4                    # Column-norm threshold that triggers solution renormalization

src/ForceFreeStates/CoreTypes.jl

@@ -131,6 +131,8 @@ gpec.toml.
   - `nstep::Int` - Maximum number of integration steps (not yet implemented)
   - `ksing::Int` - Singular surface handling parameter
   - `eulerlagrange_tolerance::Float64` - Relative tolerance for ODE integration of Euler-Lagrange equations
+  - `ode_abstol::Float64` - Absolute tolerance for the forward sweep, the Riccati outer plasma and the propagator chunks. Default `1e-8`: the OrdinaryDiffEq default of `1e-6` lets small state entries escape the relative control, so the error stops responding to `eulerlagrange_tolerance` below about `1e-8`. The Δ′ shooting solves do not use it; they scale a per-column absolute tolerance from `eulerlagrange_tolerance`.
+  - `ode_solver::String` - Explicit Runge-Kutta method for every Euler-Lagrange solve, one of `"Vern6"`, `"Vern7"`, `"Vern8"`, `"Vern9"`, `"DP8"` (the benchmarked set, [`EL_ODE_SOLVERS`](@ref)). Default `"Vern7"`: at equal tolerances it reaches the same δW and a closer Δ′ than `"Vern9"` for fewer RHS evaluations, because the ninth-order method rejects most of its steps near the rational surfaces.
   - `ucrit::Float64` - Critical value of unorm ratio to trigger solution normalization. In the standard path it triggers Gaussian reduction; in the Riccati path it triggers `renormalize_riccati_inplace!`. Default `1e4` empirically keeps max(|U₁|, |U₂|) in O(1)–O(10⁴) over the integration domain on DIII-D / Solovev sweeps; lower triggers excess renorms without accuracy gain, higher risks overflow before the next renorm.
   - `numsteps_init::Int` - Initial array size for ODE data storage
   - `numunorms_init::Int` - Initial array size for solution normalization data
@@ -168,6 +170,8 @@ gpec.toml.
     nstep::Int = typemax(Int)
     ksing::Int = -1
     eulerlagrange_tolerance::Float64 = 1e-8
+    ode_abstol::Float64 = 1e-8
+    ode_solver::String = "Vern7"
     ucrit::Float64 = 1e4
     numsteps_init::Int = 4000
     numunorms_init::Int = 100

src/ForceFreeStates/EulerLagrange.jl

@@ -320,17 +320,27 @@ Only the Riccati branch populates `propagators` / `chunks` / `S_left`, which
 for all three.
 """
 function eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, intr::ForceFreeStatesInternal)
-
-    if ctrl.integrator == "riccati"
-        ctrl.kinetic_factor > 0 && error("kinetic runs require integrator=\"forward\"; the Riccati integrator has no kinetic crossing.")
-        return riccati_eulerlagrange_integration(ctrl, equil, mats, intr)
-    elseif ctrl.integrator == "forward"
-        return forward_eulerlagrange_integration(ctrl, equil, mats, intr)
-    elseif ctrl.integrator == "galerkin"
+    ctrl.integrator == "galerkin" &&
         error("integrator = \"galerkin\" solves the Euler-Lagrange system variationally, not by ODE integration; " *
               "it is dispatched to galerkin_solve.")
+    ctrl.integrator in ("riccati", "forward") ||
+        error("Unknown integrator: $(ctrl.integrator). Expected \"forward\", \"riccati\", or \"galerkin\".")
+    ctrl.integrator == "riccati" && ctrl.kinetic_factor > 0 &&
+        error("kinetic runs require integrator=\"forward\"; the Riccati integrator has no kinetic crossing.")
+
+    # The RHS works on mpert×mpert blocks, where multithreaded BLAS costs more in synchronization
+    # than it saves; pin BLAS to one thread for the sweep and restore it afterwards.
+    blas_threads = BLAS.get_num_threads()
+    BLAS.set_num_threads(1)
+    try
+        if ctrl.integrator == "riccati"
+            return riccati_eulerlagrange_integration(ctrl, equil, mats, intr)
+        else
+            return forward_eulerlagrange_integration(ctrl, equil, mats, intr)
+        end
+    finally
+        BLAS.set_num_threads(blas_threads)
     end
-    error("Unknown integrator: $(ctrl.integrator). Expected \"forward\", \"riccati\", or \"galerkin\".")
 end
 
 """
@@ -945,7 +955,7 @@ function integrate_el_region!(
 
     cb = DiscreteCallback((u, t, integrator) -> true, segment_callback!)
     prob = ODEProblem(sing_der!, odet.u, (chunk.psi_start, chunk.psi_end), (ctrl, equil, mats, intr, odet, chunk))
-    sol = solve(prob, Vern9(); reltol=ctrl.eulerlagrange_tolerance, callback=cb, save_everystep=false, save_end=true)
+    sol = solve(prob, el_ode_algorithm(ctrl); reltol=ctrl.eulerlagrange_tolerance, abstol=ctrl.ode_abstol, callback=cb, save_everystep=false, save_end=true)
 
     # Unconditionally save the final step if the callback did not already capture it.
     # Guarantees the pre-crossing (or pre-edge) state is always stored in u_store,

src/ForceFreeStates/ForceFreeStates.jl

@@ -15,6 +15,23 @@ using FastGaussQuadrature: gausslobatto
 using QuadGK: quadgk, quadgk!
 
 import ..Equilibrium
+
+"""
+Explicit Runge-Kutta methods benchmarked for the Euler-Lagrange solves; `ode_solver` must name one.
+"""
+const EL_ODE_SOLVERS = (Vern6=Vern6(), Vern7=Vern7(), Vern8=Vern8(), Vern9=Vern9(), DP8=DP8())
+
+"""
+    el_ode_algorithm(ctrl) -> OrdinaryDiffEq algorithm
+
+The solver named by `ctrl.ode_solver`, from [`EL_ODE_SOLVERS`](@ref).
+"""
+function el_ode_algorithm(ctrl)
+    name = Symbol(ctrl.ode_solver)
+    haskey(EL_ODE_SOLVERS, name) ||
+        throw(ArgumentError("ode_solver = \"$(ctrl.ode_solver)\" is not one of the benchmarked Euler-Lagrange solvers $(keys(EL_ODE_SOLVERS))"))
+    return EL_ODE_SOLVERS[name]
+end
 import ..Utilities
 import ..Vacuum
 import ..InnerLayer