Fix downstream tolerance failures: reuse J only, always rebuild W

- Disable Jacobian reuse for mass-matrix (DAE) problems to avoid order
  reduction from stale algebraic constraint derivatives (addresses
  @gstein3m's review comment)
- Change reuse strategy from (false,false) to (false,true): reuse the
  cached Jacobian but always rebuild W = J − M/(dt·γ) with the current
  dtgamma. The Jacobian evaluation (finite-diff/AD) is the expensive
  part; W construction + LU factorization is comparatively cheap and
  keeps step control accurate.
- Remove stale OrdinaryDiffEqNonlinearSolve from runtime [deps] in
  OrdinaryDiffEqRosenbrock (fixes Aqua stale dependencies QA test)

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Author: ChrisRackauckas

lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl

@@ -12,16 +12,20 @@ get_jac_reuse(cache) = hasproperty(cache, :jac_reuse) ? cache.jac_reuse : nothin
     _rosenbrock_jac_reuse_decision(integrator, cache, dtgamma) -> Union{Nothing, NTuple{2,Bool}}
 
 Decide whether to recompute the Jacobian and/or W matrix for Rosenbrock methods.
-For W-methods (where `isWmethod(alg) == true`), implements CVODE-inspired reuse:
+For W-methods (where `isWmethod(alg) == true`) on non-DAE problems, implements
+CVODE-inspired Jacobian reuse:
 - Always recompute on first iteration
-- Recompute after step rejection (EEst > 1), since the old W wasn't good enough
+- Recompute after step rejection (EEst > 1), since the old J wasn't good enough
 - Recompute when gamma ratio changes too much: |dtgamma/last_dtgamma - 1| > 0.3
 - Recompute every `max_jac_age` accepted steps (default 50)
 - Recompute when u_modified (callback modification)
-- Otherwise reuse both J and W (including LU factorization). The LU is the
-  expensive part, so we try the old W and only recompute if the step fails.
+- Otherwise reuse J but rebuild W from the cached J and current dtgamma.
+  The Jacobian evaluation (finite-diff / AD) is the expensive part; W
+  construction and LU factorization are comparatively cheap.
 
-For strict Rosenbrock methods, returns `nothing` to delegate to `do_newJW`.
+Returns `nothing` (delegate to `do_newJW`) for strict Rosenbrock methods,
+linear problems, and mass-matrix (DAE) problems where stale Jacobians
+cause order reduction.
 """
 function _rosenbrock_jac_reuse_decision(integrator, cache, dtgamma)
     alg = OrdinaryDiffEqCore.unwrap_alg(integrator, true)
@@ -43,6 +47,14 @@ function _rosenbrock_jac_reuse_decision(integrator, cache, dtgamma)
         return nothing
     end
 
+    # Mass matrix (DAE) problems: delegate to do_newJW.
+    # Stale Jacobians cause order reduction for DAEs because the algebraic
+    # constraint derivatives must remain accurate.  See Steinebach (2024)
+    # for W-method DAE order conditions.
+    if integrator.f.mass_matrix !== I
+        return nothing
+    end
+
     # First iteration: always compute J and W.
     if integrator.iter <= 1
         return (true, true)
@@ -102,9 +114,11 @@ function _rosenbrock_jac_reuse_decision(integrator, cache, dtgamma)
         return (true, true)
     end
 
-    # Reuse both J and W. The LU factorization is expensive, so try with
-    # the old W and only recompute if the step is rejected (caught above).
-    return (false, false)
+    # Reuse J but rebuild W with the current dtgamma.  Following CVODE's
+    # approach: the Jacobian evaluation (finite-diff / AD) is expensive,
+    # while reconstructing W = J − M/(dt·γ) and its LU is comparatively
+    # cheap and keeps the step controller accurate.
+    return (false, true)
 end
 
 function calc_tderivative!(integrator, cache, dtd1, repeat_step)
@@ -893,10 +907,9 @@ function calc_rosenbrock_differentiation(integrator, cache, dtgamma, repeat_step
         jac_reuse.last_u_length = length(integrator.u)
     end
 
-    # Build and cache W, or reuse cached W (including LU factorization).
-    # Reusing old W (with old dtgamma) mirrors IIP behavior: if the old W
-    # is too inaccurate, the step will be rejected and EEst > 1 triggers
-    # a fresh J+W computation on the retry.
+    # Always rebuild W from the (possibly cached) J and the current dtgamma.
+    # The Jacobian evaluation is the expensive part; W = J − M/(dt·γ) and
+    # its factorization are comparatively cheap and keep step control accurate.
     if new_W
         W = J - mass_matrix * inv(dtgamma)
         if !isa(W, Number)

lib/OrdinaryDiffEqRosenbrock/Project.toml

@@ -16,7 +16,6 @@ MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
 OrdinaryDiffEqDifferentiation = "4302a76b-040a-498a-8c04-15b101fed76b"
-OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8"
 Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
@@ -28,7 +27,6 @@ Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
 [sources]
 OrdinaryDiffEqCore = {path = "../OrdinaryDiffEqCore"}
 OrdinaryDiffEqDifferentiation = {path = "../OrdinaryDiffEqDifferentiation"}
-OrdinaryDiffEqNonlinearSolve = {path = "../OrdinaryDiffEqNonlinearSolve"}
 
 [compat]
 ADTypes = "1.16"