Merge pull request #264 from avik-pal/ap/patch

Don't change the default termination condition
SciML · Oct 27, 2023 · 3565824 · 3565824
2 parents 8fe131a + c428bd9
commit 3565824
Show file tree

Hide file tree

Showing 11 changed files with 95 additions and 161 deletions.
diff --git a/src/broyden.jl b/src/broyden.jl
@@ -185,9 +185,8 @@ function SciMLBase.reinit!(cache::GeneralBroydenCache{iip}, u0 = cache.u; p = ca
         cache.u = u0
         cache.fu = cache.f(cache.u, p)
     end
-    termination_condition = _get_reinit_termination_condition(cache,
-        abstol,
-        reltol,
+
+    termination_condition = _get_reinit_termination_condition(cache, abstol, reltol,
         termination_condition)
 
     cache.abstol = abstol

diff --git a/src/dfsane.jl b/src/dfsane.jl
@@ -114,23 +114,18 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::DFSane, args.
     𝒹, uₙ₋₁, fuₙ, fuₙ₋₁ = copy(uₙ), copy(uₙ), copy(uₙ), copy(uₙ)
 
     if iip
-        # f = (dx, x) -> prob.f(dx, x, p)
-        # f(fuₙ₋₁, uₙ₋₁)
         prob.f(fuₙ₋₁, uₙ₋₁, p)
     else
-        # f = (x) -> prob.f(x, p)
-        fuₙ₋₁ = prob.f(uₙ₋₁, p) # f(uₙ₋₁)
+        fuₙ₋₁ = prob.f(uₙ₋₁, p)
     end
 
     f₍ₙₒᵣₘ₎ₙ₋₁ = norm(fuₙ₋₁)^nₑₓₚ
     f₍ₙₒᵣₘ₎₀ = f₍ₙₒᵣₘ₎ₙ₋₁
 
     ℋ = fill(f₍ₙₒᵣₘ₎ₙ₋₁, M)
 
-    abstol, reltol, termination_condition = _init_termination_elements(abstol,
-        reltol,
-        termination_condition,
-        T)
+    abstol, reltol, termination_condition = _init_termination_elements(abstol, reltol,
+        termination_condition, T)
 
     mode = DiffEqBase.get_termination_mode(termination_condition)
 
@@ -167,14 +162,13 @@ function perform_step!(cache::DFSaneCache{true})
 
     f(cache.fuₙ, cache.uₙ)
     f₍ₙₒᵣₘ₎ₙ = norm(cache.fuₙ)^nₑₓₚ
-    for _ in 1:(cache.alg.max_inner_iterations)
+    for jjj in 1:(cache.alg.max_inner_iterations)
         𝒸 = f̄ + η - γ * α₊^2 * f₍ₙₒᵣₘ₎ₙ₋₁
 
         f₍ₙₒᵣₘ₎ₙ ≤ 𝒸 && break
 
         α₊ = clamp(α₊^2 * f₍ₙₒᵣₘ₎ₙ₋₁ / (f₍ₙₒᵣₘ₎ₙ + (T(2) * α₊ - T(1)) * f₍ₙₒᵣₘ₎ₙ₋₁),
-            τₘᵢₙ * α₊,
-            τₘₐₓ * α₊)
+            τₘᵢₙ * α₊, τₘₐₓ * α₊)
         @. cache.uₙ = cache.uₙ₋₁ - α₋ * cache.𝒹
 
         f(cache.fuₙ, cache.uₙ)
@@ -183,8 +177,7 @@ function perform_step!(cache::DFSaneCache{true})
         f₍ₙₒᵣₘ₎ₙ .≤ 𝒸 && break
 
         α₋ = clamp(α₋^2 * f₍ₙₒᵣₘ₎ₙ₋₁ / (f₍ₙₒᵣₘ₎ₙ + (T(2) * α₋ - T(1)) * f₍ₙₒᵣₘ₎ₙ₋₁),
-            τₘᵢₙ * α₋,
-            τₘₐₓ * α₋)
+            τₘᵢₙ * α₋, τₘₐₓ * α₋)
 
         @. cache.uₙ = cache.uₙ₋₁ + α₊ * cache.𝒹
         f(cache.fuₙ, cache.uₙ)
