Improve NULL FSAL and fix du for non-FSAL

Fixes the downstream SciMLSensitivity

Author: ChrisRackauckas

lib/OrdinaryDiffEqBDF/src/dae_caches.jl

@@ -15,7 +15,7 @@ end
 end
 
 # Not FSAL
-get_fsalfirstlast(cache::DImplicitEulerCache, u) = (u, u)
+get_fsalfirstlast(cache::DImplicitEulerCache, u) = (nothing, nothing)
 
 mutable struct DImplicitEulerConstantCache{N} <: OrdinaryDiffEqConstantCache
     nlsolver::N

lib/OrdinaryDiffEqCore/src/integrators/integrator_interface.jl

@@ -59,7 +59,7 @@ end
 @inline function DiffEqBase.get_du(integrator::ODEIntegrator)
     isdiscretecache(integrator.cache) &&
         error("Derivatives are not defined for this stepper.")
-    return if isdefined(integrator, :fsallast)
+    return if isfsal(integrator.alg)
         integrator.fsallast
     else
         integrator(integrator.t, Val{1})
@@ -72,7 +72,7 @@ end
     if isdiscretecache(integrator.cache)
         out .= integrator.cache.tmp
     else
-        return if isdefined(integrator, :fsallast) &&
+        return if isfsal(integrator.alg) &&
                   !has_stiff_interpolation(integrator.alg)
             # Special stiff interpolations do not store the
             # right value in fsallast

lib/OrdinaryDiffEqFunctionMap/src/functionmap_caches.jl

@@ -3,7 +3,7 @@
     uprev::uType
     tmp::rateType
 end
-get_fsalfirstlast(cache::FunctionMapCache, u) = (cache.u, cache.uprev)
+get_fsalfirstlast(cache::FunctionMapCache, u) = (nothing, nothing)
 
 function alg_cache(alg::FunctionMap, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,

lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_perform_step.jl

@@ -33,7 +33,7 @@ end
     integrator.u = u
 end
 
-get_fsalfirstlast(cache::LowStorageRK2NCache, u) = (cache.k, cache.k)
+get_fsalfirstlast(cache::LowStorageRK2NCache, u) = (nothing, nothing)
 
 function initialize!(integrator, cache::LowStorageRK2NCache)
     @unpack k, tmp, williamson_condition = cache

lib/OrdinaryDiffEqPDIRK/src/pdirk_caches.jl

@@ -8,7 +8,7 @@
 end
 
 # Non-FSAL
-get_fsalfirstlast(cache::PDIRK44Cache, u) = (cache.u, cache.uprev)
+get_fsalfirstlast(cache::PDIRK44Cache, u) = (nothing, nothing)
 
 struct PDIRK44ConstantCache{N, TabType} <: OrdinaryDiffEqConstantCache
     nlsolver::N

lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl

@@ -2,7 +2,7 @@ abstract type RosenbrockMutableCache <: OrdinaryDiffEqMutableCache end
 abstract type RosenbrockConstantCache <: OrdinaryDiffEqConstantCache end
 
 # Fake values since non-FSAL
-get_fsalfirstlast(cache::RosenbrockMutableCache, u) = (zero(u), zero(u))
+get_fsalfirstlast(cache::RosenbrockMutableCache, u) = (nothing, nothing)
 
 ################################################################################
 

lib/OrdinaryDiffEqVerner/src/verner_caches.jl

@@ -87,7 +87,7 @@ end
 end
 
 # fake values since non-FSAL method
-get_fsalfirstlast(cache::Vern7Cache, u) = (cache.k1, cache.k2)
+get_fsalfirstlast(cache::Vern7Cache, u) = (nothing, nothing)
 
 function alg_cache(alg::Vern7, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
@@ -153,7 +153,7 @@ end
 end
 
 # fake values since non-FSAL method
-get_fsalfirstlast(cache::Vern8Cache, u) = (cache.k1, cache.k2)
+get_fsalfirstlast(cache::Vern8Cache, u) = (nothing, nothing)
 
 function alg_cache(alg::Vern8, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
@@ -227,7 +227,7 @@ end
 end
 
 # fake values since non-FSAL method
-get_fsalfirstlast(cache::Vern9Cache, u) = (cache.k1, cache.k2)
+get_fsalfirstlast(cache::Vern9Cache, u) = (nothing, nothing)
 
 function alg_cache(alg::Vern9, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,

test/interface/get_du.jl

@@ -0,0 +1,26 @@
+using OrdinaryDiffEq, OrdinaryDiffEqCore, Test
+function lorenz!(du, u, p, t)
+    du[1] = 10.0(u[2] - u[1])
+    du[2] = u[1] * (28.0 - u[3]) - u[2]
+    du[3] = u[1] * u[2] - (8 / 3) * u[3]
+end
+u0 = [1.0; 0.0; 0.0]
+tspan = (0.0, 3.0)
+condition(u,t,integrator) = t == 0.2
+cache = zeros(3)
+affect!(integrator) = cache .= get_du(integrator)
+dusave = DiscreteCallback(condition, affect!)
+affect2!(integrator) = get_du!(cache, integrator)
+dusave_inplace = DiscreteCallback(condition, affect2!)
+
+prob = ODEProblem(lorenz!, u0, tspan)
+sol = solve(prob, Tsit5(), tstops = [0.2], callback = dusave, abstol=1e-12, reltol=1e-12)
+res = copy(cache)
+
+for alg in [Vern6(), Vern7(), Vern8(), Vern9(), Rodas4(), Rodas4P(), Rodas5(), Rodas5P(), TRBDF2(), KenCarp4(), FBDF(), QNDF()]
+    sol = solve(prob, alg, tstops = [0.2], callback = dusave, abstol=1e-12, reltol=1e-12)
+    @test res ≈ cache rtol=1e-5
+
+    sol = solve(prob, alg, tstops = [0.2], callback = dusave_inplace, abstol=1e-12, reltol=1e-12)
+    @test res ≈ cache rtol=1e-5
+end

test/runtests.jl

@@ -66,6 +66,7 @@ end
         @time @safetestset "Linear Solver Split ODE Tests" include("interface/linear_solver_split_ode_test.jl")
         @time @safetestset "Sparse Diff Tests" include("interface/sparsediff_tests.jl")
         @time @safetestset "Enum Tests" include("interface/enums.jl")
+        @time @safetestset "Enum Tests" include("interface/get_du.jl")
         @time @safetestset "Mass Matrix Tests" include("interface/mass_matrix_tests.jl")
         @time @safetestset "W-Operator prototype tests" include("interface/wprototype_tests.jl")
     end