Skip to content

Commit 7f6a1e8

Browse files
Harsh SinghHarsh Singh
authored andcommitted
fix(Rosenbrock): strip_cache override and Hermite first-derivative interpolant
- Add OrdinaryDiffEqCore.strip_cache(::RosenbrockCache) to fix InterfaceI failures: generic version passes all-Nothing args but concrete-typed fields (dense, dtC, dtd, ks, interp_order) reject Nothing - Add interp_order==-1 branch to all 4 Val{1} interpolant methods to fix Regression_I BoundsError: empty-H Hermite uses k[1]=f₀,k[2]=f₁ but old else-branch accessed k[3]
1 parent 2c5651c commit 7f6a1e8

File tree

2 files changed

+53
-4
lines changed

2 files changed

+53
-4
lines changed

lib/OrdinaryDiffEqRosenbrock/src/interp_func.jl

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,3 +27,31 @@ function SciMLBase.interp_summary(
2727
"specialized 4th (Rodas6P = 5th) order \"free\" stiffness-aware interpolation" :
2828
"1st order linear"
2929
end
30+
31+
# strip_cache for RosenbrockCache: the generic OrdinaryDiffEqCore version passes
32+
# all-Nothing args, but several fields (dense, dtC, dtd, ks, interp_order) have
33+
# concrete types that don't accept Nothing. Provide a custom override that clears
34+
# only the interpolation-related fields to zero-length arrays / nothing.
35+
function OrdinaryDiffEqCore.strip_cache(cache::RosenbrockCache)
36+
return RosenbrockCache(
37+
cache.u, cache.uprev,
38+
similar(cache.dense, 0), # dense::Vector{rateType} — must not be Nothing
39+
cache.du, cache.du1, cache.du2,
40+
similar(cache.dtC, 0, 0), # dtC::Matrix{tabType} — must not be Nothing
41+
similar(cache.dtd, 0), # dtd::Vector{tabType} — must not be Nothing
42+
similar(cache.ks, 0), # ks::Vector{rateType} — must not be Nothing
43+
cache.fsalfirst, cache.fsallast, cache.dT,
44+
cache.J, cache.W,
45+
cache.tmp, cache.atmp, cache.weight,
46+
cache.tab,
47+
nothing, # tf
48+
nothing, # uf
49+
cache.linsolve_tmp, cache.linsolve,
50+
nothing, # jac_config
51+
nothing, # grad_config
52+
cache.reltol, cache.alg,
53+
cache.step_limiter!, cache.stage_limiter!,
54+
cache.interp_order, # interp_order::Int — must not be Nothing
55+
cache.jac_reuse # jac_reuse — preserve state across strip
56+
)
57+
end

lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_interpolants.jl

Lines changed: 25 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -145,7 +145,13 @@ end
145145
},
146146
idxs::Nothing, T::Type{Val{1}}, differential_vars
147147
)
148-
if !hasproperty(cache, :interp_order) || cache.interp_order == 2
148+
if hasproperty(cache, :interp_order) && cache.interp_order == -1
149+
# Hermite interpolation: k[1]=f₀, k[2]=f₁
150+
# dH/dt = 6Θ(1-Θ)(y₁-y₀)/dt + f₀(3Θ²-4Θ+1) + f₁(3Θ²-2Θ)
151+
@.. 6 * Θ * (1 - Θ) * (y₁ - y₀) / dt +
152+
k[1] * (3 * Θ^2 - 4 * Θ + 1) +
153+
k[2] * (3 * Θ^2 - 2 * Θ)
154+
elseif !hasproperty(cache, :interp_order) || cache.interp_order == 2
149155
@.. (k[1] + Θ * (-2 * k[1] + 2 * k[2] - 3 * k[2] * Θ) - y₀ + y₁) / dt
150156
elseif cache.interp_order == 4
151157
@.. (
@@ -177,7 +183,12 @@ end
177183
},
178184
idxs, T::Type{Val{1}}, differential_vars
179185
)
180-
if !hasproperty(cache, :interp_order) || cache.interp_order == 2
186+
if hasproperty(cache, :interp_order) && cache.interp_order == -1
187+
# Hermite interpolation: k[1]=f₀, k[2]=f₁
188+
@views @.. 6 * Θ * (1 - Θ) * (y₁[idxs] - y₀[idxs]) / dt +
189+
k[1][idxs] * (3 * Θ^2 - 4 * Θ + 1) +
190+
k[2][idxs] * (3 * Θ^2 - 2 * Θ)
191+
elseif !hasproperty(cache, :interp_order) || cache.interp_order == 2
181192
@views @.. (
182193
k[1][idxs] +
183194
Θ * (-2 * k[1][idxs] + 2 * k[2][idxs] - 3 * k[2][idxs] * Θ) -
@@ -215,7 +226,12 @@ end
215226
},
216227
idxs::Nothing, T::Type{Val{1}}, differential_vars
217228
)
218-
if !hasproperty(cache, :interp_order) || cache.interp_order == 2
229+
if hasproperty(cache, :interp_order) && cache.interp_order == -1
230+
# Hermite interpolation: k[1]=f₀, k[2]=f₁
231+
@.. out = 6 * Θ * (1 - Θ) * (y₁ - y₀) / dt +
232+
k[1] * (3 * Θ^2 - 4 * Θ + 1) +
233+
k[2] * (3 * Θ^2 - 2 * Θ)
234+
elseif !hasproperty(cache, :interp_order) || cache.interp_order == 2
219235
@.. out = (k[1] + Θ * (-2 * k[1] + 2 * k[2] - 3 * k[2] * Θ) - y₀ + y₁) / dt
220236
elseif cache.interp_order == 4
221237
@.. out = (
@@ -249,7 +265,12 @@ end
249265
},
250266
idxs, T::Type{Val{1}}, differential_vars
251267
)
252-
if !hasproperty(cache, :interp_order) || cache.interp_order == 2
268+
if hasproperty(cache, :interp_order) && cache.interp_order == -1
269+
# Hermite interpolation: k[1]=f₀, k[2]=f₁
270+
@views @.. out = 6 * Θ * (1 - Θ) * (y₁[idxs] - y₀[idxs]) / dt +
271+
k[1][idxs] * (3 * Θ^2 - 4 * Θ + 1) +
272+
k[2][idxs] * (3 * Θ^2 - 2 * Θ)
273+
elseif !hasproperty(cache, :interp_order) || cache.interp_order == 2
253274
@views @.. out = (
254275
k[1][idxs] +
255276
Θ *

0 commit comments

Comments
 (0)