@@ -207,7 +200,7 @@ function perform_step!(cache::DFSaneCache{true})
     # Spectral parameter bounds check
     if abs(cache.σₙ) > σₘₐₓ || abs(cache.σₙ) < σₘᵢₙ
         test_norm = sqrt(sum(abs2, cache.fuₙ₋₁))
-        cache.σₙ = clamp(1.0 / test_norm, 1, 1e5)
+        cache.σₙ = clamp(T(1) / test_norm, T(1), T(1e5))
     end
 
     # Take step
@@ -283,7 +276,7 @@ function perform_step!(cache::DFSaneCache{false})
     # Spectral parameter bounds check
     if abs(cache.σₙ) > σₘₐₓ || abs(cache.σₙ) < σₘᵢₙ
         test_norm = sqrt(sum(abs2, cache.fuₙ₋₁))
-        cache.σₙ = clamp(1.0 / test_norm, 1, 1e5)
+        cache.σₙ = clamp(T(1) / test_norm, T(1), T(1e5))
     end
 
     # Take step
@@ -337,9 +330,7 @@ function SciMLBase.reinit!(cache::DFSaneCache{iip}, u0 = cache.uₙ; p = cache.p
     T = eltype(cache.uₙ)
     cache.σₙ = T(cache.alg.σ_1)
 
-    termination_condition = _get_reinit_termination_condition(cache,
-        abstol,
-        reltol,
+    termination_condition = _get_reinit_termination_condition(cache, abstol, reltol,
         termination_condition)
 
     cache.abstol = abstol

diff --git a/src/gaussnewton.jl b/src/gaussnewton.jl
@@ -109,21 +109,17 @@ function SciMLBase.__init(prob::NonlinearLeastSquaresProblem{uType, iip}, alg_::
         JᵀJ, Jᵀf = nothing, nothing
     end
 
-    abstol, reltol, termination_condition = _init_termination_elements(abstol,
-        reltol,
-        termination_condition,
-        eltype(u); mode = NLSolveTerminationMode.AbsNorm)
+    abstol, reltol, termination_condition = _init_termination_elements(abstol, reltol,
+        termination_condition, eltype(u); mode = NLSolveTerminationMode.AbsNorm)
 
     mode = DiffEqBase.get_termination_mode(termination_condition)
 
     storage = mode ∈ DiffEqBase.SAFE_TERMINATION_MODES ? NLSolveSafeTerminationResult() :
               nothing
 
     return GaussNewtonCache{iip}(f, alg, u, copy(u), fu1, fu2, zero(fu1), du, p, uf,
-        linsolve, J,
-        JᵀJ, Jᵀf, jac_cache, false, maxiters, internalnorm, ReturnCode.Default, abstol,
-        reltol,
-        prob, NLStats(1, 0, 0, 0, 0), storage, termination_condition)
+        linsolve, J, JᵀJ, Jᵀf, jac_cache, false, maxiters, internalnorm, ReturnCode.Default,
+        abstol, reltol, prob, NLStats(1, 0, 0, 0, 0), storage, termination_condition)
 end
 
 function perform_step!(cache::GaussNewtonCache{true})
@@ -149,10 +145,7 @@ function perform_step!(cache::GaussNewtonCache{true})
     @. u = u - du
     f(cache.fu_new, u, p)
 
-    (termination_condition(cache.fu_new .- cache.fu1,
-        cache.u,
-        u_prev,
-        cache.abstol,
+    (termination_condition(cache.fu_new .- cache.fu1, cache.u, u_prev, cache.abstol,
         cache.reltol) ||
      termination_condition(cache.fu_new, cache.u, u_prev, cache.abstol, cache.reltol)) &&
         (cache.force_stop = true)
@@ -219,9 +212,8 @@ function SciMLBase.reinit!(cache::GaussNewtonCache{iip}, u0 = cache.u; p = cache
         cache.u = u0
         cache.fu1 = cache.f(cache.u, p)
     end
-    termination_condition = _get_reinit_termination_condition(cache,
-        abstol,
-        reltol,
+
+    termination_condition = _get_reinit_termination_condition(cache, abstol, reltol,
         termination_condition)
 
     cache.abstol = abstol

diff --git a/src/klement.jl b/src/klement.jl
@@ -238,10 +238,9 @@ function SciMLBase.reinit!(cache::GeneralKlementCache{iip}, u0 = cache.u; p = ca
         cache.fu = cache.f(cache.u, p)
     end
 
-    termination_condition = _get_reinit_termination_condition(cache,
-        abstol,
-        reltol,
+    termination_condition = _get_reinit_termination_condition(cache, abstol, reltol,
         termination_condition)
+
     cache.abstol = abstol
     cache.reltol = reltol
     cache.termination_condition = termination_condition

diff --git a/src/levenberg.jl b/src/levenberg.jl
@@ -185,10 +185,8 @@ function SciMLBase.__init(prob::Union{NonlinearProblem{uType, iip},
         v = similar(du)
     end
 
-    abstol, reltol, termination_condition = _init_termination_elements(abstol,
-        reltol,
-        termination_condition,
-        eltype(u); mode = NLSolveTerminationMode.AbsNorm)
+    abstol, reltol, termination_condition = _init_termination_elements(abstol, reltol,
+        termination_condition, eltype(u); mode = NLSolveTerminationMode.AbsNorm)
 
     λ = convert(eltype(u), alg.damping_initial)
     λ_factor = convert(eltype(u), alg.damping_increase_factor)

diff --git a/src/pseudotransient.jl b/src/pseudotransient.jl
@@ -2,10 +2,13 @@
     PseudoTransient(; concrete_jac = nothing, linsolve = nothing,
         precs = DEFAULT_PRECS, alpha_initial = 1e-3, adkwargs...)
 
-An implementation of PseudoTransient method that is used to solve steady state problems in an accelerated manner. It uses an adaptive time-stepping to
-integrate an initial value of nonlinear problem until sufficient accuracy in the desired steady-state is achieved to switch over to Newton's method and
-gain a rapid convergence. This implementation specifically uses "switched evolution relaxation" SER method. For detail information about the time-stepping and algorithm,
-please see the paper: [Coffey, Todd S. and Kelley, C. T. and Keyes, David E. (2003), Pseudotransient Continuation and Differential-Algebraic Equations,
+An implementation of PseudoTransient method that is used to solve steady state problems in
+an accelerated manner. It uses an adaptive time-stepping to integrate an initial value of
+nonlinear problem until sufficient accuracy in the desired steady-state is achieved to
+switch over to Newton's method and gain a rapid convergence. This implementation
+specifically uses "switched evolution relaxation" SER method. For detail information about
+the time-stepping and algorithm, please see the paper:
+[Coffey, Todd S. and Kelley, C. T. and Keyes, David E. (2003), Pseudotransient Continuation and Differential-Algebraic Equations,
 SIAM Journal on Scientific Computing,25, 553-569.](https://doi.org/10.1137/S106482750241044X)
 
 ### Keyword Arguments
@@ -48,7 +51,7 @@ function PseudoTransient(; concrete_jac = nothing, linsolve = nothing,
     return PseudoTransient{_unwrap_val(concrete_jac)}(ad, linsolve, precs, alpha_initial)
 end
 
-@concrete mutable struct PseudoTransientCache{iip}
+@concrete mutable struct PseudoTransientCache{iip} <: AbstractNonlinearSolveCache{iip}
     f
     alg
     u
@@ -75,14 +78,10 @@ end
     tc_storage
 end
 
-isinplace(::PseudoTransientCache{iip}) where {iip} = iip
-
 function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg_::PseudoTransient,
-    args...;
-    alias_u0 = false, maxiters = 1000, abstol = nothing, reltol = nothing,
+    args...; alias_u0 = false, maxiters = 1000, abstol = nothing, reltol = nothing,
     termination_condition = nothing, internalnorm = DEFAULT_NORM,
-    linsolve_kwargs = (;),
-    kwargs...) where {uType, iip}
+    linsolve_kwargs = (;), kwargs...) where {uType, iip}
     alg = get_concrete_algorithm(alg_, prob)
 
     @unpack f, u0, p = prob
@@ -99,9 +98,7 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg_::PseudoTransi
     res_norm = internalnorm(fu1)
 
     abstol, reltol, termination_condition = _init_termination_elements(abstol,
-        reltol,
-        termination_condition,
-        eltype(u))
+        reltol, termination_condition, eltype(u))
 
     mode = DiffEqBase.get_termination_mode(termination_condition)
 
@@ -111,8 +108,7 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg_::PseudoTransi
     return PseudoTransientCache{iip}(f, alg, u, copy(u), fu1, fu2, du, p, alpha, res_norm,
         uf,
         linsolve, J, jac_cache, false, maxiters, internalnorm, ReturnCode.Default, abstol,
-        reltol,
-        prob, NLStats(1, 0, 0, 0, 0), termination_condition, storage)
+        reltol, prob, NLStats(1, 0, 0, 0, 0), termination_condition, storage)
 end
 
 function perform_step!(cache::PseudoTransientCache{true})
@@ -176,22 +172,6 @@ function perform_step!(cache::PseudoTransientCache{false})
     return nothing
 end
 
-function SciMLBase.solve!(cache::PseudoTransientCache)
-    while !cache.force_stop && cache.stats.nsteps < cache.maxiters
-        perform_step!(cache)
-        cache.stats.nsteps += 1
-    end
-
-    if cache.stats.nsteps == cache.maxiters
-        cache.retcode = ReturnCode.MaxIters
-    else
-        cache.retcode = ReturnCode.Success
-    end
-
-    return SciMLBase.build_solution(cache.prob, cache.alg, cache.u, cache.fu1;
-        cache.retcode, cache.stats)
-end
-
 function SciMLBase.reinit!(cache::PseudoTransientCache{iip}, u0 = cache.u; p = cache.p,
     alpha_new,
     abstol = cache.abstol, reltol = cache.reltol,
@@ -207,9 +187,7 @@ function SciMLBase.reinit!(cache::PseudoTransientCache{iip}, u0 = cache.u; p = c
         cache.fu1 = cache.f(cache.u, p)
     end
 
-    termination_condition = _get_reinit_termination_condition(cache,
-        abstol,
-        reltol,
+    termination_condition = _get_reinit_termination_condition(cache, abstol, reltol,
         termination_condition)
 
     cache.alpha = convert(eltype(cache.u), alpha_new)

diff --git a/src/raphson.jl b/src/raphson.jl
@@ -179,10 +179,9 @@ function SciMLBase.reinit!(cache::NewtonRaphsonCache{iip}, u0 = cache.u; p = cac
         cache.fu1 = cache.f(cache.u, p)
     end
 
-    termination_condition = _get_reinit_termination_condition(cache,
-        abstol,
-        reltol,
+    termination_condition = _get_reinit_termination_condition(cache, abstol, reltol,
         termination_condition)
+
     cache.abstol = abstol
     cache.reltol = reltol
     cache.termination_condition = termination_condition

diff --git a/src/trustRegion.jl b/src/trustRegion.jl
@@ -736,9 +736,8 @@ function SciMLBase.reinit!(cache::TrustRegionCache{iip}, u0 = cache.u; p = cache
         cache.u = u0
         cache.fu = cache.f(cache.u, p)
     end
-    termination_condition = _get_reinit_termination_condition(cache,
-        abstol,
-        reltol,
+
+    termination_condition = _get_reinit_termination_condition(cache, abstol, reltol,
         termination_condition)
 
     cache.abstol = abstol

diff --git a/src/utils.jl b/src/utils.jl
@@ -215,46 +215,44 @@ function _get_tolerance(η, tc_η, ::Type{T}) where {T}
     return T(ifelse(η !== nothing, η, ifelse(tc_η !== nothing, tc_η, fallback_η)))
 end
 
-function _init_termination_elements(abstol,
-    reltol,
-    termination_condition,
-    ::Type{T}; mode = NLSolveTerminationMode.NLSolveDefault) where {T}
+function _init_termination_elements(abstol, reltol, termination_condition,
+    ::Type{T}; mode = NLSolveTerminationMode.AbsNorm) where {T}
     if termination_condition !== nothing
-        abstol !== nothing ?
-        (abstol != termination_condition.abstol ?
-         error("Incompatible absolute tolerances found. The tolerances supplied as the keyword argument and the one supplied in the termination condition should be same.") :
-         nothing) : nothing
-        reltol !== nothing ?
-        (reltol != termination_condition.abstol ?
-         error("Incompatible relative tolerances found. The tolerances supplied as the keyword argument and the one supplied in the termination condition should be same.") :
-         nothing) : nothing
+        if abstol !== nothing && abstol != termination_condition.abstol
+            error("Incompatible absolute tolerances found. The tolerances supplied as the \
+                keyword argument and the one supplied in the termination condition should \
+                be same.")
+        end
+        if reltol !== nothing && reltol != termination_condition.reltol
+            error("Incompatible relative tolerances found. The tolerances supplied as the \
+                keyword argument and the one supplied in the termination condition should \
+                be same.")
+        end
         abstol = _get_tolerance(abstol, termination_condition.abstol, T)
         reltol = _get_tolerance(reltol, termination_condition.reltol, T)
         return abstol, reltol, termination_condition
     else
         abstol = _get_tolerance(abstol, nothing, T)
         reltol = _get_tolerance(reltol, nothing, T)
-        termination_condition = NLSolveTerminationCondition(mode;
-            abstol,
-            reltol)
+        termination_condition = NLSolveTerminationCondition(mode; abstol, reltol)
         return abstol, reltol, termination_condition
     end
 end
 
 function _get_reinit_termination_condition(cache, abstol, reltol, termination_condition)
     if termination_condition != cache.termination_condition
-        if abstol != cache.abstol
-            if abstol != termination_condition.abstol
-                error("Incompatible absolute tolerances found. The tolerances supplied as the keyword argument and the one supplied in the termination condition should be same.")
-            end
+        if abstol != cache.abstol && abstol != termination_condition.abstol
+            error("Incompatible absolute tolerances found. The tolerances supplied as the \
+                   keyword argument and the one supplied in the termination condition \
+                   should be same.")
         end
 
-        if reltol != cache.reltol
-            if reltol != termination_condition.reltol
-                error("Incompatible absolute tolerances found. The tolerances supplied as the keyword argument and the one supplied in the termination condition should be same.")
-            end
+        if reltol != cache.reltol && reltol != termination_condition.reltol
+            error("Incompatible absolute tolerances found. The tolerances supplied as the \
+                   keyword argument and the one supplied in the termination condition \
+                   should be same.")
         end
-        termination_condition
+        return termination_condition
     else
         # Build the termination_condition with new abstol and reltol
         return NLSolveTerminationCondition{

diff --git a/test/23_test_problems.jl b/test/23_test_problems.jl
@@ -11,8 +11,9 @@ function test_on_library(problems, dicts, alg_ops, broken_tests, ϵ = 1e-4)
         @testset "$idx: $(dict["title"])" begin
             for alg in alg_ops
                 try
-                    sol = solve(nlprob, alg, abstol = 1e-18, reltol = 1e-18)
+                    sol = solve(nlprob, alg)
                     problem(res, sol.u, nothing)
+
                     broken = idx in broken_tests[alg] ? true : false
                     @test norm(res)≤ϵ broken=broken
                 catch