Debug instability

Author: termi-official

lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl

@@ -12,7 +12,7 @@ Algorithm 2 from the classical work by Eisenstat and Walker (1996) as described
     safeguard_threshold
 end
 
-function EisenstatWalkerForcing2(; η₀ = 0.5, ηₘₐₓ = 0.99, γ = 0.9, α = 2, safeguard = true, safeguard_threshold = 0.1)
+function EisenstatWalkerForcing2(; η₀ = 0.5, ηₘₐₓ = 0.9, γ = 0.9, α = 2, safeguard = true, safeguard_threshold = 0.1)
     EisenstatWalkerForcing2(η₀, ηₘₐₓ, γ, α, safeguard, safeguard_threshold)
 end
 
@@ -32,6 +32,7 @@ function pre_step_forcing!(cache::EisenstatWalkerForcing2Cache, descend_cache::N
     # On the first iteration we initialize η with the default initial value and stop.
     if iter == 0
         cache.η = cache.p.η₀
+        @SciMLMessage("Eisenstat-Walker initial iteration to η=$(cache.η).", cache.verbosity, :linear_verbosity)
         LinearSolve.update_tolerances!(descend_cache.lincache; reltol=cache.η)
         return nothing
     end
@@ -56,7 +57,7 @@ function pre_step_forcing!(cache::EisenstatWalkerForcing2Cache, descend_cache::N
     # Far away from the root we also need to respect η ∈ [0,1)
     cache.η = clamp(cache.η, 0.0, cache.p.ηₘₐₓ)
 
-    @SciMLMessage("Eisenstat-Walker update to η=$(cache.η).", cache.verbosity, :linear_verbosity)
+    @SciMLMessage("Eisenstat-Walker iter $iter update to η=$(cache.η).", cache.verbosity, :linear_verbosity)
 
     # Communicate new relative tolerance to linear solve
     LinearSolve.update_tolerances!(descend_cache.lincache; reltol=cache.η)
@@ -69,8 +70,9 @@ end
 function post_step_forcing!(cache::EisenstatWalkerForcing2Cache, J, u, fu, δu, iter)
     # Cache previous residual norm
     cache.rnorm_prev = cache.rnorm
-    # And update the current one
     cache.rnorm = cache.internalnorm(fu)
+
+    # @SciMLMessage("Eisenstat-Walker sanity check: $(cache.internalnorm(fu + J*δu)) ≤ $(cache.η * cache.internalnorm(fu)).", cache.verbosity, :linear_verbosity)
 end
 
 

lib/NonlinearSolveFirstOrder/src/levenberg_marquardt.jl

@@ -57,7 +57,8 @@ function LevenbergMarquardt(;
         vjp_autodiff,
         jvp_autodiff,
         name = :LevenbergMarquardt,
-        concrete_jac = Val(true)
+        concrete_jac = Val(true),
+
     )
 end
 

lib/NonlinearSolveFirstOrder/src/raphson.jl

@@ -1,7 +1,8 @@
 """
     NewtonRaphson(;
         concrete_jac = nothing, linsolve = nothing, linesearch = missing,
-        autodiff = nothing, vjp_autodiff = nothing, jvp_autodiff = nothing
+        autodiff = nothing, vjp_autodiff = nothing, jvp_autodiff = nothing,
+        forcing = nothing,
     )
 
 An advanced NewtonRaphson implementation with support for efficient handling of sparse
@@ -10,13 +11,15 @@ for large-scale and numerically-difficult nonlinear systems.
 """
 function NewtonRaphson(;
         concrete_jac = nothing, linsolve = nothing, linesearch = missing,
-        autodiff = nothing, vjp_autodiff = nothing, jvp_autodiff = nothing
+        autodiff = nothing, vjp_autodiff = nothing, jvp_autodiff = nothing,
+        forcing = nothing,
 )
     return GeneralizedFirstOrderAlgorithm(;
         linesearch,
         descent = NewtonDescent(; linsolve),
         autodiff, vjp_autodiff, jvp_autodiff,
         concrete_jac,
+        forcing,
         name = :NewtonRaphson
     )
 end

lib/NonlinearSolveFirstOrder/src/solve.jl

@@ -314,6 +314,10 @@ function InternalAPI.step!(
     δu, descent_intermediates = descent_result.δu, descent_result.extras
 
     if descent_result.success
+        if has_forcing
+            post_step_forcing!(cache.forcing_cache, J, cache.u, cache.fu, δu, cache.nsteps)
+        end
+
         cache.make_new_jacobian = true
         if cache.globalization isa Val{:LineSearch}
             @static_timeit cache.timer "linesearch" begin
@@ -360,10 +364,6 @@ function InternalAPI.step!(
                    are (:LineSearch, :TrustRegion, :None)")
         end
         NonlinearSolveBase.check_and_update!(cache, cache.fu, cache.u, cache.u_cache)
-
-        if has_forcing
-            post_step_forcing!(cache.forcing_cache, J, cache.u, cache.fu, δu, cache.nsteps)
-        end
     else
         α = false
         cache.make_new_jacobian = false

lib/NonlinearSolveFirstOrder/src/trust_region.jl

@@ -6,7 +6,7 @@
         shrink_threshold::Real = 1 // 4, expand_threshold::Real = 3 // 4,
         shrink_factor::Real = 1 // 4, expand_factor::Real = 2 // 1,
         max_shrink_times::Int = 32,
-        vjp_autodiff = nothing, autodiff = nothing, jvp_autodiff = nothing
+        vjp_autodiff = nothing, autodiff = nothing, jvp_autodiff = nothing,
     )
 
 An advanced TrustRegion implementation with support for efficient handling of sparse
@@ -29,7 +29,7 @@ function TrustRegion(;
         shrink_threshold::Real = 1 // 4, expand_threshold::Real = 3 // 4,
         shrink_factor::Real = 1 // 4, expand_factor::Real = 2 // 1,
         max_shrink_times::Int = 32,
-        autodiff = nothing, vjp_autodiff = nothing, jvp_autodiff = nothing
+        autodiff = nothing, vjp_autodiff = nothing, jvp_autodiff = nothing,
 )
     descent = Dogleg(; linsolve)
     trustregion = GenericTrustRegionScheme(;