Replace NLsolve.jl with NonlinearSolve.jl TrustRegion

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"

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"

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
 
 """

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

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

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)

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"

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,