NSA: gate iter-1 convergence on inner solver retcode (#3817)

Author: singhharsh1708

lib/OrdinaryDiffEqNonlinearSolve/src/OrdinaryDiffEqNonlinearSolve.jl

@@ -14,7 +14,8 @@ import DiffEqBase: OrdinaryDiffEqTag, calculate_residuals, calculate_residuals!,
     BrownFullBasicInit, ShampineCollocationInit
 import ConstructionBase
 import PreallocationTools: DiffCache, get_tmp
-using SimpleNonlinearSolve: SimpleTrustRegion, SimpleGaussNewton
+using SimpleNonlinearSolve: SimpleTrustRegion, SimpleGaussNewton,
+    AbstractSimpleNonlinearSolveAlgorithm
 using NonlinearSolve: FastShortcutNonlinearPolyalg, FastShortcutNLLSPolyalg, NewtonRaphson,
     HomotopySweep, step!
 # The operator Jacobian path is implemented in NonlinearSolveBase and needs its own floor.

lib/OrdinaryDiffEqNonlinearSolve/src/nlsolve.jl

@@ -1,5 +1,21 @@
 @inline eps_around_one(θ::T) where {T} = 100sqrt(eps(one(θ)))
 
+# The `iter == 1 && ndz < 1e-5` shortcut below assumes a plain (modified-)Newton
+# contraction: a tiny first-iterate displacement means the initial guess was
+# already the solution. For `NonlinearSolveAlg` backed by a globalized inner
+# solver (e.g. TrustRegion), a tiny displacement can instead mean the inner
+# solver took a cautious/rejected trial step and has not decided anything yet —
+# accepting that as convergence applies zero correction and silently produces
+# the wrong stage value (#3817). Only the inner cache's own retcode
+# distinguishes "genuinely done already" from "hasn't started converging".
+# Iterations beyond the first use `η`/`θ` computed from actual step history and
+# are not vulnerable to this false positive, so they are left untouched.
+_nsa_inner_converged(nlsolver) = true
+function _nsa_inner_converged(nlsolver::NLSolver{<:NonlinearSolveAlg})
+    nlsolver.alg.alg isa AbstractSimpleNonlinearSolveAlgorithm && return true
+    return SciMLBase.successful_retcode(nlsolver.cache.cache.retcode)
+end
+
 """
     nlsolve!(nlsolver::AbstractNLSolver, integrator)
 
@@ -136,7 +152,7 @@ function nlsolve!(
                 )
             )
         )
-        if (iter == 1 && ndz < 1.0e-5) ||
+        if (iter == 1 && ndz < 1.0e-5 && _nsa_inner_converged(nlsolver)) ||
                 (check_η_convergence && η >= zero(η) && η * ndz < κ)
             @SciMLMessage(
                 lazy"Newton iteration converged in $(iter) iterations: η = $(η), ndz = $(ndz)",

lib/OrdinaryDiffEqNonlinearSolve/test/nsa_jacobian_reuse_tests.jl

@@ -1,7 +1,7 @@
 using OrdinaryDiffEqBDF, OrdinaryDiffEqSDIRK
 using OrdinaryDiffEqNonlinearSolve
 using OrdinaryDiffEqNonlinearSolve: NonlinearSolveAlg
-using NonlinearSolve: NewtonRaphson
+using NonlinearSolve: NewtonRaphson, TrustRegion
 using ADTypes, LinearAlgebra, SciMLBase
 using Test
 
@@ -48,3 +48,30 @@ refsol = solve(prob, FBDF(); reltol = 1.0e-12, abstol = 1.0e-14)
         @test JAC_CALLS[] >= 1
     end
 end
+
+@testset "globalized inner solver does not converge on a rejected step" begin
+    # A TrustRegion step that is rejected/truncated leaves the iterate nearly
+    # unmoved, which the outer displacement test alone reads as convergence — the
+    # stage is then accepted with no correction applied (#3817). Driven at a dt
+    # that is far too coarse for the Robertson transient, so the inner solves
+    # genuinely cannot converge and the failure must be reported rather than
+    # silently absorbed.
+    nsa_tr = NonlinearSolveAlg(TrustRegion(; autodiff = AutoForwardDiff()))
+    nsa_nr = NonlinearSolveAlg(NewtonRaphson(; autodiff = AutoForwardDiff()))
+    hard = ODEProblem(f, [1.0, 0.0, 0.0], (0.0, 1.0e3), [0.04, 3.0e7, 1.0e4])
+
+    sol_tr = solve(hard, FBDF(nlsolve = nsa_tr); dt = 1.0, adaptive = false)
+    sol_nr = solve(hard, FBDF(nlsolve = nsa_nr); dt = 1.0, adaptive = false)
+
+    # The state must not be reported as a successful solve while frozen at u0.
+    @test !(SciMLBase.successful_retcode(sol_tr) && sol_tr.u[end] == hard.u0)
+    # TrustRegion must reach the same verdict as the non-globalized inner solver.
+    @test SciMLBase.successful_retcode(sol_tr) == SciMLBase.successful_retcode(sol_nr)
+end
+
+@testset "TrustRegion matches NewtonRaphson when the solves do converge" begin
+    nsa_tr = NonlinearSolveAlg(TrustRegion(; autodiff = AutoForwardDiff()))
+    sol = solve(prob, FBDF(nlsolve = nsa_tr); reltol = 1.0e-8, abstol = 1.0e-10)
+    @test SciMLBase.successful_retcode(sol)
+    @test norm(sol.u[end] .- refsol.u[end]) / norm(refsol.u[end]) < 1.0e-4
+end