diff --git a/docs/src/api.md b/docs/src/api.md index bf734f7621..1f80de680a 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -122,8 +122,8 @@ Keyword arguments of `TearingProblem` are the `[SLAYER]` deck section's procedur (scan mode and Q-domain, coupling mode, critical-Δ convention, extraction filters); kinetic profiles come from `profile_file`, relative to the run directory. -Kinetic runs (`kinetic_factor > 0`) need the `[KineticForces]` profiles and remain -TOML-driven. +Kinetic runs (`kinetic_factor > 0`) take `Forward()` and the `"fixed"` kinetic source; +the `"calculated"` source needs the `[KineticForces]` profiles and remains TOML-driven. ## Entry points diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl index 731f27f3b2..42c891e124 100644 --- a/src/ForceFreeStates/CoreTypes.jl +++ b/src/ForceFreeStates/CoreTypes.jl @@ -138,7 +138,7 @@ gpec.toml. - `rational_offset_fraction::Float64` - Distance beyond last rational surface (normalised ∈ [0,1) in units of 1/n). Only used when `truncate_at_rational_offset` is true. Fortran STRIDE convention is 0.2 (truncate 20 % of one rational-surface spacing above the last surface), retained here. - `frobenius_order::Int` - Order of singular layer (Frobenius) expansion at rational surfaces. Default 6 (Fortran STRIDE convention for Δ' calculations; lower values trade accuracy for speed). - `qhigh::Float64` - Integration terminated at q limit determined by minimum of qhigh and qa from equil - - `kinetic_source::String` - Kinetic matrix source: "fixed" (X-shaped test matrices scaled by kinetic_factor relative to ideal matrix Frobenius norms; Ak, Dk, Hk Hermitian, Bk, Ck, Ek non-Hermitian), "calculated" (the drift-kinetic matrices of the KineticForces model, evaluated on the `[KineticForces]` profiles; needs `integrator = "forward"`) + - `kinetic_source::String` - Kinetic matrix source: "fixed" (X-shaped test matrices scaled by kinetic_factor relative to ideal matrix Frobenius norms; Ak, Dk, Hk Hermitian, Bk, Ck, Ek non-Hermitian), "calculated" (the drift-kinetic matrices of the KineticForces model, evaluated on the `[KineticForces]` profiles; needs `integrator = "forward"` and is deck-driven through `main`) - `kinetic_factor::Float64` - Dimensionless scaling factor for kinetic matrices. Zero (the default) disables the kinetic path; any positive value enables it and scales the kinetic matrices: when kinetic_source="fixed", scales X-shaped test matrices relative to ideal matrix norms; when kinetic_source="calculated", applied as uniform post-hoc multiplier to W and T components. - `qlow::Float64` - Integration terminated at q limit determined by minimum of qlow and q0 from equil - `dW_edge_scan_start::Float64` - If less than psilim, records a dW(ψ) diagnostic scan over [dW_edge_scan_start, psilim] on odet.edge_scan. The integration domain (psilim) is always controlled by qhigh / psihigh and is not modified by this scan (unless `truncate_at_dW_peak=true`, see caveats below). diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index ef1d94cf4d..db085067c1 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -361,7 +361,7 @@ 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.") + ctrl.kinetic_factor > 0 && error("integrator = \"riccati\" does not support kinetic runs (kinetic_factor > 0); use integrator = \"forward\".") return riccati_eulerlagrange_integration(ctrl, equil, mats, intr) elseif ctrl.integrator == "forward" return forward_eulerlagrange_integration(ctrl, equil, mats, intr) @@ -901,8 +901,15 @@ const COLUMN_ABSTOL_CALLBACK = DiscreteCallback((u, t, integrator) -> true, refr Vern9 solve of an Euler-Lagrange problem returning the end state, with the per-column abstol of `column_abstol!`. A custom `callback` must refresh that abstol itself, after any reduction or renormalization of the state. +Errors when the solver stops short of the interval's end (e.g. the step size underflows), naming the ψ reached. """ -solve_el(prob, rtol; callback=COLUMN_ABSTOL_CALLBACK) = solve(prob, Vern9(); reltol=rtol, abstol=column_abstol(prob.u0, rtol), callback, save_everystep=false, save_end=true) +function solve_el(prob, rtol; callback=COLUMN_ABSTOL_CALLBACK) + sol = solve(prob, Vern9(); reltol=rtol, abstol=column_abstol(prob.u0, rtol), callback, save_everystep=false, save_end=true) + SciMLBase.successful_retcode(sol) || + error("Euler-Lagrange integration aborted at ψ = $(@sprintf("%.6f", sol.t[end])) before ψ = $(@sprintf("%.6f", prob.tspan[2])) " * + "(retcode $(sol.retcode)). Loosen eulerlagrange_tolerance or refine the ψ grid; for kinetic runs, lower kinetic_factor.") + return sol +end """ integrate_el_region!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, intr::ForceFreeStatesInternal, chunk::IntegrationChunk) diff --git a/src/ForceFreeStates/Riccati/Driver.jl b/src/ForceFreeStates/Riccati/Driver.jl index f3e5d4254a..941144b008 100644 --- a/src/ForceFreeStates/Riccati/Driver.jl +++ b/src/ForceFreeStates/Riccati/Driver.jl @@ -262,7 +262,7 @@ function _assemble_propagators_serially!(odet::OdeState, propagators::Vector{Chu end if chunk.needs_crossing - ctrl.kinetic_factor > 0 && error("kinetic_factor > 0 not implemented yet in Riccati!") + ctrl.kinetic_factor > 0 && error("integrator = \"riccati\" does not support kinetic runs (kinetic_factor > 0); use integrator = \"forward\".") # State is (S, I) from the renorm above — well-conditioned at the surface's left boundary. push!(S_at_surface_left, copy(odet.u[:, :, 1])) riccati_cross_ideal_singular_surf!(odet, ctrl, equil, mats, intr, chunk.ising) diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 6af667a0a4..f864604cb3 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -849,7 +849,7 @@ end dir_path=".", debug=DebugSettings(), kwargs...) The perturbed-plasma Euler-Lagrange problem posed on an equilibrium: the extremization of -the perturbed potential energy whose solutions are the force-free (and, via the TOML path, +the perturbed potential energy whose solutions are the force-free (or, with `kinetic_factor > 0`, kinetic) perturbed states. This is the WHAT of a stability solve; the integrator passed to [`solve`](@ref) is the HOW. A `PlasmaEquilibrium` hosts many possible problems — this type names this one, so `solve` stays unambiguous as other problem classes appear. Inner-layer @@ -903,9 +903,9 @@ Solve the perturbed-plasma [`EulerLagrangeProblem`](@ref) with the formalism `al a `gpec.toml` run of `main` does and produces the same result object. The second form is sugar building the problem from an equilibrium and the problem keywords in one call. -Knobs owned by `alg` are rejected as `ForceFreeStatesControl` keywords. Kinetic -runs are TOML-driven this cycle: `kinetic_factor > 0` needs the `[KineticForces]` profiles -and errors here. +Knobs owned by `alg` are rejected as `ForceFreeStatesControl` keywords. Kinetic runs +(`kinetic_factor > 0`) take `Forward()` and the `"fixed"` kinetic source; the `"calculated"` +source needs the `[KineticForces]` profiles and stays TOML-driven this cycle. ```julia eq = PlasmaEquilibrium("input.geqdsk"; jac_type="hamada") @@ -922,8 +922,10 @@ function solve(prob::EulerLagrangeProblem, alg::ForceFreeStates.AbstractIntegrat ForceFreeStates._apply_alg!(ctrl_kwargs, alg) ctrl = ForceFreeStatesControl(; ctrl_kwargs...) - ctrl.kinetic_factor > 0 && - error("kinetic runs (kinetic_factor > 0) need the [KineticForces] profiles and are TOML-driven; run them through `main`") + ctrl.kinetic_factor > 0 && ctrl.kinetic_source == "calculated" && + error("kinetic_source = \"calculated\" needs the [KineticForces] profiles and is TOML-driven; run it through `main`") + ctrl.kinetic_factor > 0 && !(alg isa ForceFreeStates.Forward) && + error("$(nameof(typeof(alg)))() does not support kinetic runs (kinetic_factor > 0); use Forward().") intr = ForceFreeStatesInternal(; dir_path=prob.dir_path) intr.wall_settings = prob.wall diff --git a/test/runtests_eulerlagrange.jl b/test/runtests_eulerlagrange.jl index e160032af2..2a693a1e12 100644 --- a/test/runtests_eulerlagrange.jl +++ b/test/runtests_eulerlagrange.jl @@ -640,3 +640,12 @@ end @test !odet.du_store_populated end end + +@testset "solve_el reports an aborted integration" begin + FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates + # du/dψ = u² blows up at ψ = 1, so the solver stops short of the interval's end. + blowup = FFS.ODEProblem((du, u, p, t) -> (du .= u .^ 2), ones(ComplexF64, 1, 1, 2), (0.0, 2.0)) + @test_throws "Euler-Lagrange integration aborted" FFS.solve_el(blowup, 1e-8) + calm = FFS.ODEProblem((du, u, p, t) -> (du .= -u), ones(ComplexF64, 1, 1, 2), (0.0, 1.0)) + @test FFS.solve_el(calm, 1e-8).u[end] ≈ fill(exp(-1), 1, 1, 2) rtol = 1e-6 +end diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl index bc77acf044..6e8d2aa393 100644 --- a/test/runtests_solve_api.jl +++ b/test/runtests_solve_api.jl @@ -76,6 +76,27 @@ using TOML end end + @testset "fixed kinetic solve matches the TOML-driven run" begin + # Large enough that the fixed kinetic matrices visibly shift the free-boundary energies. + kf = 1e-3 + mktempdir() do dir + reference = GPEC.main([_stage_deck(dir; ffs_overrides=Dict{String,Any}("kinetic_factor" => kf))]).ffs + api = solve(equil, Forward(); nn=1, wall=wall, dir_path=dir, ffs_kwargs..., kinetic_factor=kf, kinetic_source="fixed") + + @test api.mats.kinetic !== nothing + @test api.free_boundary.et[1] ≈ reference.free_boundary.et[1] rtol = 1e-12 + @test api.diagnostics.nzero == reference.diagnostics.nzero + @test api.mlow == reference.mlow + @test api.mhigh == reference.mhigh + @test api.psilim ≈ reference.psilim rtol = 1e-12 + @test api.kinetic.kmsing == reference.kinetic.kmsing + @test [s.psifac for s in api.kinetic.kinsing] ≈ [s.psifac for s in reference.kinetic.kinsing] + + # Too large a factor for this coarse deck: the integration aborts near ψ = 0.4 and says so. + @test_throws "Euler-Lagrange integration aborted" solve(equil, Forward(); nn=1, wall=wall, dir_path=dir, ffs_kwargs..., kinetic_factor=0.1) + end + end + @testset "galerkin solve returns a gal result with a unified delta_prime" begin mktempdir() do dir api = solve(equil, Galerkin(; nx=32); nn=1, wall=wall, dir_path=dir, ffs_kwargs...) @@ -186,7 +207,9 @@ using TOML end @testset "rejected keyword combinations" begin - @test_throws ErrorException solve(equil, Forward(); nn=1, dir_path=".", ffs_kwargs..., kinetic_factor=0.5) + @test_throws "kinetic_source = \"calculated\"" solve(equil, Forward(); nn=1, dir_path=".", ffs_kwargs..., kinetic_factor=0.5, kinetic_source="calculated") + @test_throws "does not support kinetic runs" solve(equil, Riccati(); nn=1, dir_path=".", ffs_kwargs..., kinetic_factor=1e-9) + @test_throws "does not support kinetic runs" solve(equil, Galerkin(); nn=1, dir_path=".", ffs_kwargs..., kinetic_factor=1e-9) @test_throws ErrorException solve(equil, Forward(); nn=1, dir_path=".", ffs_kwargs..., integrator="riccati") @test_throws ErrorException solve(equil, Riccati(); nn=1, dir_path=".", ffs_kwargs..., nchunks=8) @test_throws ErrorException solve(equil, Forward(); nn=1, dir_path=".", ffs_kwargs..., nn_low=2)