Let a caller ask solve! to throw on a failed solve

CHANGELOG.md

@@ -1,5 +1,16 @@
 # Changelog
 
+## Unreleased
+
+### Added
+
+- `solve!(...; throw_on_fail=true)` throws a `SolveFailure` when the circulation
+  loop missed the solver's tolerances or the coefficients it assembled are not
+  finite — `ForwardDiff.Dual` partials included — so a caller that cannot use a
+  failed solve gets an exception rather than a `VSMSolution` it has to inspect.
+  The default is off: a post-stall solve that misses the tolerances still returns
+  its `solver_status == FAILURE` solution.
+
 ## VortexStepMethod v5.0.0 2026-09-07
 
 ### Added

docs/src/private_functions.md

@@ -19,6 +19,7 @@ local_lift_slope!
 apply_artificial_viscosity!
 frozen_wake!
 calc_forces!
+finite_full
 calculate_cl
 calculate_cd
 calculate_cm

docs/src/types.md

@@ -50,4 +50,5 @@ BodyAerodynamics
 ```@docs
 Solver
 VSMSolution
+SolveFailure
 ```

src/VortexStepMethod.jl

@@ -32,6 +32,7 @@ export slice_args, preview_args
 export ObjWing, Section, Wing, refine!, reinit!
 export BodyAerodynamics
 export Solver, VSMSolution, linearize, solve, solve!, solve_base!, calc_forces!
+export SolveFailure
 export calculate_results
 export add_section!, set_va!, section_pitch_rate
 export calculate_projected_area, calculate_span

src/solver.jl

@@ -223,9 +223,32 @@ function Solver(body_aero, settings::VSMSettings)
     )
 end
 
+"""
+    SolveFailure(msg)
+
+Thrown by `solve!(...; throw_on_fail=true)` when the circulation loop missed the
+solver's tolerances, or when the coefficients it assembled are not finite.
+"""
+struct SolveFailure <: Exception
+    msg::String
+end
+
+Base.showerror(io::IO, failure::SolveFailure) = print(io, failure.msg)
+
+"""
+    finite_full(x) -> Bool
+
+`true` if `x` is finite. A `ForwardDiff.Dual` needs every partial finite too, so
+a non-finite *derivative* of a solve is caught as well as a non-finite value.
+"""
+finite_full(x::Real) = isfinite(x)
+finite_full(x::ForwardDiff.Dual) =
+    isfinite(ForwardDiff.value(x)) && all(isfinite, ForwardDiff.partials(x))
+
 """
     solve!(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=solver.sol.gamma_distribution; 
-          log=false, reference_point=solver.reference_point, moment_frac=0.1)
+          log=false, reference_point=solver.reference_point, moment_frac=0.1,
+          throw_on_fail=false)
 
 Main solving routine for the aerodynamic model. Reference point is in the kite body (KB) frame.
 This version is modifying the `solver.sol` struct and is faster than the `solve` function which returns
@@ -240,12 +263,16 @@ a dictionary.
 - log=false: If true, print the number of iterations and other info.
 - reference_point=solver.reference_point
 - moment_frac=0.1: X-coordinate of normalized panel around which the moment distribution should be calculated.
+- throw_on_fail=false: If true, throw a [`SolveFailure`](@ref) instead of returning a
+  solution that missed the solver's tolerances or carries a non-finite coefficient.
 
 # Returns
-The solution of type [`VSMSolution`](@ref)
+The solution of type [`VSMSolution`](@ref), whose `solver_status` is `FAILURE` when the
+circulation loop missed the solver's tolerances.
 """
 function solve!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics, gamma_distribution=solver.sol.gamma_distribution;
-        log=false, reference_point=solver.reference_point, moment_frac=0.1) where {P, U, T}
+        log=false, reference_point=solver.reference_point, moment_frac=0.1,
+        throw_on_fail=false) where {P, U, T}
 
     # calculate intermediate result
     solve_base!(solver, body_aero, gamma_distribution; log)
@@ -256,7 +283,13 @@ function solve!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics, gamma_dist
         solver.sol.gamma_distribution = gamma_new
     end
 
-    return calc_forces!(solver, body_aero; reference_point, moment_frac)
+    sol = calc_forces!(solver, body_aero; reference_point, moment_frac)
+    throw_on_fail || return sol
+    solver.lr.converged || throw(SolveFailure(
+        "VSM solve did not converge in $(solver.max_iterations) iterations."))
+    all(finite_full, sol.force_coeffs) && all(finite_full, sol.moment_coeffs) ||
+        throw(SolveFailure("VSM solve assembled non-finite coefficients."))
+    return sol
 end
 
 """

test/solver/test_solver.jl

@@ -1,5 +1,6 @@
 using VortexStepMethod
 using VortexStepMethod.AirfoilAero: lei_poly_coeffs
+using ForwardDiff
 using LinearAlgebra
 using Test
 if !@isdefined(test_data_path)
@@ -197,3 +198,28 @@ end
     sol = solve!(solver_on, body_aero)
     @test all(isfinite, sol.gamma_distribution)
 end
+
+@testset "solve! reports a solve that missed the tolerances" begin
+    body_aero = BodyAerodynamics([poststall_wing])
+    solver = Solver(body_aero; solver_type=LOOP, aerodynamic_model_type=VSM,
+        max_iterations=1)
+    set_va!(body_aero, [10.0, 0.0, 0.0])
+
+    sol = solve!(solver, body_aero)
+    @test !solver.lr.converged
+    @test sol.solver_status == FAILURE
+
+    @test_throws SolveFailure solve!(solver, body_aero; throw_on_fail=true)
+    @test_throws "did not converge in 1 iterations" solve!(solver, body_aero;
+        throw_on_fail=true)
+
+    converged = Solver(body_aero; solver_type=LOOP, aerodynamic_model_type=VSM)
+    @test solve!(converged, body_aero; throw_on_fail=true) isa VSMSolution
+end
+
+@testset "finite_full sees a Dual's partials, not just its value" begin
+    @test VortexStepMethod.finite_full(1.0)
+    @test !VortexStepMethod.finite_full(NaN)
+    @test VortexStepMethod.finite_full(ForwardDiff.Dual(1.0, 2.0))
+    @test !VortexStepMethod.finite_full(ForwardDiff.Dual(1.0, Inf))
+end