Fix resize, algorithm-switch, and OOP W-caching for Jacobian reuse

Three robustness fixes for Rosenbrock-W Jacobian reuse:

1. Resize detection: Track last_u_length in JacReuseState. After a
   callback resize!, u_modified is cleared before perform_step! runs,
   so we check length(u) directly to detect dimension changes.

2. CompositeAlgorithm switch detection: Track last_step_iter to detect
   gaps in integrator.iter, which indicate another algorithm (e.g.,
   Vern8) ran between Rosenbrock steps. Forces J recomputation on
   switch-back to avoid using a stale Jacobian.

3. OOP W caching: Honor the new_W=false flag by caching and reusing
   the factorized W (LU), matching IIP behavior. Previously the OOP
   path always rebuilt W from stale J, which masked inaccuracy and
   prevented the self-correcting rejection feedback loop. This fixes
   AutoVern8(Rosenbrock23()) OOP hitting MaxIters on stiff problems.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

Author: ChrisRackauckas

lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl

@@ -63,11 +63,26 @@ function _rosenbrock_jac_reuse_decision(integrator, cache, dtgamma)
         return (true, true)
     end
 
+    # Detect algorithm switch in CompositeAlgorithm: if integrator.iter jumped
+    # by more than 1 since our last Rosenbrock step, another algorithm ran in
+    # between and the cached Jacobian is evaluated at a stale u.
+    if jac_reuse.last_step_iter != 0 && integrator.iter > jac_reuse.last_step_iter + 1
+        return (true, true)
+    end
+
     # Callback modification: recompute
     if integrator.u_modified
         return (true, true)
     end
 
+    # Resize detection: if u changed length since last J computation,
+    # the cached LU factorization has wrong dimensions.
+    # (u_modified is already cleared by reeval_internals_due_to_modification!
+    #  before perform_step! runs, so we need this explicit check.)
+    if length(integrator.u) != jac_reuse.last_u_length && jac_reuse.last_u_length != 0
+        return (true, true)
+    end
+
     # Previous step was rejected (EEst > 1): the old W wasn't good enough.
     # Recompute everything since we're retrying with a different dt anyway.
     if integrator.EEst > 1
@@ -792,8 +807,12 @@ function calc_rosenbrock_differentiation!(integrator, cache, dtd1, dtgamma, repe
         # _rosenbrock_jac_reuse_decision). This tracks the dtgamma at the
         # last J computation for the gamma ratio heuristic.
         jac_reuse = get_jac_reuse(cache)
-        if jac_reuse !== nothing && new_jac
-            jac_reuse.pending_dtgamma = dtgamma
+        if jac_reuse !== nothing
+            jac_reuse.last_step_iter = integrator.iter
+            if new_jac
+                jac_reuse.pending_dtgamma = dtgamma
+                jac_reuse.last_u_length = length(integrator.u)
+            end
         end
     end
     # If the Jacobian is not updated, we won't have to update ∂/∂t either.
@@ -829,6 +848,9 @@ function calc_rosenbrock_differentiation(integrator, cache, dtgamma, repeat_step
 
     new_jac, new_W = newJW
 
+    # Track iteration for algorithm-switch detection in CompositeAlgorithm
+    jac_reuse.last_step_iter = integrator.iter
+
     # For complex W types (operators), delegate to standard calc_W
     if cache.W isa StaticWOperator || cache.W isa WOperator ||
             cache.W isa AbstractSciMLOperator
@@ -841,10 +863,14 @@ function calc_rosenbrock_differentiation(integrator, cache, dtgamma, repeat_step
     mass_matrix = integrator.f.mass_matrix
     update_coefficients!(mass_matrix, integrator.uprev, integrator.p, integrator.t)
 
-    # Safety: if cached_J is nothing (e.g. first use after algorithm switch),
+    # Safety: if cached_J or cached_W is nothing (e.g. first use after algorithm switch),
     # force a fresh computation regardless of the decision.
     if !new_jac && jac_reuse.cached_J === nothing
         new_jac = true
+        new_W = true
+    end
+    if !new_W && jac_reuse.cached_W === nothing
+        new_W = true
     end
 
     if new_jac
@@ -864,14 +890,23 @@ function calc_rosenbrock_differentiation(integrator, cache, dtgamma, repeat_step
     # committed as last_dtgamma on the next accepted step.
     if new_jac
         jac_reuse.pending_dtgamma = dtgamma
+        jac_reuse.last_u_length = length(integrator.u)
     end
 
-    # Build W from (possibly cached) J
-    W = J - mass_matrix * inv(dtgamma)
-    if !isa(W, Number)
-        W = DiffEqBase.default_factorize(W)
+    # 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.
+    if new_W
+        W = J - mass_matrix * inv(dtgamma)
+        if !isa(W, Number)
+            W = DiffEqBase.default_factorize(W)
+        end
+        integrator.stats.nw += 1
+        jac_reuse.cached_W = W
+    else
+        W = jac_reuse.cached_W
     end
-    integrator.stats.nw += 1
 
     return dT, W
 end

lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl

@@ -16,6 +16,8 @@ Fields:
 - `cached_J`: Cached Jacobian for OOP reuse (type-erased for flexibility)
 - `cached_dT`: Cached time derivative for OOP reuse
 - `cached_W`: Cached factorized W for OOP reuse (LU factorization is expensive)
+- `last_u_length`: Length of u when J/W were last computed (detects resize!)
+- `last_step_iter`: `integrator.iter` at last Rosenbrock step (detects algorithm switches)
 """
 mutable struct JacReuseState{T}
     last_dtgamma::T
@@ -25,10 +27,12 @@ mutable struct JacReuseState{T}
     cached_J::Any
     cached_dT::Any
     cached_W::Any
+    last_u_length::Int
+    last_step_iter::Int
 end
 
 function JacReuseState(dtgamma::T) where {T}
-    return JacReuseState{T}(dtgamma, dtgamma, 0, 50, nothing, nothing, nothing)
+    return JacReuseState{T}(dtgamma, dtgamma, 0, 50, nothing, nothing, nothing, 0, 0)
 end
 
 # Fake values since non-FSAL

lib/OrdinaryDiffEqRosenbrock/test/jacobian_reuse_test.jl

@@ -124,6 +124,8 @@ strict_rosenbrock = [
         @test jr.cached_J === nothing
         @test jr.cached_dT === nothing
         @test jr.cached_W === nothing
+        @test jr.last_u_length == 0
+        @test jr.last_step_iter == 0
     end
 
     # ========================================================================