diff --git a/CHANGELOG.md b/CHANGELOG.md index 6bd63200..4ae8f69f 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index 9658d880..cd8f3db7 100644 --- a/docs/src/private_functions.md +++ b/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 diff --git a/docs/src/types.md b/docs/src/types.md index cb1c9388..21a4d425 100644 --- a/docs/src/types.md +++ b/docs/src/types.md @@ -50,4 +50,5 @@ BodyAerodynamics ```@docs Solver VSMSolution +SolveFailure ``` diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl index bef3f873..0521ddfc 100644 --- a/src/VortexStepMethod.jl +++ b/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 diff --git a/src/solver.jl b/src/solver.jl index c0e35282..b52161b6 100644 --- a/src/solver.jl +++ b/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 """ diff --git a/test/solver/test_solver.jl b/test/solver/test_solver.jl index 338587c0..1dfbd78a 100644 --- a/test/solver/test_solver.jl +++ b/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