From eb90218c44e72ef07540b1a1a0a2c53f2ae1878e Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Thu, 19 Feb 2026 04:54:36 -0500 Subject: [PATCH] Replace NLsolve.jl with NonlinearSolve.jl TrustRegion for steady-state solving MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit NLsolve.jl's trust region dogleg solver is replaced by NonlinearSolve.jl's TrustRegion() in both KPS3 and KPS4 find_steady_state! functions. This removes NLsolve and NLSolversBase as direct dependencies while NonlinearSolve was already a transitive dependency via OrdinaryDiffEqNonlinearSolve. The custom finite-difference Jacobian (make_jac) is retained because the model's MVector{3, Float64} internals cannot hold ForwardDiff.Dual numbers. NLsolve's autoscale=true is not needed — TrustRegion works without it and actually produces slightly better results on the l=392 test case. All 452 tests pass (441 pass + 11 pre-existing broken). Co-Authored-By: Chris Rackauckas --- Project.toml | 6 ++---- examples_3d/Project.toml | 2 -- src/KPS3.jl | 21 ++++++++++++++------- src/KPS4.jl | 21 ++++++++++++++------- src/KiteModels.jl | 3 ++- src/init.jl | 10 +++++++--- test/Project.toml | 1 - test/create_sys_image.jl | 4 ++-- 8 files changed, 41 insertions(+), 27 deletions(-) diff --git a/Project.toml b/Project.toml index 7a6891590..70020ccbf 100644 --- a/Project.toml +++ b/Project.toml @@ -18,9 +18,8 @@ KiteUtils = "90980105-b163-44e5-ba9f-8b1c83bb0533" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" -NLSolversBase = "d41bc354-129a-5804-8e4c-c37616107c6c" -NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" +SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" OrdinaryDiffEqBDF = "6ad6398a-0878-4a85-9266-38940aa047c8" OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8" @@ -48,9 +47,8 @@ KiteUtils = "0.11.1" LinearAlgebra = "1.10, 1.11, 1.12" LinearSolve = "~2.39.0, 3" Logging = "1.11.0" -NLSolversBase = "7.10.0" -NLsolve = "4.5" NonlinearSolve = "4.11" +SciMLBase = "2" OrdinaryDiffEqBDF = "1.5.0" OrdinaryDiffEqCore = "1.23.0" OrdinaryDiffEqNonlinearSolve = "1.6.1" diff --git a/examples_3d/Project.toml b/examples_3d/Project.toml index cd6e1c9be..0d7d0793f 100644 --- a/examples_3d/Project.toml +++ b/examples_3d/Project.toml @@ -11,7 +11,6 @@ KiteUtils = "90980105-b163-44e5-ba9f-8b1c83bb0533" KiteViewers = "2e593061-95e7-45e4-95f4-df0491f2e601" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" -NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" OrdinaryDiffEqBDF = "6ad6398a-0878-4a85-9266-38940aa047c8" OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf" @@ -40,7 +39,6 @@ KitePodModels = "0.3.4" KiteUtils = "0.8.1" KiteViewers = "0.4.18" ModelingToolkit = "~9.42.0" -NLsolve = "4.5" OrdinaryDiffEqBDF = "1.1.1" OrdinaryDiffEqCore = "1.4.0" OrdinaryDiffEqSDIRK = "1.1.0" diff --git a/src/KPS3.jl b/src/KPS3.jl index 827e2495d..a0813a6fe 100644 --- a/src/KPS3.jl +++ b/src/KPS3.jl @@ -510,14 +510,21 @@ function find_steady_state_inner(s::KPS3, X, prn=false; delta=0.0) return nothing end - if prn println("\nStarted function test_nlsolve...") end - jac! = make_jac(test_initial_condition!, length(X)) - results = nlsolve(test_initial_condition!, jac!, X, xtol=1e-6, ftol=1e-6, autoscale=true, iterations=1000) - if prn println("\nresult: $results") end - if !converged(results) - @warn "find_steady_state!: solver did not converge! (f_converged=$(results.f_converged), x_converged=$(results.x_converged), iterations=$(results.iterations))" + if prn println("\nStarted NonlinearSolve...") end + jac_2arg! = make_jac(test_initial_condition!, length(X)) + + # Wrap for NonlinearSolve's 3-arg (out, u, p) convention + f_nl!(F, x, p) = test_initial_condition!(F, x) + jac_nl!(J, x, p) = jac_2arg!(J, x) + + nf = NonlinearFunction(f_nl!; jac=jac_nl!) + prob = NonlinearProblem(nf, X) + sol = solve(prob, TrustRegion(); abstol=1e-6, reltol=1e-6, maxiters=1000) + if prn println("\nresult retcode: $(sol.retcode)") end + if !SciMLBase.successful_retcode(sol.retcode) + @warn "find_steady_state!: solver did not converge! retcode=$(sol.retcode)" end - results.zero + sol.u end """ diff --git a/src/KPS4.jl b/src/KPS4.jl index be12463f5..23800ad97 100644 --- a/src/KPS4.jl +++ b/src/KPS4.jl @@ -688,15 +688,22 @@ function find_steady_state!(s::KPS4; prn=false, delta = 0.001, stiffness_factor= iter += 1 return nothing end - if prn println("\nStarted function test_nlsolve...") end + if prn println("\nStarted NonlinearSolve...") end X00 = zeros(SimFloat, 2*(s.set.segments+KITE_PARTICLES-1)+2) - jac! = make_jac(test_initial_condition!, length(X00)) - results = nlsolve(test_initial_condition!, jac!, X00, autoscale=true, xtol=4e-7, ftol=4e-7, iterations=s.set.max_iter) - if prn println("\nresult: $results") end - if !converged(results) - @warn "find_steady_state!: solver did not converge! (f_converged=$(results.f_converged), x_converged=$(results.x_converged), iterations=$(results.iterations))" + jac_2arg! = make_jac(test_initial_condition!, length(X00)) + + # Wrap for NonlinearSolve's 3-arg (out, u, p) convention + f_nl!(F, x, p) = test_initial_condition!(F, x) + jac_nl!(J, x, p) = jac_2arg!(J, x) + + nf = NonlinearFunction(f_nl!; jac=jac_nl!) + prob = NonlinearProblem(nf, X00) + sol = solve(prob, TrustRegion(); abstol=4e-7, reltol=4e-7, maxiters=s.set.max_iter) + if prn println("\nresult retcode: $(sol.retcode)") end + if !SciMLBase.successful_retcode(sol.retcode) + @warn "find_steady_state!: solver did not converge! retcode=$(sol.retcode)" end - y0, yd0 = init(s, results.zero; upwind_dir) + y0, yd0 = init(s, sol.u; upwind_dir) set_v_wind_ground!(s, calc_height(s), s.set.v_wind; upwind_dir) residual!(res, yd0, y0, s, 0.0) y0, yd0 diff --git a/src/KiteModels.jl b/src/KiteModels.jl index 206cbbdf8..57b40dbee 100644 --- a/src/KiteModels.jl +++ b/src/KiteModels.jl @@ -14,9 +14,10 @@ module KiteModels using PrecompileTools: @setup_workload, @compile_workload using Logging -using Dierckx, Interpolations, StaticArrays, LinearAlgebra, Statistics, Parameters, NLsolve, +using Dierckx, Interpolations, StaticArrays, LinearAlgebra, Statistics, Parameters, DocStringExtensions, OrdinaryDiffEqCore, OrdinaryDiffEqBDF, OrdinaryDiffEqNonlinearSolve, NonlinearSolve, Suppressor +import SciMLBase import Sundials using Reexport, Pkg using KiteUtils diff --git a/src/init.jl b/src/init.jl index 2440f0d00..5c4b56a1d 100644 --- a/src/init.jl +++ b/src/init.jl @@ -9,9 +9,13 @@ const PRE_STRESS = 0.9998 # Multiplier for the initial spring lengths. Create an explicit Jacobian function using central finite differences. -Works around NLSolversBase >= 7.10 producing subnormal Jacobian entries via -DifferentiationInterface, which cause NaN in NLsolve autoscale (subnormal -column norms squared underflow to zero). +Returns a two-argument `jac!(J, x)` closure. Callers wrap it for +NonlinearSolve's three-argument `jac!(J, x, p)` convention. + +The custom Jacobian is needed because the model's `MVector{3, Float64}` +internals cannot hold `ForwardDiff.Dual` numbers, and NonlinearSolve's +built-in `AutoFiniteDiff` uses different step sizing without subnormal +flushing, which causes convergence failures. """ function make_jac(f!, n_vars) _F1 = zeros(SimFloat, n_vars) diff --git a/test/Project.toml b/test/Project.toml index 5d1616b57..a74ef4845 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -16,7 +16,6 @@ KitePodModels = "9de5dc81-f971-414a-927b-652b2f41c539" KiteUtils = "90980105-b163-44e5-ba9f-8b1c83bb0533" LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" OrdinaryDiffEqBDF = "6ad6398a-0878-4a85-9266-38940aa047c8" OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" diff --git a/test/create_sys_image.jl b/test/create_sys_image.jl index b0f503c3b..246671400 100644 --- a/test/create_sys_image.jl +++ b/test/create_sys_image.jl @@ -9,7 +9,7 @@ if ! ("PackageCompiler" ∈ keys(Pkg.project().dependencies)) end @info "Loading packages ..." using AtmosphericModels, Colors, ControlSystemsBase, DSP, Dierckx, DiscretePIDs, - DocStringExtensions, JLD2, KitePodModels, KiteUtils, LinearAlgebra, NLsolve, + DocStringExtensions, JLD2, KitePodModels, KiteUtils, LinearAlgebra, NonlinearSolve, OrdinaryDiffEqBDF, OrdinaryDiffEqCore, OrdinaryDiffEqNonlinearSolve, Parameters, REPL, StaticArrays, StatsBase, Sundials, WinchModels using BenchmarkTools, Documenter, PackageCompiler @@ -18,7 +18,7 @@ using BenchmarkTools, Documenter, PackageCompiler push!(LOAD_PATH,joinpath(pwd(),"src")) PackageCompiler.create_sysimage( - [:Dierckx, :StaticArrays, :Parameters, :NLsolve, :DocStringExtensions, :Sundials, :KiteUtils, + [:Dierckx, :StaticArrays, :Parameters, :DocStringExtensions, :Sundials, :KiteUtils, :KitePodModels, :AtmosphericModels, :OrdinaryDiffEqCore, :OrdinaryDiffEqBDF, :WinchModels, :OrdinaryDiffEqNonlinearSolve, :StatsBase, :PackageCompiler, :BenchmarkTools, :Documenter, :JLD2, :Colors, :REPL, :NonlinearSolve, :DSP, :DiscretePIDs,