diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index f107680fa1..db085067c1 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -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/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 9a2e5ad7d4..db4d2fdfa4 100644 --- a/test/runtests_solve_api.jl +++ b/test/runtests_solve_api.jl @@ -91,6 +91,9 @@ using TOML @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