Merge pull request #240 from SciML/default_cases

Change default in-place specialization
SciML · Oct 17, 2023 · 08f7f1d · 08f7f1d
2 parents d327107 + 143c634
commit 08f7f1d
Show file tree

Hide file tree

Showing 2 changed files with 86 additions and 4 deletions.
diff --git a/src/default.jl b/src/default.jl
@@ -138,8 +138,8 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::FastShortcutN
     )
 end
 
-function SciMLBase.__solve(prob::NonlinearProblem{uType, iip}, alg::FastShortcutNonlinearPolyalg, args...; 
-    kwargs...) where {uType, iip}
+function SciMLBase.__solve(prob::NonlinearProblem{uType, false}, alg::FastShortcutNonlinearPolyalg, args...; 
+    kwargs...) where {uType}
 
     adkwargs = alg.adkwargs
     linsolve = alg.linsolve
@@ -190,6 +190,58 @@ function SciMLBase.__solve(prob::NonlinearProblem{uType, iip}, alg::FastShortcut
 
 end
 
+function SciMLBase.__solve(prob::NonlinearProblem{uType, true}, alg::FastShortcutNonlinearPolyalg, args...; 
+    kwargs...) where {uType}
+
+    adkwargs = alg.adkwargs
+    linsolve = alg.linsolve
+    precs = alg.precs
+
+    sol1 = SciMLBase.__solve(prob, NewtonRaphson(;linsolve, precs, adkwargs...), args...; kwargs...)
+    if SciMLBase.successful_retcode(sol1)
+        return SciMLBase.build_solution(prob, alg, sol1.u, sol1.resid;
+                                        sol1.retcode, sol1.stats)
+    end
+
+    sol2 = SciMLBase.__solve(prob, NewtonRaphson(;linsolve, precs, linesearch=BackTracking(), adkwargs...), args...; kwargs...)
+    if SciMLBase.successful_retcode(sol2)
+        return SciMLBase.build_solution(prob, alg, sol2.u, sol2.resid;
+                                        sol2.retcode, sol2.stats)
+    end
+
+    sol3 = SciMLBase.__solve(prob, TrustRegion(;linsolve, precs, adkwargs...), args...; kwargs...)
+    if SciMLBase.successful_retcode(sol3)
+        return SciMLBase.build_solution(prob, alg, sol3.u, sol3.resid;
+                                        sol3.retcode, sol3.stats)
+    end
+
+    sol4 = SciMLBase.__solve(prob,  TrustRegion(;linsolve, precs, radius_update_scheme = RadiusUpdateSchemes.Bastin, adkwargs...), args...; kwargs...)
+    if SciMLBase.successful_retcode(sol4)
+        return SciMLBase.build_solution(prob, alg, sol4.u, sol4.resid;
+                                        sol4.retcode, sol4.stats)
+    end
+
+    resids = (sol1.resid, sol2.resid, sol3.resid, sol4.resid)
+    minfu, idx = findmin(DEFAULT_NORM, resids)
+
+    if idx == 1
+        SciMLBase.build_solution(prob, alg, sol1.u, sol1.resid;
+                                        sol1.retcode, sol1.stats)
+    elseif idx == 2
+        SciMLBase.build_solution(prob, alg, sol2.u, sol2.resid;
+                                        sol2.retcode, sol2.stats)
+    elseif idx == 3
+        SciMLBase.build_solution(prob, alg, sol3.u, sol3.resid;
+                                sol3.retcode, sol3.stats)
+    elseif idx == 4
+        SciMLBase.build_solution(prob, alg, sol4.u, sol4.resid;
+                                sol4.retcode, sol4.stats)
+    else
+        error("Unreachable reached, 박정석")
+    end
+
+end
+
 ## General shared polyalg functions
 
 function perform_step!(cache::Union{RobustMultiNewtonCache, FastShortcutNonlinearPolyalgCache})

diff --git a/test/polyalgs.jl b/test/polyalgs.jl
@@ -1,8 +1,38 @@
-using NonlinearSolve
+using NonlinearSolve, Test
 
 f(u, p) = u .* u .- 2
 u0 = [1.0, 1.0]
 probN = NonlinearProblem(f, u0)
 @time solver = solve(probN, abstol = 1e-9)
 @time solver = solve(probN, RobustMultiNewton(), abstol = 1e-9)
-@time solver = solve(probN, FastShortcutNonlinearPolyalg(), abstol = 1e-9)
+@time solver = solve(probN, FastShortcutNonlinearPolyalg(), abstol = 1e-9)
+
+# https://github.com/SciML/NonlinearSolve.jl/issues/153
+
+function f(du, u, p)
+    s1, s1s2, s2 = u
+    k1, c1, Δt = p
+
+    du[1] = -0.25 * c1 * k1 * s1 * s2
+    du[2] = 0.25 * c1 * k1 * s1 * s2
+    du[3] = -0.25 * c1 * k1 * s1 * s2
+end
+
+prob = NonlinearProblem(f, [2.0,2.0,2.0], [1.0, 2.0, 2.5])
+sol = solve(prob)
+@test SciMLBase.successful_retcode(sol)
+
+# https://github.com/SciML/NonlinearSolve.jl/issues/187
+
+ff(u, p) = 0.5/1.5*log.(u./(1.0.-u)) .- 2.0*u .+1.0
+
+uspan = (0.02, 0.1)
+prob = IntervalNonlinearProblem(ff, uspan)
+sol = solve(prob)
+@test SciMLBase.successful_retcode(sol)
+
+u0 = 0.06
+p = 2.0
+prob = NonlinearProblem(ff, u0, p)
+solver = solve(prob)
+@test SciMLBase.successful_retcode(sol)