From 6dacbdd288617080fc85fcd8763847dfa72497d5 Mon Sep 17 00:00:00 2001 From: "Chris Rackauckas (Claude)" Date: Tue, 9 Jun 2026 22:41:53 -0400 Subject: [PATCH 1/2] Fix Float32/static-array mass-matrix OOP Rosenbrock JacReuseState W typing Fixes #3721. Two failure modes in the OOP Jacobian-reuse path for static-array mass-matrix problems: 1. JacReuseState's cached_W slot was seeded with build_J_W's W, whose eltype follows the state (e.g. Float32), but calc_W promotes the W eltype with typeof(dtgamma) (e.g. Float64 when dt promotes the time type). The fresh-Jacobian branch's `jac_reuse.cached_W = W` then throws, since StaticWOperator has no eltype-converting convert. 2. The cached-J reuse branch rebuilt W via default_factorize, producing a StaticArrays.LU instead of the StaticWOperator the slot holds, so any OOP static mass-matrix solve with a rejected step threw even in homogeneous Float64. Fix: add build_jac_reuse_W_seed, which builds the cached_W seed with the same concrete type calc_W produces at solve time, seed the OOP alg_caches with it and with zero(dt * gamma) for the dtgamma fields, and make the reuse branch mirror calc_W's StaticWOperator branch. OrdinaryDiffEqDifferentiation 3.2.0 -> 3.2.1, OrdinaryDiffEqRosenbrock 2.3.0 -> 2.3.1 (requires the new function, so the compat lower bound is raised accordingly). Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Fable 5 --- .../Project.toml | 2 +- .../src/derivative_utils.jl | 47 +++++++++++-- lib/OrdinaryDiffEqRosenbrock/Project.toml | 8 ++- .../src/OrdinaryDiffEqRosenbrock.jl | 1 + .../src/rosenbrock_caches.jl | 29 +++++--- .../test/ode_rosenbrock_tests.jl | 70 +++++++++++++++++++ 6 files changed, 138 insertions(+), 19 deletions(-) diff --git a/lib/OrdinaryDiffEqDifferentiation/Project.toml b/lib/OrdinaryDiffEqDifferentiation/Project.toml index a0533874e3a..2d27b2667c9 100644 --- a/lib/OrdinaryDiffEqDifferentiation/Project.toml +++ b/lib/OrdinaryDiffEqDifferentiation/Project.toml @@ -1,5 +1,5 @@ name = "OrdinaryDiffEqDifferentiation" -version = "3.2.0" +version = "3.2.1" uuid = "4302a76b-040a-498a-8c04-15b101fed76b" authors = ["Chris Rackauckas ", "Yingbo Ma "] diff --git a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl index 5b2263ea38c..c562fd29447 100644 --- a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl +++ b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl @@ -937,12 +937,17 @@ 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) - # Rebuild W from cached J and current dtgamma. + # Rebuild W from cached J and current dtgamma, mirroring calc_W's branches + # so W (and hence the cached_W slot) keeps the same concrete type. # The Jacobian evaluation is the expensive part; W = J − M/(dt·γ) and # its factorization are comparatively cheap and keep step control accurate. - W = J - mass_matrix * inv(dtgamma) - if !isa(W, Number) - W = DiffEqBase.default_factorize(W) + if cache.W isa StaticWOperator + W = StaticWOperator(J - mass_matrix * inv(dtgamma)) + else + W = J - mass_matrix * inv(dtgamma) + if !isa(W, Number) + W = DiffEqBase.default_factorize(W) + end end integrator.stats.nw += 1 jac_reuse.cached_W = W @@ -1184,6 +1189,40 @@ function build_J_W( return J, W end +""" + build_jac_reuse_W_seed(W, J, u, mass_matrix, dtgamma) + +Build the seed value for `JacReuseState`'s `cached_W` slot in OOP caches. + +The seed's *type* must match what `calc_W` produces at solve time, which +promotes the W eltype with `typeof(dtgamma)` and `eltype(mass_matrix)` — e.g. +a `Float32` state solved over a `Float64` time span yields a `Float64`-eltype +W. Seeding with `build_J_W`'s state-eltype `W` would make the +`jac_reuse.cached_W = W` assignment in `calc_rosenbrock_differentiation` throw +for wrapper types without an eltype-converting `convert` (`StaticWOperator`). + +`dtgamma` must be a prototype carrying the solve-time `dt * gamma` type; only +its type is used, so the seed's values are garbage. That is safe because the +"first use" guard in `_rosenbrock_jac_reuse_decision` guarantees the seed is +overwritten before it is ever read. +""" +function build_jac_reuse_W_seed(W, J, u, mass_matrix, dtgamma) + invdtgamma′ = inv(oneunit(dtgamma)) + if W isa StaticWOperator + # callinv = false skips inverting the garbage-valued seed while still + # producing the same concrete type as calc_W's + # `StaticWOperator(J - mass_matrix * inv(dtgamma))`. + return StaticWOperator(J - mass_matrix * invdtgamma′, false) + elseif W isa WOperator + return WOperator{false}(mass_matrix, dtgamma, J, u, W.jacvec) + elseif W isa AbstractSciMLOperator + return W + else + _W = J - mass_matrix * invdtgamma′ + return _W isa Number ? _W : ArrayInterface.lu_instance(_W) + end +end + build_uf(alg, nf, t, p, ::Val{true}) = UJacobianWrapper(nf, t, p) build_uf(alg, nf, t, p, ::Val{false}) = UDerivativeWrapper(nf, t, p) diff --git a/lib/OrdinaryDiffEqRosenbrock/Project.toml b/lib/OrdinaryDiffEqRosenbrock/Project.toml index d36305fe33c..2468836a8f1 100644 --- a/lib/OrdinaryDiffEqRosenbrock/Project.toml +++ b/lib/OrdinaryDiffEqRosenbrock/Project.toml @@ -1,6 +1,6 @@ name = "OrdinaryDiffEqRosenbrock" uuid = "43230ef6-c299-4910-a778-202eb28ce4ce" -version = "2.3.0" +version = "2.3.1" authors = ["ParamThakkar123 "] [deps] @@ -46,7 +46,7 @@ MacroTools = "0.5.6" MuladdMacro = "0.2.1" ODEProblemLibrary = "1" OrdinaryDiffEqCore = "4" -OrdinaryDiffEqDifferentiation = "3" +OrdinaryDiffEqDifferentiation = "3.2.1" OrdinaryDiffEqNonlinearSolve = "2" OrdinaryDiffEqRosenbrockTableaus = "2.1" Pkg = "1" @@ -56,6 +56,7 @@ Random = "<0.0.1, 1" RecursiveArrayTools = "4" Reexport = "1.2.2" SafeTestsets = "0.1" +StaticArrays = "1.9.18" SciMLBase = "3.1.0" Test = "<0.0.1, 1" julia = "1.10" @@ -68,7 +69,8 @@ OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["DiffEqDevTools", "Random", "OrdinaryDiffEqNonlinearSolve", "SafeTestsets", "Test", "ODEProblemLibrary", "Enzyme", "Pkg"] +test = ["DiffEqDevTools", "Random", "OrdinaryDiffEqNonlinearSolve", "SafeTestsets", "StaticArrays", "Test", "ODEProblemLibrary", "Enzyme", "Pkg"] diff --git a/lib/OrdinaryDiffEqRosenbrock/src/OrdinaryDiffEqRosenbrock.jl b/lib/OrdinaryDiffEqRosenbrock/src/OrdinaryDiffEqRosenbrock.jl index 50153bacbc4..9b13e8b8cd9 100644 --- a/lib/OrdinaryDiffEqRosenbrock/src/OrdinaryDiffEqRosenbrock.jl +++ b/lib/OrdinaryDiffEqRosenbrock/src/OrdinaryDiffEqRosenbrock.jl @@ -33,6 +33,7 @@ using OrdinaryDiffEqDifferentiation: TimeDerivativeWrapper, TimeGradientWrapper, build_jac_config, issuccess_W, jacobian2W!, resize_jac_config!, resize_grad_config!, calc_W, calc_rosenbrock_differentiation!, build_J_W, + build_jac_reuse_W_seed, UJacobianWrapper, dolinsolve, WOperator, resize_J_W! using OrdinaryDiffEqDifferentiation: calc_rosenbrock_differentiation diff --git a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl index 095a3795700..8f0b9ef43a4 100644 --- a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl +++ b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl @@ -11,7 +11,7 @@ Jacobian recomputations when conditions allow it. `T` is the element type of `dtgamma` at solve time — usually `Float64`, but `ForwardDiff.Dual` (including nested Duals under `ForwardDiff.hessian`) when the solve is being differentiated. Callers must construct this with a zero -value of the right type (typically `zero(dt)`) so the same cache can hold +value of the right type (typically `zero(dt * gamma)`) so the same cache can hold Dual-valued dtgammas without any `ForwardDiff.value` unwrapping — unconditional unwrapping is unsafe under nested AD because it collapses a still-active inner derivative into a primal. @@ -342,13 +342,16 @@ function alg_cache( linprob = nothing #LinearProblem(W,copy(u); u0=copy(u)) linsolve = nothing #init(linprob,alg.linsolve,alias_A=true,alias_b=true) tab = Rosenbrock23Tableau(constvalue(uBottomEltypeNoUnits)) - # Seed JacReuseState with `zero(dt)` rather than a `constvalue`-stripped - # type: under nested ForwardDiff (e.g. `hessian`), dt is a Dual-of-Dual and - # dtgamma inherits that full type; a Float64 field would reject the assign. + # Seed JacReuseState with `zero(dt * tab.d)` rather than a + # `constvalue`-stripped type: this carries the exact type dtgamma has at + # solve time, including Dual-of-Dual under nested ForwardDiff + # (e.g. `hessian`), where a Float64 field would reject the assign. + dtgamma_seed = zero(dt * tab.d) return Rosenbrock23ConstantCache( tab.c₃₂, tab.d, tf, uf, J, W, linsolve, alg_autodiff(alg), _make_jac_reuse_state( - zero(dt), alg.max_jac_age, J, zero(rate_prototype), W + dtgamma_seed, alg.max_jac_age, J, zero(rate_prototype), + build_jac_reuse_W_seed(W, J, u, f.mass_matrix, dtgamma_seed) ) ) end @@ -378,12 +381,14 @@ function alg_cache( linprob = nothing #LinearProblem(W,copy(u); u0=copy(u)) linsolve = nothing #init(linprob,alg.linsolve,alias_A=true,alias_b=true) tab = Rosenbrock32Tableau(constvalue(uBottomEltypeNoUnits)) - # See the Rosenbrock23 OOP alg_cache above for why we pass `zero(dt)` here - # rather than a `constvalue`-stripped type. + # See the Rosenbrock23 OOP alg_cache above for why we pass `zero(dt * tab.d)` + # here rather than a `constvalue`-stripped type. + dtgamma_seed = zero(dt * tab.d) return Rosenbrock32ConstantCache( tab.c₃₂, tab.d, tf, uf, J, W, linsolve, alg_autodiff(alg), _make_jac_reuse_state( - zero(dt), alg.max_jac_age, J, zero(rate_prototype), W + dtgamma_seed, alg.max_jac_age, J, zero(rate_prototype), + build_jac_reuse_W_seed(W, J, u, f.mass_matrix, dtgamma_seed) ) ) end @@ -475,14 +480,16 @@ function alg_cache( else interp_order = H_rows end - # Seed JacReuseState with `zero(dt)` so its dtgamma fields carry the full - # (possibly ForwardDiff.Dual) type dtgamma will have at solve time. + # Seed JacReuseState with `zero(dt * tab.gamma)` so its dtgamma fields carry + # the full (possibly ForwardDiff.Dual) type dtgamma will have at solve time. + dtgamma_seed = zero(dt * tab.gamma) return RosenbrockCombinedConstantCache( tf, uf, tab, J, W, linsolve, alg_autodiff(alg), interp_order, _make_jac_reuse_state( - zero(dt), alg.max_jac_age, J, zero(rate_prototype), W + dtgamma_seed, alg.max_jac_age, J, zero(rate_prototype), + build_jac_reuse_W_seed(W, J, u, f.mass_matrix, dtgamma_seed) ) ) end diff --git a/lib/OrdinaryDiffEqRosenbrock/test/ode_rosenbrock_tests.jl b/lib/OrdinaryDiffEqRosenbrock/test/ode_rosenbrock_tests.jl index 2e466aca7eb..937eeb9f1b5 100644 --- a/lib/OrdinaryDiffEqRosenbrock/test/ode_rosenbrock_tests.jl +++ b/lib/OrdinaryDiffEqRosenbrock/test/ode_rosenbrock_tests.jl @@ -429,3 +429,73 @@ end @test err_interp < 100 * max(err_knot, 1.0e-4) end end + +# https://github.com/SciML/OrdinaryDiffEq.jl/issues/3721 +# OOP static-array mass-matrix solves were failing in two ways: +# 1. With a Float32 state and a Float64 time type (e.g. `dt = 0.1` promoting +# tspan), calc_W produces a Float64-eltype StaticWOperator that cannot be +# assigned into JacReuseState's Float32-seeded `cached_W` slot. +# 2. On rejected-step retries the OOP reuse branch rebuilt W via +# `default_factorize`, producing a `StaticArrays.LU` instead of the +# `StaticWOperator` type of the `cached_W` slot. +using StaticArrays +using OrdinaryDiffEqNonlinearSolve: BrownFullBasicInit + +@testset "Static-array mass-matrix OOP JacReuseState (#3721)" begin + function rober_static(u, p, t) + y₁, y₂, y₃ = u + k₁, k₂, k₃ = p + return @SVector [ + -k₁ * y₁ + k₃ * y₂ * y₃, + k₁ * y₁ - k₂ * y₂^2 - k₃ * y₂ * y₃, + y₁ + y₂ + y₃ - 1, + ] + end + + # Float64 reference + M64 = @SMatrix [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 0.0] + prob64 = ODEProblem( + ODEFunction(rober_static, mass_matrix = M64), + @SVector([1.0, 0.0, 0.0]), (0.0, 1.0e5), (0.04, 3.0e7, 1.0e4) + ) + ref = solve( + prob64, Rosenbrock23(), abstol = 1.0e-10, reltol = 1.0e-10, + initializealg = BrownFullBasicInit() + ) + + # Failure mode 1: Float32 state with Float64 dt (promotes the time type). + M32 = @SMatrix [1.0f0 0.0f0 0.0f0; 0.0f0 1.0f0 0.0f0; 0.0f0 0.0f0 0.0f0] + prob32 = ODEProblem( + ODEFunction(rober_static, mass_matrix = M32), + @SVector([1.0f0, 0.0f0, 0.0f0]), (0.0f0, 1.0f5), (0.04f0, 3.0f7, 1.0f4) + ) + sol32 = solve( + prob32, Rosenbrock23(), dt = 0.1, abstol = 1.0f-5, reltol = 1.0f-5, + initializealg = BrownFullBasicInit() + ) + @test SciMLBase.successful_retcode(sol32) + @test norm(sol32.u[end] - ref.u[end]) < 1.0e-4 + + # Failure mode 2: homogeneous Float64 at a tolerance that produces step + # rejections, whose retries take the cached-J reuse branch. nreject > 0 is + # asserted so this keeps covering that branch if step control changes. + sol64 = solve( + prob64, Rosenbrock23(), abstol = 1.0e-5, reltol = 1.0e-5, + initializealg = BrownFullBasicInit() + ) + @test SciMLBase.successful_retcode(sol64) + @test sol64.stats.nreject > 0 + @test norm(sol64.u[end] - ref.u[end]) < 1.0e-4 + + # The Rodas-family OOP constant cache goes through the same cached_W + # seeding path; Rodas23W and ROS34PW2 are W-methods so they additionally + # exercise the Jacobian-reuse decision logic on the static path. + for alg in (Rodas4(), Rodas23W(), ROS34PW2()) + sol_rodas = solve( + prob32, alg, dt = 0.1, abstol = 1.0f-5, reltol = 1.0f-5, + initializealg = BrownFullBasicInit() + ) + @test SciMLBase.successful_retcode(sol_rodas) + @test norm(sol_rodas.u[end] - ref.u[end]) < 1.0e-4 + end +end From 98bf3cbf6604786194cc25f6e544762c367e304a Mon Sep 17 00:00:00 2001 From: "Chris Rackauckas (Claude)" Date: Wed, 10 Jun 2026 06:13:43 -0400 Subject: [PATCH 2/2] Promote W to its solve-time type inside build_J_W Per review: fix the W typing at the source instead of patching the JacReuseState seed afterwards. build_J_W now computes a dtgamma prototype (dt's full possibly-Dual time type promoted with the constvalue-stripped tableau eltype) and constructs every non-DAE W with the type calc_W produces at solve time: - OOP static: StaticWOperator(J - mass_matrix * inv(oneunit(dtgamma))) - OOP dense: lu_instance(J - mass_matrix * inv(oneunit(dtgamma))) - WOperator: gamma slot seeded with the dtgamma prototype instead of promote(t, dt)[2] / raw dt This removes the build_jac_reuse_W_seed helper; the OOP Rosenbrock alg_caches pass build_J_W's W straight through again. It also fixes the same latent seed/solve-time W type mismatch for Newton-based OOP caches that store build_J_W's W and later assign calc_W results. DAE branches are unchanged: their W is not built from J - M/(dt*gamma). Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Fable 5 --- .../src/derivative_utils.jl | 78 +++++++++---------- .../src/OrdinaryDiffEqRosenbrock.jl | 1 - .../src/rosenbrock_caches.jl | 9 +-- 3 files changed, 40 insertions(+), 48 deletions(-) diff --git a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl index c562fd29447..700789bcacf 100644 --- a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl +++ b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl @@ -1072,6 +1072,21 @@ function update_W!( return nothing end +""" + _dtgamma_prototype(t, dt, uEltypeNoUnits) + +A value carrying the type `dtgamma = dt * γ` has at solve time: `dt`'s full +(possibly Dual) time type promoted with the tableau eltype +(`constvalue(uEltypeNoUnits)`, which `constvalue` strips of Duals the same way +tableau constructors do). Non-`Number` eltypes (e.g. array-of-array states) +have no tableau-eltype contribution to promote with, so `dt`'s type is used +as-is. +""" +@inline function _dtgamma_prototype(t, dt, ::Type{T}) where {T <: Number} + return promote(t, dt)[2] * one(constvalue(T)) +end +@inline _dtgamma_prototype(t, dt, ::Type{T}) where {T} = promote(t, dt)[2] + function build_J_W( alg, u, uprev, p, t, dt, f::F, jac_config, ::Type{uEltypeNoUnits}, ::Val{IIP} @@ -1079,6 +1094,16 @@ function build_J_W( # TODO - make J, W AbstractSciMLOperators (lazily defined with scimlops functionality) # TODO - if jvp given, make it SciMLOperators.FunctionOperator # TODO - make mass matrix a SciMLOperator so it can be updated with time. Default to IdentityOperator + # + # W must already carry the type `calc_W` produces at solve time: there + # W = J - mass_matrix * inv(dtgamma) promotes the eltype past the state + # eltype (e.g. a Float32 state over a Float64 tspan gives a Float64 W), + # and WOperator's gamma slot holds dtgamma itself. OOP caches store W + # concretely (e.g. JacReuseState.cached_W, Newton caches), so a + # state-eltype W here would reject calc_W's result on the first + # assignment. + dtgamma_prototype = _dtgamma_prototype(t, dt, uEltypeNoUnits) + invdtgamma_prototype = inv(oneunit(dtgamma_prototype)) islin, isode = islinearfunction(f, alg) if isdefined(f, :W_prototype) && (f.W_prototype isa AbstractSciMLOperator) # We use W_prototype when it is provided as a SciMLOperator, and in this case we require jac_prototype to be a SciMLOperator too. @@ -1093,10 +1118,10 @@ function build_J_W( @assert SciMLBase.has_jac(f) "f needs to have an associated jacobian" J = MatrixOperator(J; update_func! = f.jac) end - W = WOperator{IIP}(f.mass_matrix, promote(t, dt)[2], J, _vec(u)) + W = WOperator{IIP}(f.mass_matrix, dtgamma_prototype, J, _vec(u)) elseif islin J = isode ? f.f : f.f1.f # unwrap the Jacobian accordingly - W = WOperator{IIP}(f.mass_matrix, dt, J, _vec(u)) + W = WOperator{IIP}(f.mass_matrix, dtgamma_prototype, J, _vec(u)) elseif IIP && f.jac_prototype !== nothing && concrete_jac(alg) === nothing && (alg.linsolve === nothing || LinearSolve.needs_concrete_A(alg.linsolve)) @@ -1114,7 +1139,7 @@ function build_J_W( jacvec = JVPCache(f, copy(u), u, p, t, autodiff = alg_autodiff(alg)) J = jacvec - W = WOperator{IIP}(f.mass_matrix, promote(t, dt)[2], J, _vec(u), jacvec) + W = WOperator{IIP}(f.mass_matrix, dtgamma_prototype, J, _vec(u), jacvec) elseif alg.linsolve !== nothing && !LinearSolve.needs_concrete_A(alg.linsolve) || concrete_jac(alg) !== nothing && concrete_jac(alg) # The linear solver does not need a concrete Jacobian, but the user has @@ -1141,11 +1166,14 @@ function build_J_W( deepcopy(f.jac_prototype) end W = if J isa StaticMatrix - StaticWOperator(J, false) + # callinv = false skips inverting the seed-valued matrix while + # producing the same concrete type as calc_W's + # `StaticWOperator(J - mass_matrix * inv(dtgamma))`. + StaticWOperator(J - f.mass_matrix * invdtgamma_prototype, false) else jacvec = JVPCache(f, copy(u), u, p, t, autodiff = alg_autodiff(alg)) - WOperator{IIP}(f.mass_matrix, promote(t, dt)[2], J, _vec(u), jacvec) + WOperator{IIP}(f.mass_matrix, dtgamma_prototype, J, _vec(u), jacvec) end else J = if !IIP && SciMLBase.has_jac(f) @@ -1171,6 +1199,8 @@ function build_J_W( deepcopy(f.jac_prototype) end W = if alg isa DAEAlgorithm + # DAE W is not built from J - M/(dt·γ) (calc_W returns J or a + # cj-scaled combination), so no dtgamma promotion applies. if IIP similar(J) elseif J isa StaticMatrix @@ -1181,48 +1211,14 @@ function build_J_W( elseif IIP similar(J) elseif J isa StaticMatrix - StaticWOperator(J, false) + StaticWOperator(J - f.mass_matrix * invdtgamma_prototype, false) else - ArrayInterface.lu_instance(J) + ArrayInterface.lu_instance(J - f.mass_matrix * invdtgamma_prototype) end end return J, W end -""" - build_jac_reuse_W_seed(W, J, u, mass_matrix, dtgamma) - -Build the seed value for `JacReuseState`'s `cached_W` slot in OOP caches. - -The seed's *type* must match what `calc_W` produces at solve time, which -promotes the W eltype with `typeof(dtgamma)` and `eltype(mass_matrix)` — e.g. -a `Float32` state solved over a `Float64` time span yields a `Float64`-eltype -W. Seeding with `build_J_W`'s state-eltype `W` would make the -`jac_reuse.cached_W = W` assignment in `calc_rosenbrock_differentiation` throw -for wrapper types without an eltype-converting `convert` (`StaticWOperator`). - -`dtgamma` must be a prototype carrying the solve-time `dt * gamma` type; only -its type is used, so the seed's values are garbage. That is safe because the -"first use" guard in `_rosenbrock_jac_reuse_decision` guarantees the seed is -overwritten before it is ever read. -""" -function build_jac_reuse_W_seed(W, J, u, mass_matrix, dtgamma) - invdtgamma′ = inv(oneunit(dtgamma)) - if W isa StaticWOperator - # callinv = false skips inverting the garbage-valued seed while still - # producing the same concrete type as calc_W's - # `StaticWOperator(J - mass_matrix * inv(dtgamma))`. - return StaticWOperator(J - mass_matrix * invdtgamma′, false) - elseif W isa WOperator - return WOperator{false}(mass_matrix, dtgamma, J, u, W.jacvec) - elseif W isa AbstractSciMLOperator - return W - else - _W = J - mass_matrix * invdtgamma′ - return _W isa Number ? _W : ArrayInterface.lu_instance(_W) - end -end - build_uf(alg, nf, t, p, ::Val{true}) = UJacobianWrapper(nf, t, p) build_uf(alg, nf, t, p, ::Val{false}) = UDerivativeWrapper(nf, t, p) diff --git a/lib/OrdinaryDiffEqRosenbrock/src/OrdinaryDiffEqRosenbrock.jl b/lib/OrdinaryDiffEqRosenbrock/src/OrdinaryDiffEqRosenbrock.jl index 9b13e8b8cd9..50153bacbc4 100644 --- a/lib/OrdinaryDiffEqRosenbrock/src/OrdinaryDiffEqRosenbrock.jl +++ b/lib/OrdinaryDiffEqRosenbrock/src/OrdinaryDiffEqRosenbrock.jl @@ -33,7 +33,6 @@ using OrdinaryDiffEqDifferentiation: TimeDerivativeWrapper, TimeGradientWrapper, build_jac_config, issuccess_W, jacobian2W!, resize_jac_config!, resize_grad_config!, calc_W, calc_rosenbrock_differentiation!, build_J_W, - build_jac_reuse_W_seed, UJacobianWrapper, dolinsolve, WOperator, resize_J_W! using OrdinaryDiffEqDifferentiation: calc_rosenbrock_differentiation diff --git a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl index 8f0b9ef43a4..d78b1ce65e4 100644 --- a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl +++ b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl @@ -350,8 +350,7 @@ function alg_cache( return Rosenbrock23ConstantCache( tab.c₃₂, tab.d, tf, uf, J, W, linsolve, alg_autodiff(alg), _make_jac_reuse_state( - dtgamma_seed, alg.max_jac_age, J, zero(rate_prototype), - build_jac_reuse_W_seed(W, J, u, f.mass_matrix, dtgamma_seed) + dtgamma_seed, alg.max_jac_age, J, zero(rate_prototype), W ) ) end @@ -387,8 +386,7 @@ function alg_cache( return Rosenbrock32ConstantCache( tab.c₃₂, tab.d, tf, uf, J, W, linsolve, alg_autodiff(alg), _make_jac_reuse_state( - dtgamma_seed, alg.max_jac_age, J, zero(rate_prototype), - build_jac_reuse_W_seed(W, J, u, f.mass_matrix, dtgamma_seed) + dtgamma_seed, alg.max_jac_age, J, zero(rate_prototype), W ) ) end @@ -488,8 +486,7 @@ function alg_cache( tab, J, W, linsolve, alg_autodiff(alg), interp_order, _make_jac_reuse_state( - dtgamma_seed, alg.max_jac_age, J, zero(rate_prototype), - build_jac_reuse_W_seed(W, J, u, f.mass_matrix, dtgamma_seed) + dtgamma_seed, alg.max_jac_age, J, zero(rate_prototype), W ) ) end