NSA: don't count rejected inner trial steps as convergence (#3817) (#3893)

A globalized inner solver (TrustRegion and friends) can reject its trial
step: `step!` returns with the iterate exactly unmoved and the cache not
terminated, having only shrunk its trust region for the next attempt.
The outer convergence test read that zero displacement as a perfect
solve - `ndz < 1e-5` on the first iteration, `eta*ndz = 0 < kappa` on
later ones - and accepted the stage with no correction applied. #4020
surfaces inner caches that *terminated* unsuccessfully; a rejected step
terminates nothing, so it slipped through both branches.

`nlsolve!` now skips the convergence and divergence bookkeeping for an
iteration whose `NonlinearSolveAlg` step left the iterate unmoved
without terminating: no new iterate exists to judge, and a zero `ndz`
would also poison the next iteration's theta. The inner solver retries
with its shrunken radius; if it never moves, the loop runs out and the
step is rejected as unconverged. A cache that terminated at zero
displacement converged exactly (failures already return `Inf` from
`compute_step!`), so that case still counts as convergence.

ROBER with FBDF and a TrustRegion inner solver:

- fixed dt = 1.0 (transient unresolvable): was ReturnCode.Success with
  the state frozen at u0; now ConvergenceFailure, matching NewtonRaphson.
- adaptive, default tolerances: relerr 1.7e-2 -> 1.8e-5 (NLNewton 2.1e-5,
  NSA NewtonRaphson 4.0e-6). Same for RobustMultiNewton.
- adaptive, reltol=1e-8: relerr 2.4e-8 -> 1.6e-9 (NewtonRaphson 1.3e-9).

NewtonRaphson inner solves are unchanged step-for-step (their iterates
always move or terminate), and NLNewton never enters the predicate.

Author: singhharsh1708

lib/OrdinaryDiffEqNonlinearSolve/src/nlsolve.jl

@@ -1,5 +1,19 @@
 @inline eps_around_one(θ::T) where {T} = 100sqrt(eps(one(θ)))
 
+# A globalized inner solver (TrustRegion and friends) can reject its trial step: `step!`
+# returns with the iterate exactly unmoved and the cache not terminated, having only shrunk
+# its trust region for the next attempt. The zero displacement such a `step!` leaves behind
+# is not convergence evidence — the convergence tests below would read it as a perfect solve
+# (`ndz < 1e-5` on the first iteration, `η·ndz = 0 < κ` on later ones) and accept the stage
+# with no correction applied (#3817). A cache that instead *terminated* at zero displacement
+# either converged exactly (which genuinely is convergence) or failed, and failures were
+# already turned into `Inf` by `compute_step!`, so termination is the discriminator between
+# "already at the root" and "has not decided anything yet".
+_rejected_trial_step(nlsolver, ndz) = false
+function _rejected_trial_step(nlsolver::NLSolver{<:NonlinearSolveAlg}, ndz)
+    return iszero(ndz) && NonlinearSolveBase.not_terminated(nlsolver.cache.cache)
+end
+
 """
     nlsolve!(nlsolver::AbstractNLSolver, integrator)
 
@@ -70,6 +84,19 @@ function nlsolve!(
             break
         end
 
+        if _rejected_trial_step(nlsolver, ndz)
+            @SciMLMessage(
+                lazy"Inner nonlinear solver rejected its trial step (iter = $(iter)); the unmoved iterate is not treated as convergence",
+                integrator.opts.verbose, :newton_convergence
+            )
+            # No new iterate exists to judge: skip the convergence and divergence
+            # bookkeeping (a zero `ndz` would also poison the next iteration's `θ`)
+            # and let the inner solver retry with its shrunken trust region. If it
+            # never moves, the loop runs out and the step fails as unconverged.
+            ndz = ndzprev
+            continue
+        end
+
         has_prev_θ = hasfield(NL, :prev_θ)
         prev_θ = has_prev_θ ? nlsolver.prev_θ : one(ndz)
 

lib/OrdinaryDiffEqNonlinearSolve/test/nsa_jacobian_reuse_tests.jl

@@ -1,7 +1,7 @@
 using OrdinaryDiffEqBDF, OrdinaryDiffEqSDIRK, OrdinaryDiffEqRosenbrock
 using OrdinaryDiffEqNonlinearSolve
 using OrdinaryDiffEqNonlinearSolve: NonlinearSolveAlg
-using NonlinearSolve: NewtonRaphson
+using NonlinearSolve: NewtonRaphson, TrustRegion
 using ADTypes, LinearAlgebra, SciMLBase
 using Test
 
@@ -65,3 +65,38 @@ end
     @test err(sol) < 3 * err(ref)
     @test sol.stats.naccept < 1.25 * ref.stats.naccept
 end
+
+@testset "globalized inner solver does not converge on a rejected step" begin
+    # A TrustRegion trial step that is rejected leaves the iterate exactly unmoved
+    # without terminating the inner cache; the outer displacement test alone reads
+    # that as convergence and accepts the stage with no correction applied (#3817).
+    # Driven at a dt 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
+
+    # At default tolerances the rejected-step false convergence used to cost three
+    # orders of magnitude (relerr 1.7e-2, where NLNewton gives 2.1e-5); with rejected
+    # trial steps excluded from the convergence test it lands back at tolerance level
+    # (measured 1.8e-5).
+    sol_def = solve(prob, FBDF(nlsolve = nsa_tr))
+    @test SciMLBase.successful_retcode(sol_def)
+    @test norm(sol_def.u[end] .- refsol.u[end]) / norm(refsol.u[end]) < 1.0e-3
+end