diff --git a/Project.toml b/Project.toml index 2b03aaaa35..6007276696 100644 --- a/Project.toml +++ b/Project.toml @@ -13,6 +13,7 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18" FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a" +FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" FunctionWrappersWrappers = "77dc65aa-8811-40c2-897b-53d922fa7daf" diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index 4cefb0aec1..92be5a0aeb 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -28,6 +28,8 @@ using LinearSolve, SimpleNonlinearSolve using LineSearches +import FillArrays: Trues + # Interfaces import DiffEqBase: solve!, step!, initialize!, isadaptive diff --git a/src/dense/generic_dense.jl b/src/dense/generic_dense.jl index a6205f41d6..c3ee351470 100644 --- a/src/dense/generic_dense.jl +++ b/src/dense/generic_dense.jl @@ -90,7 +90,7 @@ end DiffEqBase.addsteps!(integrator) if !(integrator.cache isa CompositeCache) val = ode_interpolant(Θ, integrator.dt, integrator.uprev, integrator.u, - integrator.k, integrator.cache, idxs, deriv) + integrator.k, integrator.cache, idxs, deriv, integrator.differential_vars) else val = composite_ode_interpolant(Θ, integrator, integrator.cache.caches, integrator.cache.current, idxs, deriv) @@ -107,7 +107,7 @@ end if $i == current return ode_interpolant(Θ, integrator.dt, integrator.uprev, integrator.u, integrator.k, caches[$i], idxs, - deriv) + deriv, integrator.differential_vars) end end) end @@ -122,11 +122,11 @@ end DiffEqBase.addsteps!(integrator) if !(integrator.cache isa CompositeCache) ode_interpolant!(val, Θ, integrator.dt, integrator.uprev, integrator.u, - integrator.k, integrator.cache, idxs, deriv) + integrator.k, integrator.cache, idxs, deriv, integrator.differential_vars) else ode_interpolant!(val, Θ, integrator.dt, integrator.uprev, integrator.u, integrator.k, integrator.cache.caches[integrator.cache.current], - idxs, deriv) + idxs, deriv, integrator.differential_vars) end end @@ -210,7 +210,7 @@ end DiffEqBase.addsteps!(integrator) if !(integrator.cache isa CompositeCache) ode_interpolant!(val, Θ, integrator.t - integrator.tprev, integrator.uprev2, - integrator.uprev, integrator.k, integrator.cache, idxs, deriv) + integrator.uprev, integrator.k, integrator.cache, idxs, deriv, integrator.differential_vars) else composite_ode_extrapolant!(val, Θ, integrator, integrator.cache.caches, integrator.cache.current, idxs, deriv) @@ -226,7 +226,7 @@ end if $i == current return ode_interpolant!(val, Θ, integrator.t - integrator.tprev, integrator.uprev2, integrator.uprev, - integrator.k, caches[$i], idxs, deriv) + integrator.k, caches[$i], idxs, deriv, integrator.differential_vars) end end) end @@ -241,7 +241,7 @@ end DiffEqBase.addsteps!(integrator) if !(integrator.cache isa CompositeCache) ode_interpolant(Θ, integrator.t - integrator.tprev, integrator.uprev2, - integrator.uprev, integrator.k, integrator.cache, idxs, deriv) + integrator.uprev, integrator.k, integrator.cache, idxs, deriv, integrator.differential_vars) else composite_ode_extrapolant(Θ, integrator, integrator.cache.caches, integrator.cache.current, idxs, deriv) @@ -257,7 +257,7 @@ end if $i == current return ode_interpolant(Θ, integrator.t - integrator.tprev, integrator.uprev2, integrator.uprev, - integrator.k, caches[$i], idxs, deriv) + integrator.k, caches[$i], idxs, deriv, integrator.differential_vars) end end) end @@ -270,46 +270,46 @@ end function _evaluate_interpolant(f, Θ, dt, timeseries, i₋, i₊, cache, idxs, - deriv, ks, ts, p) + deriv, ks, ts, p, differential_vars) _ode_addsteps!(ks[i₊], ts[i₋], timeseries[i₋], timeseries[i₊], dt, f, p, cache) # update the kcurrent return ode_interpolant(Θ, dt, timeseries[i₋], timeseries[i₊], ks[i₊], - cache, idxs, deriv) + cache, idxs, deriv, differential_vars) end function evaluate_composite_cache(f, Θ, dt, timeseries, i₋, i₊, caches::Tuple{C1, C2, Vararg}, idxs, - deriv, ks, ts, p, cacheid) where {C1, C2} + deriv, ks, ts, p, cacheid, differential_vars) where {C1, C2} if (cacheid -= 1) != 0 return evaluate_composite_cache(f, Θ, dt, timeseries, i₋, i₊, Base.tail(caches), idxs, - deriv, ks, ts, p, cacheid) + deriv, ks, ts, p, cacheid, differential_vars) end _evaluate_interpolant(f, Θ, dt, timeseries, i₋, i₊, first(caches), idxs, - deriv, ks, ts, p) + deriv, ks, ts, p, differential_vars) end function evaluate_composite_cache(f, Θ, dt, timeseries, i₋, i₊, caches::Tuple{C}, idxs, - deriv, ks, ts, p, _) where {C} + deriv, ks, ts, p, _, differential_vars) where {C} _evaluate_interpolant(f, Θ, dt, timeseries, i₋, i₊, only(caches), idxs, - deriv, ks, ts, p) + deriv, ks, ts, p, differential_vars) end function evaluate_interpolant(f, Θ, dt, timeseries, i₋, i₊, cache, idxs, - deriv, ks, ts, id, p) + deriv, ks, ts, id, p, differential_vars) if cache isa (FunctionMapCache) || cache isa FunctionMapConstantCache return ode_interpolant(Θ, dt, timeseries[i₋], timeseries[i₊], 0, cache, idxs, - deriv) + deriv, differential_vars) elseif !id.dense return linear_interpolant(Θ, dt, timeseries[i₋], timeseries[i₊], idxs, deriv) elseif cache isa CompositeCache return evaluate_composite_cache(f, Θ, dt, timeseries, i₋, i₊, cache.caches, idxs, - deriv, ks, ts, p, id.alg_choice[i₊]) + deriv, ks, ts, p, id.alg_choice[i₊], differential_vars) else return _evaluate_interpolant(f, Θ, dt, timeseries, i₋, i₊, cache, idxs, - deriv, ks, ts, p) + deriv, ks, ts, p, differential_vars) end end @@ -321,7 +321,7 @@ times ts (sorted), with values timeseries and derivatives ks """ function ode_interpolation(tvals, id::I, idxs, deriv::D, p, continuity::Symbol = :left) where {I, D} - @unpack ts, timeseries, ks, f, cache = id + @unpack ts, timeseries, ks, f, cache, differential_vars = id @inbounds tdir = sign(ts[end] - ts[1]) idx = sortperm(tvals, rev = tdir < 0) # start the search thinking it's ts[1]-ts[2] @@ -344,7 +344,7 @@ function ode_interpolation(tvals, id::I, idxs, deriv::D, p, dt = ts[i₊] - ts[i₋] Θ = iszero(dt) ? oneunit(t) / oneunit(dt) : (t - ts[i₋]) / dt evaluate_interpolant(f, Θ, dt, timeseries, i₋, i₊, cache, idxs, - deriv, ks, ts, id, p) + deriv, ks, ts, id, p, differential_vars) end invpermute!(vals, idx) DiffEqArray(vals, tvals) @@ -358,7 +358,7 @@ times ts (sorted), with values timeseries and derivatives ks """ function ode_interpolation!(vals, tvals, id::I, idxs, deriv::D, p, continuity::Symbol = :left) where {I, D} - @unpack ts, timeseries, ks, f, cache = id + @unpack ts, timeseries, ks, f, cache, differential_vars = id @inbounds tdir = sign(ts[end] - ts[1]) idx = sortperm(tvals, rev = tdir < 0) @@ -392,10 +392,10 @@ function ode_interpolation!(vals, tvals, id::I, idxs, deriv::D, p, if cache isa (FunctionMapCache) || cache isa FunctionMapConstantCache if eltype(vals) <: AbstractArray ode_interpolant!(vals[j], Θ, dt, timeseries[i₋], timeseries[i₊], 0, cache, - idxs, deriv) + idxs, deriv, differential_vars) else vals[j] = ode_interpolant(Θ, dt, timeseries[i₋], timeseries[i₊], 0, cache, - idxs, deriv) + idxs, deriv, differential_vars) end elseif !id.dense if eltype(vals) <: AbstractArray @@ -414,20 +414,20 @@ function ode_interpolation!(vals, tvals, id::I, idxs, deriv::D, p, cache_i₊) # update the kcurrent if eltype(vals) <: AbstractArray ode_interpolant!(vals[j], Θ, dt, timeseries[i₋], timeseries[i₊], ks[i₊], - cache_i₊, idxs, deriv) + cache_i₊, idxs, deriv, differential_vars) else vals[j] = ode_interpolant(Θ, dt, timeseries[i₋], timeseries[i₊], ks[i₊], - cache_i₊, idxs, deriv) + cache_i₊, idxs, deriv, differential_vars) end else _ode_addsteps!(ks[i₊], ts[i₋], timeseries[i₋], timeseries[i₊], dt, f, p, cache) # update the kcurrent if eltype(vals) <: AbstractArray ode_interpolant!(vals[j], Θ, dt, timeseries[i₋], timeseries[i₊], ks[i₊], - cache, idxs, deriv) + cache, idxs, deriv, differential_vars) else vals[j] = ode_interpolant(Θ, dt, timeseries[i₋], timeseries[i₊], ks[i₊], - cache, idxs, deriv) + cache, idxs, deriv, differential_vars) end end end @@ -443,7 +443,7 @@ times ts (sorted), with values timeseries and derivatives ks """ function ode_interpolation(tval::Number, id::I, idxs, deriv::D, p, continuity::Symbol = :left) where {I, D} - @unpack ts, timeseries, ks, f, cache = id + @unpack ts, timeseries, ks, f, cache, differential_vars = id @inbounds tdir = sign(ts[end] - ts[1]) if continuity === :left @@ -464,19 +464,19 @@ function ode_interpolation(tval::Number, id::I, idxs, deriv::D, p, if cache isa (FunctionMapCache) || cache isa FunctionMapConstantCache val = ode_interpolant(Θ, dt, timeseries[i₋], timeseries[i₊], 0, cache, idxs, - deriv) + deriv, differential_vars) elseif !id.dense val = linear_interpolant(Θ, dt, timeseries[i₋], timeseries[i₊], idxs, deriv) elseif cache isa CompositeCache _ode_addsteps!(ks[i₊], ts[i₋], timeseries[i₋], timeseries[i₊], dt, f, p, cache.caches[id.alg_choice[i₊]]) # update the kcurrent val = ode_interpolant(Θ, dt, timeseries[i₋], timeseries[i₊], ks[i₊], - cache.caches[id.alg_choice[i₊]], idxs, deriv) + cache.caches[id.alg_choice[i₊]], idxs, deriv, differential_vars) else _ode_addsteps!(ks[i₊], ts[i₋], timeseries[i₋], timeseries[i₊], dt, f, p, cache) # update the kcurrent val = ode_interpolant(Θ, dt, timeseries[i₋], timeseries[i₊], ks[i₊], cache, - idxs, deriv) + idxs, deriv, differential_vars) end end @@ -491,7 +491,7 @@ times ts (sorted), with values timeseries and derivatives ks """ function ode_interpolation!(out, tval::Number, id::I, idxs, deriv::D, p, continuity::Symbol = :left) where {I, D} - @unpack ts, timeseries, ks, f, cache = id + @unpack ts, timeseries, ks, f, cache, differential_vars = id @inbounds tdir = sign(ts[end] - ts[1]) if continuity === :left @@ -512,19 +512,19 @@ function ode_interpolation!(out, tval::Number, id::I, idxs, deriv::D, p, if cache isa (FunctionMapCache) || cache isa FunctionMapConstantCache ode_interpolant!(out, Θ, dt, timeseries[i₋], timeseries[i₊], 0, cache, idxs, - deriv) + deriv, differential_vars) elseif !id.dense linear_interpolant!(out, Θ, dt, timeseries[i₋], timeseries[i₊], idxs, deriv) elseif cache isa CompositeCache _ode_addsteps!(ks[i₊], ts[i₋], timeseries[i₋], timeseries[i₊], dt, f, p, cache.caches[id.alg_choice[i₊]]) # update the kcurrent ode_interpolant!(out, Θ, dt, timeseries[i₋], timeseries[i₊], ks[i₊], - cache.caches[id.alg_choice[i₊]], idxs, deriv) + cache.caches[id.alg_choice[i₊]], idxs, deriv, differential_vars) else _ode_addsteps!(ks[i₊], ts[i₋], timeseries[i₋], timeseries[i₊], dt, f, p, cache) # update the kcurrent ode_interpolant!(out, Θ, dt, timeseries[i₋], timeseries[i₊], ks[i₊], cache, - idxs, deriv) + idxs, deriv, differential_vars) end end @@ -554,15 +554,15 @@ end """ ode_interpolant and ode_interpolant! dispatch """ -function ode_interpolant(Θ, dt, y₀, y₁, k, cache, idxs, T::Type{Val{TI}}) where {TI} - _ode_interpolant(Θ, dt, y₀, y₁, k, cache, idxs, T) +function ode_interpolant(Θ, dt, y₀, y₁, k, cache, idxs, T::Type{Val{TI}}, differential_vars) where {TI} + _ode_interpolant(Θ, dt, y₀, y₁, k, cache, idxs, T, differential_vars) end function ode_interpolant(Θ, dt, y₀, y₁, k, cache::OrdinaryDiffEqMutableCache, idxs, - T::Type{Val{TI}}) where {TI} + T::Type{Val{TI}}, differential_vars) where {TI} if idxs isa Number || y₀ isa Union{Number, SArray} # typeof(y₀) can be these if saveidxs gives a single value - _ode_interpolant(Θ, dt, y₀, y₁, k, cache, idxs, T) + _ode_interpolant(Θ, dt, y₀, y₁, k, cache, idxs, T, differential_vars) elseif idxs isa Nothing if y₁ isa Array{<:Number} out = similar(y₁, eltype(first(y₁) * oneunit(Θ))) @@ -570,7 +570,7 @@ function ode_interpolant(Θ, dt, y₀, y₁, k, cache::OrdinaryDiffEqMutableCach else out = oneunit(Θ) .* y₁ end - _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache, idxs, T) + _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache, idxs, T, differential_vars) else if y₁ isa Array{<:Number} out = similar(y₁, eltype(first(y₁) * oneunit(Θ)), axes(idxs)) @@ -580,24 +580,81 @@ function ode_interpolant(Θ, dt, y₀, y₁, k, cache::OrdinaryDiffEqMutableCach else out = oneunit(Θ) .* y₁[idxs] end - _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache, idxs, T) + _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache, idxs, T, differential_vars) end end -function ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache, idxs, T::Type{Val{TI}}) where {TI} - _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache, idxs, T) +function ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache, idxs, T::Type{Val{TI}}, differential_vars) where {TI} + _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache, idxs, T, differential_vars) end ##################### Hermite Interpolants +const HERMITE_CASE_NOT_DEFINED_MESSAGE = """ + Hermite interpolation is not defined in this case. The Hermite interpolation + fallback only supports diagonal mass matrices. If you have a DAE with a + non-diagonal mass matrix, then the dense output is not supported with this + ODE solver. Either use a method which has a specialized interpolation, + such as Rodas5P, or use `dense=false` + + You can find the list of available DAE solvers with their documented interpolations at: + https://docs.sciml.ai/DiffEqDocs/stable/solvers/dae_solve/ + """ + +struct HermiteInterpolationNonDiagonalError <: Exception end + +function Base.showerror(io::IO, e::HermiteInterpolationNonDiagonalError) + print(io, HERMITE_CASE_NOT_DEFINED_MESSAGE) + println(io, TruncatedStacktraces.VERBOSE_MSG) +end + # If no dispatch found, assume Hermite -function _ode_interpolant(Θ, dt, y₀, y₁, k, cache, idxs, T::Type{Val{TI}}) where {TI} - hermite_interpolant(Θ, dt, y₀, y₁, k, Val{cache isa OrdinaryDiffEqMutableCache}, - idxs, T) +function _ode_interpolant(Θ, dt, y₀, y₁, k, cache, idxs, T::Type{Val{TI}}, differential_vars) where {TI} + differential_vars isa DifferentialVarsUndefined && throw(HermiteInterpolationNonDiagonalError()) + + differential_vars = if differential_vars === nothing + if y₀ isa Number + differential_vars = true + elseif idxs === nothing + differential_vars = Trues(size(y₀)) + elseif idxs isa Number + differential_vars = true + else + differential_vars = Trues(size(idxs)) + end + elseif idxs isa Number + differential_vars[idxs] + elseif idxs === nothing + differential_vars + else + @view differential_vars[idxs] + end + + hermite_interpolant(Θ, dt, y₀, y₁, k, Val{cache isa OrdinaryDiffEqMutableCache}, idxs, T, differential_vars) end -function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache, idxs, T::Type{Val{TI}}) where {TI} - hermite_interpolant!(out, Θ, dt, y₀, y₁, k, idxs, T) +function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache, idxs, T::Type{Val{TI}}, differential_vars) where {TI} + differential_vars isa DifferentialVarsUndefined && throw(HermiteInterpolationNonDiagonalError()) + + differential_vars = if differential_vars === nothing + if y₀ isa Number + differential_vars = true + elseif idxs === nothing + differential_vars = Trues(size(out)) + elseif idxs isa Number + differential_vars = true + else + differential_vars = Trues(size(idxs)) + end + elseif idxs isa Number + differential_vars[idxs] + elseif idxs === nothing + differential_vars + else + @view differential_vars[idxs] + end + + hermite_interpolant!(out, Θ, dt, y₀, y₁, k, idxs, T, differential_vars) end """ @@ -606,65 +663,65 @@ Hairer Norsett Wanner Solving Ordinary Differential Euations I - Nonstiff Proble Herimte Interpolation, chosen if no other dispatch for ode_interpolant """ @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{false}}, idxs::Nothing, - T::Type{Val{0}}) # Default interpolant is Hermite + T::Type{Val{0}}, differential_vars) # Default interpolant is Hermite #@.. broadcast=false (1-Θ)*y₀+Θ*y₁+Θ*(Θ-1)*((1-2Θ)*(y₁-y₀)+(Θ-1)*dt*k[1] + Θ*dt*k[2]) @inbounds (1 - Θ) * y₀ + Θ * y₁ + - Θ * (Θ - 1) * ((1 - 2Θ) * (y₁ - y₀) + (Θ - 1) * dt * k[1] + Θ * dt * k[2]) + differential_vars .* (Θ * (Θ - 1) * ((1 - 2Θ) * (y₁ - y₀) + (Θ - 1) * dt * k[1] + Θ * dt * k[2])) end @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{true}}, idxs::Nothing, - T::Type{Val{0}}) # Default interpolant is Hermite + T::Type{Val{0}}, differential_vars) # Default interpolant is Hermite #@.. broadcast=false (1-Θ)*y₀+Θ*y₁+Θ*(Θ-1)*((1-2Θ)*(y₁-y₀)+(Θ-1)*dt*k[1] + Θ*dt*k[2]) @inbounds @.. broadcast=false (1 - Θ)*y₀+Θ*y₁+ - Θ*(Θ-1)*((1 - 2Θ)*(y₁ - y₀)+(Θ-1)*dt*k[1]+Θ*dt*k[2]) + differential_vars*Θ*(Θ-1)*((1 - 2Θ)*(y₁ - y₀)+(Θ-1)*dt*k[1]+Θ*dt*k[2]) end @muladd function hermite_interpolant(Θ, dt, y₀::Array, y₁, k, ::Type{Val{true}}, - idxs::Nothing, T::Type{Val{0}}) # Default interpolant is Hermite + idxs::Nothing, T::Type{Val{0}}, differential_vars) # Default interpolant is Hermite out = similar(y₀) @inbounds @simd ivdep for i in eachindex(y₀) out[i] = (1 - Θ) * y₀[i] + Θ * y₁[i] + - Θ * (Θ - 1) * + differential_vars[i] * Θ * (Θ - 1) * ((1 - 2Θ) * (y₁[i] - y₀[i]) + (Θ - 1) * dt * k[1][i] + Θ * dt * k[2][i]) end end -@muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, cache, idxs, T::Type{Val{0}}) # Default interpolant is Hermite +@muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, cache, idxs, T::Type{Val{0}}, differential_vars) # Default interpolant is Hermite # return @.. broadcast=false (1-Θ)*y₀[idxs]+Θ*y₁[idxs]+Θ*(Θ-1)*((1-2Θ)*(y₁[idxs]-y₀[idxs])+(Θ-1)*dt*k[1][idxs] + Θ*dt*k[2][idxs]) return (1 - Θ) * y₀[idxs] + Θ * y₁[idxs] + - Θ * (Θ - 1) * + differential_vars .* (Θ * (Θ - 1) * ((1 - 2Θ) * (y₁[idxs] - y₀[idxs]) + (Θ - 1) * dt * k[1][idxs] + - Θ * dt * k[2][idxs]) + Θ * dt * k[2][idxs])) end -@muladd function hermite_interpolant!(out, Θ, dt, y₀, y₁, k, idxs::Nothing, T::Type{Val{0}}) # Default interpolant is Hermite +@muladd function hermite_interpolant!(out, Θ, dt, y₀, y₁, k, idxs::Nothing, T::Type{Val{0}}, differential_vars) # Default interpolant is Hermite @inbounds @.. broadcast=false out=(1 - Θ) * y₀ + Θ * y₁ + - Θ * (Θ - 1) * + differential_vars * Θ * (Θ - 1) * ((1 - 2Θ) * (y₁ - y₀) + (Θ - 1) * dt * k[1] + Θ * dt * k[2]) end @muladd function hermite_interpolant!(out::Array, Θ, dt, y₀, y₁, k, idxs::Nothing, - T::Type{Val{0}}) # Default interpolant is Hermite + T::Type{Val{0}}, differential_vars) # Default interpolant is Hermite @inbounds @simd ivdep for i in eachindex(out) out[i] = (1 - Θ) * y₀[i] + Θ * y₁[i] + - Θ * (Θ - 1) * + differential_vars[i] * Θ * (Θ - 1) * ((1 - 2Θ) * (y₁[i] - y₀[i]) + (Θ - 1) * dt * k[1][i] + Θ * dt * k[2][i]) end out end -@muladd function hermite_interpolant!(out, Θ, dt, y₀, y₁, k, idxs, T::Type{Val{0}}) # Default interpolant is Hermite +@muladd function hermite_interpolant!(out, Θ, dt, y₀, y₁, k, idxs, T::Type{Val{0}}, differential_vars) # Default interpolant is Hermite @views @.. broadcast=false out=(1 - Θ) * y₀[idxs] + Θ * y₁[idxs] + - Θ * (Θ - 1) * + differential_vars * Θ * (Θ - 1) * ((1 - 2Θ) * (y₁[idxs] - y₀[idxs]) + (Θ - 1) * dt * k[1][idxs] + Θ * dt * k[2][idxs]) end -@muladd function hermite_interpolant!(out::Array, Θ, dt, y₀, y₁, k, idxs, T::Type{Val{0}}) # Default interpolant is Hermite +@muladd function hermite_interpolant!(out::Array, Θ, dt, y₀, y₁, k, idxs, T::Type{Val{0}}, differential_vars) # Default interpolant is Hermite @inbounds for (j, i) in enumerate(idxs) out[j] = (1 - Θ) * y₀[i] + Θ * y₁[i] + - Θ * (Θ - 1) * + differential_vars[j] * Θ * (Θ - 1) * ((1 - 2Θ) * (y₁[i] - y₀[i]) + (Θ - 1) * dt * k[1][i] + Θ * dt * k[2][i]) end out @@ -674,62 +731,69 @@ end Herimte Interpolation, chosen if no other dispatch for ode_interpolant """ @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{false}}, idxs::Nothing, - T::Type{Val{1}}) # Default interpolant is Hermite + T::Type{Val{1}}, differential_vars) # Default interpolant is Hermite #@.. broadcast=false k[1] + Θ*(-4*dt*k[1] - 2*dt*k[2] - 6*y₀ + Θ*(3*dt*k[1] + 3*dt*k[2] + 6*y₀ - 6*y₁) + 6*y₁)/dt - @inbounds k[1] + + @inbounds (.!differential_vars).*(y₁ - y₀)/dt + differential_vars .*( + k[1] + Θ * (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ + - Θ * (3 * dt * k[1] + 3 * dt * k[2] + 6 * y₀ - 6 * y₁) + 6 * y₁) / dt + Θ * (3 * dt * k[1] + 3 * dt * k[2] + 6 * y₀ - 6 * y₁) + 6 * y₁) / dt) end @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{true}}, idxs::Nothing, - T::Type{Val{1}}) # Default interpolant is Hermite - @inbounds @.. broadcast=false k[1]+Θ * (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ + + T::Type{Val{1}}, differential_vars) # Default interpolant is Hermite + @inbounds @.. broadcast=false !differential_vars*((y₁ - y₀)/dt)+differential_vars*( + k[1]+Θ * (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ + Θ * (3 * dt * k[1] + 3 * dt * k[2] + 6 * y₀ - 6 * y₁) + - 6 * y₁) / dt + 6 * y₁) / dt) end -@muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, cache, idxs, T::Type{Val{1}}) # Default interpolant is Hermite +@muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, cache, idxs, T::Type{Val{1}}, differential_vars) # Default interpolant is Hermite # return @.. broadcast=false k[1][idxs] + Θ*(-4*dt*k[1][idxs] - 2*dt*k[2][idxs] - 6*y₀[idxs] + Θ*(3*dt*k[1][idxs] + 3*dt*k[2][idxs] + 6*y₀[idxs] - 6*y₁[idxs]) + 6*y₁[idxs])/dt - return k[1][idxs] + + return (.!differential_vars).*((y₁[idxs] - y₀[idxs])/dt)+differential_vars.*( + k[1][idxs] + Θ * (-4 * dt * k[1][idxs] - 2 * dt * k[2][idxs] - 6 * y₀[idxs] + Θ * (3 * dt * k[1][idxs] + 3 * dt * k[2][idxs] + 6 * y₀[idxs] - 6 * y₁[idxs]) + - 6 * y₁[idxs]) / dt + 6 * y₁[idxs]) / dt) end -@muladd function hermite_interpolant!(out, Θ, dt, y₀, y₁, k, idxs::Nothing, T::Type{Val{1}}) # Default interpolant is Hermite - @inbounds @.. broadcast=false out=k[1] + +@muladd function hermite_interpolant!(out, Θ, dt, y₀, y₁, k, idxs::Nothing, T::Type{Val{1}}, differential_vars) # Default interpolant is Hermite + @inbounds @.. broadcast=false out=!differential_vars*((y₁ - y₀)/dt)+differential_vars*( + k[1] + Θ * (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ + Θ * (3 * dt * k[1] + 3 * dt * k[2] + 6 * y₀ - 6 * y₁) + - 6 * y₁) / dt + 6 * y₁) / dt) end @muladd function hermite_interpolant!(out::Array, Θ, dt, y₀, y₁, k, idxs::Nothing, - T::Type{Val{1}}) # Default interpolant is Hermite + T::Type{Val{1}}, differential_vars) # Default interpolant is Hermite @inbounds @simd ivdep for i in eachindex(out) - out[i] = k[1][i] + + out[i] = !differential_vars[i]*((y₁[i] - y₀[i])/dt)+differential_vars[i]*( + k[1][i] + Θ * (-4 * dt * k[1][i] - 2 * dt * k[2][i] - 6 * y₀[i] + Θ * (3 * dt * k[1][i] + 3 * dt * k[2][i] + 6 * y₀[i] - 6 * y₁[i]) + - 6 * y₁[i]) / dt + 6 * y₁[i]) / dt) end out end -@muladd function hermite_interpolant!(out, Θ, dt, y₀, y₁, k, idxs, T::Type{Val{1}}) # Default interpolant is Hermite - @views @.. broadcast=false out=k[1][idxs] + +@muladd function hermite_interpolant!(out, Θ, dt, y₀, y₁, k, idxs, T::Type{Val{1}}, differential_vars) # Default interpolant is Hermite + @views @.. broadcast=false out=!differential_vars*((y₁ - y₀)/dt)+differential_vars*( + k[1][idxs] + Θ * (-4 * dt * k[1][idxs] - 2 * dt * k[2][idxs] - 6 * y₀[idxs] + Θ * (3 * dt * k[1][idxs] + 3 * dt * k[2][idxs] + - 6 * y₀[idxs] - 6 * y₁[idxs]) + 6 * y₁[idxs]) / dt + 6 * y₀[idxs] - 6 * y₁[idxs]) + 6 * y₁[idxs]) / dt) end -@muladd function hermite_interpolant!(out::Array, Θ, dt, y₀, y₁, k, idxs, T::Type{Val{1}}) # Default interpolant is Hermite +@muladd function hermite_interpolant!(out::Array, Θ, dt, y₀, y₁, k, idxs, T::Type{Val{1}}, differential_vars) # Default interpolant is Hermite @inbounds for (j, i) in enumerate(idxs) - out[j] = k[1][i] + + out[j] = !differential_vars[j]*((y₁[i] - y₀[i])/dt)+differential_vars[j]*( + k[1][i] + Θ * (-4 * dt * k[1][i] - 2 * dt * k[2][i] - 6 * y₀[i] + Θ * (3 * dt * k[1][i] + 3 * dt * k[2][i] + 6 * y₀[i] - 6 * y₁[i]) + - 6 * y₁[i]) / dt + 6 * y₁[i]) / dt) end out end @@ -738,62 +802,57 @@ end Herimte Interpolation, chosen if no other dispatch for ode_interpolant """ @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{false}}, idxs::Nothing, - T::Type{Val{2}}) # Default interpolant is Hermite + T::Type{Val{2}}, differential_vars) # Default interpolant is Hermite #@.. broadcast=false (-4*dt*k[1] - 2*dt*k[2] - 6*y₀ + Θ*(6*dt*k[1] + 6*dt*k[2] + 12*y₀ - 12*y₁) + 6*y₁)/(dt*dt) - @inbounds (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ + + @inbounds differential_vars .* (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ + Θ * (6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - 12 * y₁) + 6 * y₁) / (dt * dt) end @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{true}}, idxs::Nothing, - T::Type{Val{2}}) # Default interpolant is Hermite + T::Type{Val{2}}, differential_vars) # Default interpolant is Hermite #@.. broadcast=false (-4*dt*k[1] - 2*dt*k[2] - 6*y₀ + Θ*(6*dt*k[1] + 6*dt*k[2] + 12*y₀ - 12*y₁) + 6*y₁)/(dt*dt) - @inbounds @.. broadcast=false (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ + + @inbounds @.. broadcast=false differential_vars * (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ + Θ * (6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - 12 * y₁) + 6 * y₁)/(dt * dt) end -@muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, cache, idxs, T::Type{Val{2}}) # Default interpolant is Hermite +@muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, cache, idxs, T::Type{Val{2}}, differential_vars) # Default interpolant is Hermite #out = similar(y₀,axes(idxs)) #@views @.. broadcast=false out = (-4*dt*k[1][idxs] - 2*dt*k[2][idxs] - 6*y₀[idxs] + Θ*(6*dt*k[1][idxs] + 6*dt*k[2][idxs] + 12*y₀[idxs] - 12*y₁[idxs]) + 6*y₁[idxs])/(dt*dt) - @views out = (-4 * dt * k[1][idxs] - 2 * dt * k[2][idxs] - 6 * y₀[idxs] + + @views out = differential_vars .* (-4 * dt * k[1][idxs] - 2 * dt * k[2][idxs] - 6 * y₀[idxs] + Θ * (6 * dt * k[1][idxs] + 6 * dt * k[2][idxs] + 12 * y₀[idxs] - 12 * y₁[idxs]) + 6 * y₁[idxs]) / (dt * dt) out end -@muladd function hermite_interpolant!(out, Θ, dt, y₀, y₁, k, idxs::Nothing, T::Type{Val{2}}) # Default interpolant is Hermite - @inbounds @.. broadcast=false out=(-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ + +@muladd function hermite_interpolant!(out, Θ, dt, y₀, y₁, k, idxs::Nothing, T::Type{Val{2}}, differential_vars) # Default interpolant is Hermite + @inbounds @.. broadcast=false out= differential_vars * (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ + Θ * (6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - 12 * y₁) + 6 * y₁) / (dt * dt) end @muladd function hermite_interpolant!(out::Array, Θ, dt, y₀, y₁, k, idxs::Nothing, - T::Type{Val{2}}) # Default interpolant is Hermite + T::Type{Val{2}}, differential_vars) # Default interpolant is Hermite @inbounds @simd ivdep for i in eachindex(out) - out[i] = (-4 * dt * k[1][i] - 2 * dt * k[2][i] - 6 * y₀[i] + + out[i] = differential_vars[i] * (-4 * dt * k[1][i] - 2 * dt * k[2][i] - 6 * y₀[i] + Θ * (6 * dt * k[1][i] + 6 * dt * k[2][i] + 12 * y₀[i] - 12 * y₁[i]) + 6 * y₁[i]) / (dt * dt) end out end -@muladd function hermite_interpolant!(out, Θ, dt, y₀, y₁, k, idxs, T::Type{Val{2}}) # Default interpolant is Hermite - @views @.. broadcast=false out=(-4 * dt * k[1][idxs] - 2 * dt * k[2][idxs] - +@muladd function hermite_interpolant!(out, Θ, dt, y₀, y₁, k, idxs, T::Type{Val{2}}, differential_vars) # Default interpolant is Hermite + @views @.. broadcast=false out=differential_vars * (-4 * dt * k[1][idxs] - 2 * dt * k[2][idxs] - 6 * y₀[idxs] + Θ * (6 * dt * k[1][idxs] + 6 * dt * k[2][idxs] + 12 * y₀[idxs] - 12 * y₁[idxs]) + 6 * y₁[idxs]) / (dt * dt) end -@muladd function hermite_interpolant!(out::Array, Θ, dt, y₀, y₁, k, idxs, T::Type{Val{2}}) # Default interpolant is Hermite - @views @.. broadcast=false out=(-4 * dt * k[1][idxs] - 2 * dt * k[2][idxs] - - 6 * y₀[idxs] + - Θ * (6 * dt * k[1][idxs] + 6 * dt * k[2][idxs] + - 12 * y₀[idxs] - 12 * y₁[idxs]) + 6 * y₁[idxs]) / - (dt * dt) +@muladd function hermite_interpolant!(out::Array, Θ, dt, y₀, y₁, k, idxs, T::Type{Val{2}}, differential_vars) # Default interpolant is Hermite @inbounds for (j, i) in enumerate(idxs) - out[j] = (-4 * dt * k[1][i] - 2 * dt * k[2][i] - 6 * y₀[i] + + out[j] = differential_vars[j] * (-4 * dt * k[1][i] - 2 * dt * k[2][i] - 6 * y₀[i] + Θ * (6 * dt * k[1][i] + 6 * dt * k[2][i] + 12 * y₀[i] - 12 * y₁[i]) + 6 * y₁[i]) / (dt * dt) end @@ -804,31 +863,31 @@ end Herimte Interpolation, chosen if no other dispatch for ode_interpolant """ @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{false}}, idxs::Nothing, - T::Type{Val{3}}) # Default interpolant is Hermite + T::Type{Val{3}}, differential_vars) # Default interpolant is Hermite #@.. broadcast=false (6*dt*k[1] + 6*dt*k[2] + 12*y₀ - 12*y₁)/(dt*dt*dt) - @inbounds (6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - 12 * y₁) / (dt * dt * dt) + @inbounds differential_vars .* (6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - 12 * y₁) / (dt * dt * dt) end @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{true}}, idxs::Nothing, - T::Type{Val{3}}) # Default interpolant is Hermite + T::Type{Val{3}}, differential_vars) # Default interpolant is Hermite #@.. broadcast=false (6*dt*k[1] + 6*dt*k[2] + 12*y₀ - 12*y₁)/(dt*dt*dt) - @inbounds @.. broadcast=false (6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - + @inbounds @.. broadcast=false differential_vars * (6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - 12 * y₁)/(dt * dt * dt) end -@muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, cache, idxs, T::Type{Val{3}}) # Default interpolant is Hermite +@muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, cache, idxs, T::Type{Val{3}}, differential_vars) # Default interpolant is Hermite #out = similar(y₀,axes(idxs)) #@views @.. broadcast=false out = (6*dt*k[1][idxs] + 6*dt*k[2][idxs] + 12*y₀[idxs] - 12*y₁[idxs])/(dt*dt*dt) - @views out = (6 * dt * k[1][idxs] + 6 * dt * k[2][idxs] + 12 * y₀[idxs] - + @views out = differential_vars .* (6 * dt * k[1][idxs] + 6 * dt * k[2][idxs] + 12 * y₀[idxs] - 12 * y₁[idxs]) / (dt * dt * dt) out end -@muladd function hermite_interpolant!(out, Θ, dt, y₀, y₁, k, idxs::Nothing, T::Type{Val{3}}) # Default interpolant is Hermite - @inbounds @.. broadcast=false out=(6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - 12 * y₁) / +@muladd function hermite_interpolant!(out, Θ, dt, y₀, y₁, k, idxs::Nothing, T::Type{Val{3}}, differential_vars) # Default interpolant is Hermite + @inbounds @.. broadcast=false out=differential_vars * (6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - 12 * y₁) / (dt * dt * dt) #for i in eachindex(out) # out[i] = (6*dt*k[1][i] + 6*dt*k[2][i] + 12*y₀[i] - 12*y₁[i])/(dt*dt*dt) @@ -837,22 +896,22 @@ end end @muladd function hermite_interpolant!(out::Array, Θ, dt, y₀, y₁, k, idxs::Nothing, - T::Type{Val{3}}) # Default interpolant is Hermite + T::Type{Val{3}}, differential_vars) # Default interpolant is Hermite @inbounds @simd ivdep for i in eachindex(out) - out[i] = (6 * dt * k[1][i] + 6 * dt * k[2][i] + 12 * y₀[i] - 12 * y₁[i]) / + out[i] = differential_vars[i] * (6 * dt * k[1][i] + 6 * dt * k[2][i] + 12 * y₀[i] - 12 * y₁[i]) / (dt * dt * dt) end out end -@muladd function hermite_interpolant!(out, Θ, dt, y₀, y₁, k, idxs, T::Type{Val{3}}) # Default interpolant is Hermite +@muladd function hermite_interpolant!(out, Θ, dt, y₀, y₁, k, idxs, T::Type{Val{3}}, differential_vars) # Default interpolant is Hermite @views @.. broadcast=false out=(6 * dt * k[1][idxs] + 6 * dt * k[2][idxs] + 12 * y₀[idxs] - 12 * y₁[idxs]) / (dt * dt * dt) end -@muladd function hermite_interpolant!(out::Array, Θ, dt, y₀, y₁, k, idxs, T::Type{Val{3}}) # Default interpolant is Hermite +@muladd function hermite_interpolant!(out::Array, Θ, dt, y₀, y₁, k, idxs, T::Type{Val{3}}, differential_vars) # Default interpolant is Hermite @inbounds for (j, i) in enumerate(idxs) - out[j] = (6 * dt * k[1][i] + 6 * dt * k[2][i] + 12 * y₀[i] - 12 * y₁[i]) / + out[j] = differential_vars[j] * (6 * dt * k[1][i] + 6 * dt * k[2][i] + 12 * y₀[i] - 12 * y₁[i]) / (dt * dt * dt) end out diff --git a/src/dense/interpolants.jl b/src/dense/interpolants.jl index 722d685130..6e1cfbc6cf 100644 --- a/src/dense/interpolants.jl +++ b/src/dense/interpolants.jl @@ -1,12 +1,12 @@ @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{FunctionMapConstantCache, FunctionMapCache}, - idxs, T::Type{Val{0}}) + idxs, T::Type{Val{0}}, differential_vars::Nothing) y₀ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{FunctionMapConstantCache, FunctionMapCache}, - idxs, T::Type{Val{0}}) + idxs, T::Type{Val{0}}, differential_vars::Nothing) recursivecopy!(out, y₀) end @@ -21,13 +21,13 @@ Hairer Norsett Wanner Solving Ordinary Differential Euations I - Nonstiff Proble end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::DP5ConstantCache, idxs::Nothing, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @dp5pre0 @inbounds y₀ + dt * (k[1] * b10 + k[2] * b20 + k[3] * b30 + k[4] * b40) end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::DP5Cache, idxs::Nothing, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @dp5pre0 @inbounds @.. broadcast=false y₀+dt * (k[1] * b10 + k[2] * b20 + k[3] * b30 + k[4] * b40) @@ -35,7 +35,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{DP5ConstantCache, DP5Cache}, idxs, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @dp5pre0 @views @.. broadcast=false y₀[idxs]+dt * (k[1][idxs] * b10 + k[2][idxs] * b20 + k[3][idxs] * b30 + k[4][idxs] * b40) @@ -43,7 +43,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{DP5ConstantCache, DP5Cache}, idxs::Nothing, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @dp5pre0 @inbounds @.. broadcast=false out=y₀ + dt * @@ -53,7 +53,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{DP5ConstantCache, DP5Cache}, idxs, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @dp5pre0 @views @.. broadcast=false out=y₀[idxs] + dt * @@ -70,14 +70,14 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{DP5ConstantCache, DP5Cache}, idxs::Nothing, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars::Nothing) @dp5pre1 @inbounds @.. broadcast=false k[1]+k[2]*b20diff+k[3]*b30diff+k[4]*b40diff end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{DP5ConstantCache, DP5Cache}, idxs, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars::Nothing) @dp5pre1 @views @.. broadcast=false k[1][idxs]+k[2][idxs]*b20diff+k[3][idxs]*b30diff+ k[4][idxs]*b40diff @@ -85,7 +85,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{DP5ConstantCache, DP5Cache}, idxs::Nothing, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars::Nothing) @dp5pre1 @inbounds @.. broadcast=false out=k[1] + k[2] * b20diff + k[3] * b30diff + k[4] * b40diff @@ -94,7 +94,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{DP5ConstantCache, DP5Cache}, idxs, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars::Nothing) @dp5pre1 @views @.. broadcast=false out=k[1][idxs] + k[2][idxs] * b20diff + k[3][idxs] * b30diff + k[4][idxs] * b40diff @@ -110,7 +110,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{DP5ConstantCache, DP5Cache}, idxs::Nothing, - T::Type{Val{2}}) + T::Type{Val{2}}, differential_vars::Nothing) @dp5pre2 @inbounds @.. broadcast=false (k[2] * b20diff2 + k[3] * b30diff2 + k[4] * b40diff2)*invdt @@ -118,7 +118,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{DP5ConstantCache, DP5Cache}, idxs, - T::Type{Val{2}}) + T::Type{Val{2}}, differential_vars::Nothing) @dp5pre2 @views @.. broadcast=false (k[2][idxs] * b20diff2 + k[3][idxs] * b30diff2 + k[4][idxs] * b40diff2)*invdt @@ -126,7 +126,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{DP5ConstantCache, DP5Cache}, idxs::Nothing, - T::Type{Val{2}}) + T::Type{Val{2}}, differential_vars::Nothing) @dp5pre2 @inbounds @.. broadcast=false out=(k[2] * b20diff2 + k[3] * b30diff2 + k[4] * b40diff2) * @@ -136,7 +136,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{DP5ConstantCache, DP5Cache}, idxs, - T::Type{Val{2}}) + T::Type{Val{2}}, differential_vars::Nothing) @dp5pre2 @views @.. broadcast=false out=(k[2][idxs] * b20diff2 + k[3][idxs] * b30diff2 + k[4][idxs] * b40diff2) * invdt @@ -151,21 +151,21 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{DP5ConstantCache, DP5Cache}, idxs::Nothing, - T::Type{Val{3}}) + T::Type{Val{3}}, differential_vars::Nothing) @dp5pre3 @inbounds @.. broadcast=false (k[3] * b30diff3 + k[4] * b40diff3)*invdt2 end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{DP5ConstantCache, DP5Cache}, idxs, - T::Type{Val{3}}) + T::Type{Val{3}}, differential_vars::Nothing) @dp5pre3 @views @.. broadcast=false (k[3][idxs] * b30diff3 + k[4][idxs] * b40diff3)*invdt2 end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{DP5ConstantCache, DP5Cache}, idxs::Nothing, - T::Type{Val{3}}) + T::Type{Val{3}}, differential_vars::Nothing) @dp5pre3 @inbounds @.. broadcast=false out=(k[3] * b30diff3 + k[4] * b40diff3) * invdt2 out @@ -173,7 +173,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{DP5ConstantCache, DP5Cache}, idxs, - T::Type{Val{3}}) + T::Type{Val{3}}, differential_vars::Nothing) @dp5pre3 @views @.. broadcast=false out=(k[3][idxs] * b30diff3 + k[4][idxs] * b40diff3) * invdt2 out @@ -185,21 +185,21 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{DP5ConstantCache, DP5Cache}, idxs::Nothing, - T::Type{Val{4}}) + T::Type{Val{4}}, differential_vars::Nothing) @dp5pre4 @inbounds @.. broadcast=false k[4]*b40diff4invdt3 end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{DP5ConstantCache, DP5Cache}, idxs, - T::Type{Val{4}}) + T::Type{Val{4}}, differential_vars::Nothing) @dp5pre4 @views @.. broadcast=false k[4][idxs]*b40diff4invdt3 end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{DP5ConstantCache, DP5Cache}, idxs::Nothing, - T::Type{Val{4}}) + T::Type{Val{4}}, differential_vars::Nothing) @dp5pre4 @inbounds @.. broadcast=false out=k[4] * b40diff4invdt3 out @@ -207,7 +207,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{DP5ConstantCache, DP5Cache}, idxs, - T::Type{Val{4}}) + T::Type{Val{4}}, differential_vars::Nothing) @dp5pre4 @views @.. broadcast=false out=k[4][idxs] * b40diff4invdt3 out @@ -230,7 +230,7 @@ end SSPRK33ConstantCache, SSPRK33Cache, SSPRK43ConstantCache, SSPRK43Cache, SSPRK432ConstantCache, SSPRK432Cache}, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) @ssprkpre0 @inbounds @.. broadcast=false y₀*c00+y₁*c10+k[1]*b10dt end @@ -240,7 +240,7 @@ end SSPRK33ConstantCache, SSPRK33Cache, SSPRK43ConstantCache, SSPRK43Cache, SSPRK432ConstantCache, SSPRK432Cache}, idxs, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @ssprkpre0 @views @.. broadcast=false y₀[idxs]*c00+y₁[idxs]*c10+k[1][idxs]*b10dt end @@ -250,7 +250,7 @@ end SSPRK33ConstantCache, SSPRK33Cache, SSPRK43ConstantCache, SSPRK43Cache, SSPRK432ConstantCache, SSPRK432Cache}, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) @ssprkpre0 @inbounds @.. broadcast=false out=y₀ * c00 + y₁ * c10 + k[1] * b10dt out @@ -261,7 +261,7 @@ end SSPRK33ConstantCache, SSPRK33Cache, SSPRK43ConstantCache, SSPRK43Cache, SSPRK432ConstantCache, SSPRK432Cache}, idxs, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @ssprkpre0 @views @.. broadcast=false out=y₀[idxs] * c00 + y₁[idxs] * c10 + k[1][idxs] * b10dt out @@ -277,7 +277,7 @@ end SSPRK33ConstantCache, SSPRK33Cache, SSPRK43ConstantCache, SSPRK43Cache, SSPRK432ConstantCache, SSPRK432Cache}, - idxs::Nothing, T::Type{Val{1}}) + idxs::Nothing, T::Type{Val{1}}, differential_vars::Nothing) @ssprkpre1 @inbounds @.. broadcast=false (y₁ - y₀) * c10diffinvdt+k[1] * b10diff end @@ -287,7 +287,7 @@ end SSPRK33ConstantCache, SSPRK33Cache, SSPRK43ConstantCache, SSPRK43Cache, SSPRK432ConstantCache, SSPRK432Cache}, idxs, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars::Nothing) @ssprkpre1 @views @.. broadcast=false (y₁[idxs] - y₀[idxs]) * c10diffinvdt+k[1][idxs] * b10diff end @@ -297,7 +297,7 @@ end SSPRK33ConstantCache, SSPRK33Cache, SSPRK43ConstantCache, SSPRK43Cache, SSPRK432ConstantCache, SSPRK432Cache}, - idxs::Nothing, T::Type{Val{1}}) + idxs::Nothing, T::Type{Val{1}}, differential_vars::Nothing) @ssprkpre1 @inbounds @.. broadcast=false out=(y₁ - y₀) * c10diffinvdt + k[1] * b10diff out @@ -308,7 +308,7 @@ end SSPRK33ConstantCache, SSPRK33Cache, SSPRK43ConstantCache, SSPRK43Cache, SSPRK432ConstantCache, SSPRK432Cache}, idxs, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars::Nothing) @ssprkpre1 @views @.. broadcast=false out=(y₁[idxs] - y₀[idxs]) * c10diffinvdt + k[1][idxs] * b10diff @@ -326,7 +326,7 @@ end SSPRK33ConstantCache, SSPRK33Cache, SSPRK43ConstantCache, SSPRK43Cache, SSPRK432ConstantCache, SSPRK432Cache}, - idxs::Nothing, T::Type{Val{2}}) + idxs::Nothing, T::Type{Val{2}}, differential_vars::Nothing) @ssprkpre2 @inbounds @.. broadcast=false (y₁ - y₀) * c10diff2invdt2+k[1] * b10diff2invdt end @@ -336,7 +336,7 @@ end SSPRK33ConstantCache, SSPRK33Cache, SSPRK43ConstantCache, SSPRK43Cache, SSPRK432ConstantCache, SSPRK432Cache}, idxs, - T::Type{Val{2}}) + T::Type{Val{2}}, differential_vars::Nothing) @ssprkpre2 @views @.. broadcast=false (y₁[idxs] - y₀[idxs]) * c10diff2invdt2+k[1][idxs] * b10diff2invdt @@ -347,7 +347,7 @@ end SSPRK33ConstantCache, SSPRK33Cache, SSPRK43ConstantCache, SSPRK43Cache, SSPRK432ConstantCache, SSPRK432Cache}, - idxs::Nothing, T::Type{Val{2}}) + idxs::Nothing, T::Type{Val{2}}, differential_vars::Nothing) @ssprkpre2 @inbounds @.. broadcast=false out=(y₁ - y₀) * c10diff2invdt2 + k[1] * b10diff2invdt out @@ -358,7 +358,7 @@ end SSPRK33ConstantCache, SSPRK33Cache, SSPRK43ConstantCache, SSPRK43Cache, SSPRK432ConstantCache, SSPRK432Cache}, idxs, - T::Type{Val{2}}) + T::Type{Val{2}}, differential_vars::Nothing) @ssprkpre2 @views @.. broadcast=false out=(y₁[idxs] - y₀[idxs]) * c10diff2invdt2 + k[1][idxs] * b10diff2invdt @@ -389,7 +389,7 @@ end end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Tsit5ConstantCache, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) @tsit5pre0 #@.. broadcast=false y₀ + dt*(k[1]*b1Θ + k[2]*b2Θ + k[3]*b3Θ + k[4]*b4Θ + k[5]*b5Θ + k[6]*b6Θ + k[7]*b7Θ) return @inbounds y₀ + @@ -398,7 +398,7 @@ end end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Tsit5Cache, idxs::Nothing, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @tsit5pre0 return @inbounds @.. broadcast=false y₀+dt * (k[1] * b1Θ + k[2] * b2Θ + k[3] * b3Θ + k[4] * b4Θ + @@ -407,7 +407,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Tsit5ConstantCache, Tsit5Cache}, idxs, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @tsit5pre0 return y₀[idxs] + dt * (k[1][idxs] * b1Θ + k[2][idxs] * b2Θ + k[3][idxs] * b3Θ + @@ -416,7 +416,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Tsit5ConstantCache, Tsit5Cache}, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) @tsit5pre0 @inbounds @.. broadcast=false out=y₀ + dt * @@ -427,7 +427,7 @@ end @muladd function _ode_interpolant!(out::Array, Θ, dt, y₀, y₁, k, cache::Union{Tsit5ConstantCache, Tsit5Cache}, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) @tsit5pre0 @inbounds @simd ivdep for i in eachindex(out) out[i] = y₀[i] + @@ -439,7 +439,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Tsit5ConstantCache, Tsit5Cache}, idxs, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @tsit5pre0 @views @.. broadcast=false out=y₀[idxs] + dt * @@ -454,7 +454,7 @@ end @muladd function _ode_interpolant!(out::Array, Θ, dt, y₀, y₁, k, cache::Union{Tsit5ConstantCache, Tsit5Cache}, idxs, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @tsit5pre0 @inbounds for (j, i) in enumerate(idxs) out[j] = y₀[i] + @@ -476,7 +476,7 @@ end end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Tsit5ConstantCache, - idxs::Nothing, T::Type{Val{1}}) + idxs::Nothing, T::Type{Val{1}}, differential_vars::Nothing) @tsit5pre1 # return @.. broadcast=false k[1]*b1Θdiff + k[2]*b2Θdiff + k[3]*b3Θdiff + k[4]*b4Θdiff + k[5]*b5Θdiff + k[6]*b6Θdiff + k[7]*b7Θdiff return @inbounds k[1] * b1Θdiff + k[2] * b2Θdiff + k[3] * b3Θdiff + k[4] * b4Θdiff + @@ -484,7 +484,7 @@ end end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Tsit5Cache, idxs::Nothing, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars::Nothing) @tsit5pre1 return @inbounds @.. broadcast=false k[1]*b1Θdiff+k[2]*b2Θdiff+k[3]*b3Θdiff+ k[4]*b4Θdiff+k[5]*b5Θdiff+k[6]*b6Θdiff+k[7]*b7Θdiff @@ -492,7 +492,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Tsit5ConstantCache, Tsit5Cache}, idxs, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars::Nothing) @tsit5pre1 # return @.. broadcast=false k[1][idxs]*b1Θdiff + k[2][idxs]*b2Θdiff + k[3][idxs]*b3Θdiff + k[4][idxs]*b4Θdiff + k[5][idxs]*b5Θdiff + k[6][idxs]*b6Θdiff + k[7][idxs]*b7Θdiff return k[1][idxs] * b1Θdiff + k[2][idxs] * b2Θdiff + k[3][idxs] * b3Θdiff + @@ -502,7 +502,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Tsit5ConstantCache, Tsit5Cache}, - idxs::Nothing, T::Type{Val{1}}) + idxs::Nothing, T::Type{Val{1}}, differential_vars::Nothing) @tsit5pre1 @inbounds @.. broadcast=false out=k[1] * b1Θdiff + k[2] * b2Θdiff + k[3] * b3Θdiff + k[4] * b4Θdiff + k[5] * b5Θdiff + k[6] * b6Θdiff + @@ -515,7 +515,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Tsit5ConstantCache, Tsit5Cache}, idxs, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars::Nothing) @tsit5pre1 @views @.. broadcast=false out=k[1][idxs] * b1Θdiff + k[2][idxs] * b2Θdiff + k[3][idxs] * b3Θdiff + k[4][idxs] * b4Θdiff + @@ -541,7 +541,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Tsit5ConstantCache, Tsit5Cache}, - idxs::Nothing, T::Type{Val{2}}) + idxs::Nothing, T::Type{Val{2}}, differential_vars::Nothing) @tsit5pre2 # return @.. broadcast=false k[1]*b1Θdiff2 + k[2]*b2Θdiff2 + k[3]*b3Θdiff2 + k[4]*b4Θdiff2 + k[5]*b5Θdiff2 + k[6]*b6Θdiff2 + k[7]*b7Θdiff2 return @inbounds (k[1] * b1Θdiff2 + k[2] * b2Θdiff2 + k[3] * b3Θdiff2 + @@ -551,7 +551,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Tsit5ConstantCache, Tsit5Cache}, idxs, - T::Type{Val{2}}) + T::Type{Val{2}}, differential_vars::Nothing) @tsit5pre2 # return @.. broadcast=false k[1][idxs]*b1Θdiff2 + k[2][idxs]*b2Θdiff2 + k[3][idxs]*b3Θdiff2 + k[4][idxs]*b4Θdiff2 + k[5][idxs]*b5Θdiff2 + k[6][idxs]*b6Θdiff2 + k[7][idxs]*b7Θdiff2 return (k[1][idxs] * b1Θdiff2 + k[2][idxs] * b2Θdiff2 + k[3][idxs] * b3Θdiff2 + @@ -561,7 +561,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Tsit5ConstantCache, Tsit5Cache}, - idxs::Nothing, T::Type{Val{2}}) + idxs::Nothing, T::Type{Val{2}}, differential_vars::Nothing) @tsit5pre2 @inbounds @.. broadcast=false out=(k[1] * b1Θdiff2 + k[2] * b2Θdiff2 + k[3] * b3Θdiff2 + k[4] * b4Θdiff2 + k[5] * b5Θdiff2 + k[6] * b6Θdiff2 + @@ -574,7 +574,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Tsit5ConstantCache, Tsit5Cache}, idxs, - T::Type{Val{2}}) + T::Type{Val{2}}, differential_vars::Nothing) @tsit5pre2 @views @.. broadcast=false out=(k[1][idxs] * b1Θdiff2 + k[2][idxs] * b2Θdiff2 + k[3][idxs] * b3Θdiff2 + k[4][idxs] * b4Θdiff2 + @@ -600,7 +600,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Tsit5ConstantCache, Tsit5Cache}, - idxs::Nothing, T::Type{Val{3}}) + idxs::Nothing, T::Type{Val{3}}, differential_vars::Nothing) @tsit5pre3 # return @.. broadcast=false k[1]*b1Θdiff3 + k[2]*b2Θdiff3 + k[3]*b3Θdiff3 + k[4]*b4Θdiff3 + k[5]*b5Θdiff3 + k[6]*b6Θdiff3 + k[7]*b7Θdiff3 return @inbounds (k[1] * b1Θdiff3 + k[2] * b2Θdiff3 + k[3] * b3Θdiff3 + @@ -610,7 +610,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Tsit5ConstantCache, Tsit5Cache}, idxs, - T::Type{Val{3}}) + T::Type{Val{3}}, differential_vars::Nothing) @tsit5pre3 # return @.. broadcast=false k[1][idxs]*b1Θdiff3 + k[2][idxs]*b2Θdiff3 + k[3][idxs]*b3Θdiff3 + k[4][idxs]*b4Θdiff3 + k[5][idxs]*b5Θdiff3 + k[6][idxs]*b6Θdiff3 + k[7][idxs]*b7Θdiff3 return (k[1][idxs] * b1Θdiff3 + k[2][idxs] * b2Θdiff3 + k[3][idxs] * b3Θdiff3 + @@ -620,7 +620,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Tsit5ConstantCache, Tsit5Cache}, - idxs::Nothing, T::Type{Val{3}}) + idxs::Nothing, T::Type{Val{3}}, differential_vars::Nothing) @tsit5pre3 @inbounds @.. broadcast=false out=(k[1] * b1Θdiff3 + k[2] * b2Θdiff3 + k[3] * b3Θdiff3 + k[4] * b4Θdiff3 + k[5] * b5Θdiff3 + k[6] * b6Θdiff3 + @@ -633,7 +633,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Tsit5ConstantCache, Tsit5Cache}, idxs, - T::Type{Val{3}}) + T::Type{Val{3}}, differential_vars::Nothing) @tsit5pre3 @views @.. broadcast=false out=(k[1][idxs] * b1Θdiff3 + k[2][idxs] * b2Θdiff3 + k[3][idxs] * b3Θdiff3 + k[4][idxs] * b4Θdiff3 + @@ -659,7 +659,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Tsit5ConstantCache, Tsit5Cache}, - idxs::Nothing, T::Type{Val{4}}) + idxs::Nothing, T::Type{Val{4}}, differential_vars::Nothing) @tsit5pre4 # return @.. broadcast=false k[1]*b1Θdiff4 + k[2]*b2Θdiff4 + k[3]*b3Θdiff4 + k[4]*b4Θdiff4 + k[5]*b5Θdiff4 + k[6]*b6Θdiff4 + k[7]*b7Θdiff4 return @inbounds (k[1] * b1Θdiff4 + k[2] * b2Θdiff4 + k[3] * b3Θdiff4 + @@ -669,7 +669,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Tsit5ConstantCache, Tsit5Cache}, idxs, - T::Type{Val{4}}) + T::Type{Val{4}}, differential_vars::Nothing) @tsit5pre4 # return @.. broadcast=false k[1][idxs]*b1Θdiff4 + k[2][idxs]*b2Θdiff4 + k[3][idxs]*b3Θdiff4 + k[4][idxs]*b4Θdiff4 + k[5][idxs]*b5Θdiff4 + k[6][idxs]*b6Θdiff4 + k[7][idxs]*b7Θdiff4 return (k[1][idxs] * b1Θdiff4 + k[2][idxs] * b2Θdiff4 + k[3][idxs] * b3Θdiff4 + @@ -679,7 +679,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Tsit5ConstantCache, Tsit5Cache}, - idxs::Nothing, T::Type{Val{4}}) + idxs::Nothing, T::Type{Val{4}}, differential_vars::Nothing) @tsit5pre4 @inbounds @.. broadcast=false out=(k[1] * b1Θdiff4 + k[2] * b2Θdiff4 + k[3] * b3Θdiff4 + k[4] * b4Θdiff4 + k[5] * b5Θdiff4 + k[6] * b6Θdiff4 + @@ -692,7 +692,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Tsit5ConstantCache, Tsit5Cache}, idxs, - T::Type{Val{4}}) + T::Type{Val{4}}, differential_vars::Nothing) @tsit5pre4 @views @.. broadcast=false out=(k[1][idxs] * b1Θdiff4 + k[2][idxs] * b2Θdiff4 + k[3][idxs] * b3Θdiff4 + k[4][idxs] * b4Θdiff4 + @@ -725,7 +725,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen3ConstantCache, OwrenZen3Cache}, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) @owrenzen3pre0 @inbounds @.. broadcast=false y₀+dt * (k[1] * b1Θ + k[2] * b2Θ + k[3] * b3Θ + k[4] * b4Θ) @@ -733,7 +733,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen3ConstantCache, OwrenZen3Cache}, - idxs, T::Type{Val{0}}) + idxs, T::Type{Val{0}}, differential_vars::Nothing) @owrenzen3pre0 @views @.. broadcast=false y₀[idxs]+dt * (k[1][idxs] * b1Θ + k[2][idxs] * b2Θ + k[3][idxs] * b3Θ + @@ -742,7 +742,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen3ConstantCache, OwrenZen3Cache}, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) @owrenzen3pre0 @inbounds @.. broadcast=false out=y₀ + dt * @@ -752,7 +752,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen3ConstantCache, OwrenZen3Cache}, - idxs, T::Type{Val{0}}) + idxs, T::Type{Val{0}}, differential_vars::Nothing) @owrenzen3pre0 @views @.. broadcast=false out=y₀[idxs] + dt * @@ -771,14 +771,14 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen3ConstantCache, OwrenZen3Cache}, - idxs::Nothing, T::Type{Val{1}}) + idxs::Nothing, T::Type{Val{1}}, differential_vars::Nothing) @owrenzen3pre1 @inbounds @.. broadcast=false k[1]*b1Θdiff+k[2]*b2Θdiff+k[3]*b3Θdiff+k[4]*b4Θdiff end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen3ConstantCache, OwrenZen3Cache}, - idxs, T::Type{Val{1}}) + idxs, T::Type{Val{1}}, differential_vars::Nothing) @owrenzen3pre1 @views @.. broadcast=false k[1][idxs]*b1Θdiff+k[2][idxs]*b2Θdiff+k[3][idxs]*b3Θdiff+ k[4][idxs]*b4Θdiff @@ -786,7 +786,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen3ConstantCache, OwrenZen3Cache}, - idxs::Nothing, T::Type{Val{1}}) + idxs::Nothing, T::Type{Val{1}}, differential_vars::Nothing) @owrenzen3pre1 @inbounds @.. broadcast=false out=k[1] * b1Θdiff + k[2] * b2Θdiff + k[3] * b3Θdiff + k[4] * b4Θdiff @@ -795,7 +795,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen3ConstantCache, OwrenZen3Cache}, - idxs, T::Type{Val{1}}) + idxs, T::Type{Val{1}}, differential_vars::Nothing) @owrenzen3pre1 @views @.. broadcast=false out=k[1][idxs] * b1Θdiff + k[2][idxs] * b2Θdiff + k[3][idxs] * b3Θdiff + k[4][idxs] * b4Θdiff @@ -813,7 +813,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen3ConstantCache, OwrenZen3Cache}, - idxs::Nothing, T::Type{Val{2}}) + idxs::Nothing, T::Type{Val{2}}, differential_vars::Nothing) @owrenzen3pre2 @inbounds @.. broadcast=false (k[1] * b1Θdiff2 + k[2] * b2Θdiff2 + k[3] * b3Θdiff2 + k[4] * b4Θdiff2)*invdt @@ -821,7 +821,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen3ConstantCache, OwrenZen3Cache}, - idxs, T::Type{Val{2}}) + idxs, T::Type{Val{2}}, differential_vars::Nothing) @owrenzen3pre2 @views @.. broadcast=false (k[1][idxs] * b1Θdiff2 + k[2][idxs] * b2Θdiff2 + k[3][idxs] * b3Θdiff2 + k[4][idxs] * b4Θdiff2)*invdt @@ -829,7 +829,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen3ConstantCache, OwrenZen3Cache}, - idxs::Nothing, T::Type{Val{2}}) + idxs::Nothing, T::Type{Val{2}}, differential_vars::Nothing) @owrenzen3pre2 @inbounds @.. broadcast=false out=(k[1] * b1Θdiff2 + k[2] * b2Θdiff2 + k[3] * b3Θdiff2 + k[4] * b4Θdiff2) * invdt @@ -838,7 +838,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen3ConstantCache, OwrenZen3Cache}, - idxs, T::Type{Val{2}}) + idxs, T::Type{Val{2}}, differential_vars::Nothing) @owrenzen3pre2 @views @.. broadcast=false out=(k[1][idxs] * b1Θdiff2 + k[2][idxs] * b2Θdiff2 + k[3][idxs] * b3Θdiff2 + k[4][idxs] * b4Θdiff2) * invdt @@ -856,7 +856,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen3ConstantCache, OwrenZen3Cache}, - idxs::Nothing, T::Type{Val{3}}) + idxs::Nothing, T::Type{Val{3}}, differential_vars::Nothing) @owrenzen3pre3 @inbounds @.. broadcast=false (k[1] * b1Θdiff3 + k[2] * b2Θdiff3 + k[3] * b3Θdiff3 + k[4] * b4Θdiff3)*invdt2 @@ -864,7 +864,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen3ConstantCache, OwrenZen3Cache}, - idxs, T::Type{Val{3}}) + idxs, T::Type{Val{3}}, differential_vars::Nothing) @owrenzen3pre3 @views @.. broadcast=false (k[1][idxs] * b1Θdiff3 + k[2][idxs] * b2Θdiff3 + k[3][idxs] * b3Θdiff3 + k[4][idxs] * b4Θdiff3)*invdt2 @@ -872,7 +872,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen3ConstantCache, OwrenZen3Cache}, - idxs::Nothing, T::Type{Val{3}}) + idxs::Nothing, T::Type{Val{3}}, differential_vars::Nothing) @owrenzen3pre3 @inbounds @.. broadcast=false out=(k[1] * b1Θdiff3 + k[2] * b2Θdiff3 + k[3] * b3Θdiff3 + k[4] * b4Θdiff3) * invdt2 @@ -881,7 +881,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen3ConstantCache, OwrenZen3Cache}, - idxs, T::Type{Val{3}}) + idxs, T::Type{Val{3}}, differential_vars::Nothing) @owrenzen3pre3 @views @.. broadcast=false out=(k[1][idxs] * b1Θdiff3 + k[2][idxs] * b2Θdiff3 + k[3][idxs] * b3Θdiff3 + k[4][idxs] * b4Θdiff3) * invdt2 @@ -910,7 +910,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen4ConstantCache, OwrenZen4Cache}, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) @owrenzen4pre0 # return @.. broadcast=false y₀ + dt*(k[1]*b1Θ + k[3]*b3Θ + k[4]*b4Θ + k[5]*b5Θ + k[6]*b6Θ) return @inbounds y₀ + @@ -919,7 +919,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen4ConstantCache, OwrenZen4Cache}, - idxs, T::Type{Val{0}}) + idxs, T::Type{Val{0}}, differential_vars::Nothing) @owrenzen4pre0 # return @.. broadcast=false y₀[idxs] + dt*(k[1][idxs]*b1Θ + k[3][idxs]*b3Θ + # k[4][idxs]*b4Θ + k[5][idxs]*b5Θ + k[6][idxs]*b6Θ) @@ -930,7 +930,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen4ConstantCache, OwrenZen4Cache}, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) @owrenzen4pre0 @inbounds @.. broadcast=false out=y₀ + dt * @@ -945,7 +945,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen4ConstantCache, OwrenZen4Cache}, - idxs, T::Type{Val{0}}) + idxs, T::Type{Val{0}}, differential_vars::Nothing) @owrenzen4pre0 @inbounds @.. broadcast=false out=y₀[idxs] + dt * (k[1][idxs] * b1Θ + k[3][idxs] * b3Θ + @@ -969,7 +969,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen4ConstantCache, OwrenZen4Cache}, - idxs::Nothing, T::Type{Val{1}}) + idxs::Nothing, T::Type{Val{1}}, differential_vars::Nothing) @owrenzen4pre1 @inbounds @.. broadcast=false k[1]*b1Θdiff+k[3]*b3Θdiff+k[4]*b4Θdiff+k[5]*b5Θdiff+ k[6]*b6Θdiff @@ -977,7 +977,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen4ConstantCache, OwrenZen4Cache}, - idxs, T::Type{Val{1}}) + idxs, T::Type{Val{1}}, differential_vars::Nothing) @owrenzen4pre1 @views @.. broadcast=false k[1][idxs]*b1Θdiff+k[3][idxs]*b3Θdiff+k[4][idxs]*b4Θdiff+ k[5][idxs]*b5Θdiff+k[6][idxs]*b6Θdiff @@ -985,7 +985,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen4ConstantCache, OwrenZen4Cache}, - idxs::Nothing, T::Type{Val{1}}) + idxs::Nothing, T::Type{Val{1}}, differential_vars::Nothing) @owrenzen4pre1 @inbounds @.. broadcast=false out=k[1] * b1Θdiff + k[3] * b3Θdiff + k[4] * b4Θdiff + k[5] * b5Θdiff + k[6] * b6Θdiff @@ -994,7 +994,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen4ConstantCache, OwrenZen4Cache}, - idxs, T::Type{Val{1}}) + idxs, T::Type{Val{1}}, differential_vars::Nothing) @owrenzen4pre1 @views @.. broadcast=false out=k[1][idxs] * b1Θdiff + k[3][idxs] * b3Θdiff + k[4][idxs] * b4Θdiff + @@ -1014,7 +1014,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen4ConstantCache, OwrenZen4Cache}, - idxs::Nothing, T::Type{Val{2}}) + idxs::Nothing, T::Type{Val{2}}, differential_vars::Nothing) @owrenzen4pre2 @.. broadcast=false (k[1] * b1Θdiff2 + k[3] * b3Θdiff2 + k[4] * b4Θdiff2 + k[5] * b5Θdiff2 + k[6] * b6Θdiff2)*invdt @@ -1022,7 +1022,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen4ConstantCache, OwrenZen4Cache}, - idxs, T::Type{Val{2}}) + idxs, T::Type{Val{2}}, differential_vars::Nothing) @owrenzen4pre2 @views @.. broadcast=false (k[1][idxs] * b1Θdiff2 + k[3][idxs] * b3Θdiff2 + k[4][idxs] * b4Θdiff2 + @@ -1031,7 +1031,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen4ConstantCache, OwrenZen4Cache}, - idxs::Nothing, T::Type{Val{2}}) + idxs::Nothing, T::Type{Val{2}}, differential_vars::Nothing) @owrenzen4pre2 @inbounds @.. broadcast=false out=(k[1] * b1Θdiff2 + k[3] * b3Θdiff2 + k[4] * b4Θdiff2 + k[5] * b5Θdiff2 + k[6] * b6Θdiff2) * invdt @@ -1040,7 +1040,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen4ConstantCache, OwrenZen4Cache}, - idxs, T::Type{Val{2}}) + idxs, T::Type{Val{2}}, differential_vars::Nothing) @owrenzen4pre2 @views @.. broadcast=false out=(k[1][idxs] * b1Θdiff2 + k[3][idxs] * b3Θdiff2 + k[4][idxs] * b4Θdiff2 + @@ -1060,7 +1060,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen4ConstantCache, OwrenZen4Cache}, - idxs::Nothing, T::Type{Val{3}}) + idxs::Nothing, T::Type{Val{3}}, differential_vars::Nothing) @owrenzen4pre3 @inbounds @.. broadcast=false (k[1] * b1Θdiff3 + k[3] * b3Θdiff3 + k[4] * b4Θdiff3 + k[5] * b5Θdiff3 + k[6] * b6Θdiff3)*invdt2 @@ -1068,7 +1068,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen4ConstantCache, OwrenZen4Cache}, - idxs, T::Type{Val{3}}) + idxs, T::Type{Val{3}}, differential_vars::Nothing) @owrenzen4pre3 @views @.. broadcast=false (k[1][idxs] * b1Θdiff3 + k[3][idxs] * b3Θdiff3 + k[4][idxs] * b4Θdiff3 + @@ -1077,7 +1077,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen4ConstantCache, OwrenZen4Cache}, - idxs::Nothing, T::Type{Val{3}}) + idxs::Nothing, T::Type{Val{3}}, differential_vars::Nothing) @owrenzen4pre3 @inbounds @.. broadcast=false out=(k[1] * b1Θdiff3 + k[3] * b3Θdiff3 + k[4] * b4Θdiff3 + k[5] * b5Θdiff3 + k[6] * b6Θdiff3) * invdt2 @@ -1086,7 +1086,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen4ConstantCache, OwrenZen4Cache}, - idxs, T::Type{Val{3}}) + idxs, T::Type{Val{3}}, differential_vars::Nothing) @owrenzen4pre3 @views @.. broadcast=false out=(k[1][idxs] * b1Θdiff3 + k[3][idxs] * b3Θdiff3 + k[4][idxs] * b4Θdiff3 + @@ -1106,7 +1106,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen4ConstantCache, OwrenZen4Cache}, - idxs::Nothing, T::Type{Val{4}}) + idxs::Nothing, T::Type{Val{4}}, differential_vars::Nothing) @owrenzen4pre4 @.. broadcast=false (k[1] * b1Θdiff4 + k[3] * b3Θdiff4 + k[4] * b4Θdiff4 + k[5] * b5Θdiff4 + k[6] * b6Θdiff4)*invdt3 @@ -1114,7 +1114,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen4ConstantCache, OwrenZen4Cache}, - idxs, T::Type{Val{4}}) + idxs, T::Type{Val{4}}, differential_vars::Nothing) @owrenzen4pre4 @views @.. broadcast=false (k[1][idxs] * b1Θdiff4 + k[3][idxs] * b3Θdiff4 + k[4][idxs] * b4Θdiff4 + @@ -1123,7 +1123,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen4ConstantCache, OwrenZen4Cache}, - idxs::Nothing, T::Type{Val{4}}) + idxs::Nothing, T::Type{Val{4}}, differential_vars::Nothing) @owrenzen4pre4 @inbounds @.. broadcast=false out=(k[1] * b1Θdiff4 + k[3] * b3Θdiff4 + k[4] * b4Θdiff4 + k[5] * b5Θdiff4 + k[6] * b6Θdiff4) * invdt3 @@ -1132,7 +1132,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen4ConstantCache, OwrenZen4Cache}, - idxs, T::Type{Val{4}}) + idxs, T::Type{Val{4}}, differential_vars::Nothing) @owrenzen4pre4 @views @.. broadcast=false out=(k[1][idxs] * b1Θdiff4 + k[3][idxs] * b3Θdiff4 + k[4][idxs] * b4Θdiff4 + @@ -1164,7 +1164,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen5ConstantCache, OwrenZen5Cache}, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) @owrenzen5pre0 # return @.. broadcast=false y₀ + dt*(k[1]*b1Θ + k[3]*b3Θ + k[4]*b4Θ + k[5]*b5Θ + k[6]*b6Θ + # k[7]*b7Θ + k[8]*b8Θ) @@ -1175,7 +1175,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen5ConstantCache, OwrenZen5Cache}, - idxs, T::Type{Val{0}}) + idxs, T::Type{Val{0}}, differential_vars::Nothing) @owrenzen5pre0 # return @.. broadcast=false y₀[idxs] + dt*(k[1][idxs]*b1Θ + k[3][idxs]*b3Θ + # k[4][idxs]*b4Θ + k[5][idxs]*b5Θ + k[6][idxs]*b6Θ + @@ -1188,7 +1188,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen5ConstantCache, OwrenZen5Cache}, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) @owrenzen5pre0 @inbounds @.. broadcast=false out=y₀ + dt * @@ -1204,7 +1204,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen5ConstantCache, OwrenZen5Cache}, - idxs, T::Type{Val{0}}) + idxs, T::Type{Val{0}}, differential_vars::Nothing) @owrenzen5pre0 @views @.. broadcast=false out=y₀[idxs] + dt * (k[1][idxs] * b1Θ + k[3][idxs] * b3Θ + @@ -1230,7 +1230,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen5ConstantCache, OwrenZen5Cache}, - idxs::Nothing, T::Type{Val{1}}) + idxs::Nothing, T::Type{Val{1}}, differential_vars::Nothing) @owrenzen5pre1 return @inbounds k[1] * b1Θdiff + k[3] * b3Θdiff + k[4] * b4Θdiff + k[5] * b5Θdiff + k[6] * b6Θdiff + k[7] * b7Θdiff + k[8] * b8Θdiff @@ -1238,7 +1238,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen5ConstantCache, OwrenZen5Cache}, - idxs, T::Type{Val{1}}) + idxs, T::Type{Val{1}}, differential_vars::Nothing) @owrenzen5pre1 k[1][idxs] * b1Θdiff + k[3][idxs] * b3Θdiff + k[4][idxs] * b4Θdiff + k[5][idxs] * b5Θdiff + @@ -1247,7 +1247,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen5ConstantCache, OwrenZen5Cache}, - idxs::Nothing, T::Type{Val{1}}) + idxs::Nothing, T::Type{Val{1}}, differential_vars::Nothing) @owrenzen5pre1 @inbounds @.. broadcast=false out=k[1] * b1Θdiff + k[3] * b3Θdiff + k[4] * b4Θdiff + k[5] * b5Θdiff + k[6] * b6Θdiff + k[7] * b7Θdiff + @@ -1261,7 +1261,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen5ConstantCache, OwrenZen5Cache}, - idxs, T::Type{Val{1}}) + idxs, T::Type{Val{1}}, differential_vars::Nothing) @owrenzen5pre1 @views @.. broadcast=false out=k[1][idxs] * b1Θdiff + k[3][idxs] * b3Θdiff + k[4][idxs] * b4Θdiff + @@ -1288,7 +1288,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen5ConstantCache, OwrenZen5Cache}, - idxs::Nothing, T::Type{Val{2}}) + idxs::Nothing, T::Type{Val{2}}, differential_vars::Nothing) @owrenzen5pre2 return @inbounds (k[1] * b1Θdiff2 + k[3] * b3Θdiff2 + k[4] * b4Θdiff2 + k[5] * b5Θdiff2 + @@ -1297,7 +1297,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen5ConstantCache, OwrenZen5Cache}, - idxs, T::Type{Val{2}}) + idxs, T::Type{Val{2}}, differential_vars::Nothing) @owrenzen5pre2 (k[1][idxs] * b1Θdiff2 + k[3][idxs] * b3Θdiff2 + k[4][idxs] * b4Θdiff2 + k[5][idxs] * b5Θdiff2 + @@ -1306,7 +1306,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen5ConstantCache, OwrenZen5Cache}, - idxs::Nothing, T::Type{Val{2}}) + idxs::Nothing, T::Type{Val{2}}, differential_vars::Nothing) @owrenzen5pre2 @inbounds @.. broadcast=false out=(k[1] * b1Θdiff2 + k[3] * b3Θdiff2 + k[4] * b4Θdiff2 + k[5] * b5Θdiff2 + k[6] * b6Θdiff2 + k[7] * b7Θdiff2 + @@ -1320,7 +1320,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen5ConstantCache, OwrenZen5Cache}, - idxs, T::Type{Val{2}}) + idxs, T::Type{Val{2}}, differential_vars::Nothing) @owrenzen5pre2 @views @.. broadcast=false out=(k[1][idxs] * b1Θdiff2 + k[3][idxs] * b3Θdiff2 + k[4][idxs] * b4Θdiff2 + @@ -1347,7 +1347,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen5ConstantCache, OwrenZen5Cache}, - idxs::Nothing, T::Type{Val{3}}) + idxs::Nothing, T::Type{Val{3}}, differential_vars::Nothing) @owrenzen5pre3 return @inbounds (k[1] * b1Θdiff3 + k[3] * b3Θdiff3 + k[4] * b4Θdiff3 + k[5] * b5Θdiff3 + @@ -1356,7 +1356,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen5ConstantCache, OwrenZen5Cache}, - idxs, T::Type{Val{3}}) + idxs, T::Type{Val{3}}, differential_vars::Nothing) @owrenzen5pre3 (k[1][idxs] * b1Θdiff3 + k[3][idxs] * b3Θdiff3 + k[4][idxs] * b4Θdiff3 + k[5][idxs] * b5Θdiff3 + @@ -1365,7 +1365,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen5ConstantCache, OwrenZen5Cache}, - idxs::Nothing, T::Type{Val{3}}) + idxs::Nothing, T::Type{Val{3}}, differential_vars::Nothing) @owrenzen5pre3 @inbounds @.. broadcast=false out=(k[1] * b1Θdiff3 + k[3] * b3Θdiff3 + k[4] * b4Θdiff3 + k[5] * b5Θdiff3 + k[6] * b6Θdiff3 + k[7] * b7Θdiff3 + @@ -1379,7 +1379,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen5ConstantCache, OwrenZen5Cache}, - idxs, T::Type{Val{3}}) + idxs, T::Type{Val{3}}, differential_vars::Nothing) @owrenzen5pre3 @views @.. broadcast=false out=(k[1][idxs] * b1Θdiff3 + k[3][idxs] * b3Θdiff3 + k[4][idxs] * b4Θdiff3 + @@ -1406,7 +1406,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen5ConstantCache, OwrenZen5Cache}, - idxs::Nothing, T::Type{Val{4}}) + idxs::Nothing, T::Type{Val{4}}, differential_vars::Nothing) @owrenzen5pre4 return @inbounds (k[1] * b1Θdiff4 + k[3] * b3Θdiff4 + k[4] * b4Θdiff4 + k[5] * b5Θdiff4 + @@ -1415,7 +1415,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen5ConstantCache, OwrenZen5Cache}, - idxs, T::Type{Val{4}}) + idxs, T::Type{Val{4}}, differential_vars::Nothing) @owrenzen5pre4 (k[1][idxs] * b1Θdiff4 + k[3][idxs] * b3Θdiff4 + k[4][idxs] * b4Θdiff4 + k[5][idxs] * b5Θdiff4 + @@ -1424,7 +1424,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen5ConstantCache, OwrenZen5Cache}, - idxs::Nothing, T::Type{Val{4}}) + idxs::Nothing, T::Type{Val{4}}, differential_vars::Nothing) @owrenzen5pre4 @inbounds @.. broadcast=false out=(k[1] * b1Θdiff4 + k[3] * b3Θdiff4 + k[4] * b4Θdiff4 + k[5] * b5Θdiff4 + k[6] * b6Θdiff4 + k[7] * b7Θdiff4 + @@ -1438,7 +1438,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen5ConstantCache, OwrenZen5Cache}, - idxs, T::Type{Val{4}}) + idxs, T::Type{Val{4}}, differential_vars::Nothing) @owrenzen5pre4 @views @.. broadcast=false out=(k[1][idxs] * b1Θdiff4 + k[3][idxs] * b3Θdiff4 + k[4][idxs] * b4Θdiff4 + @@ -1465,7 +1465,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen5ConstantCache, OwrenZen5Cache}, - idxs::Nothing, T::Type{Val{5}}) + idxs::Nothing, T::Type{Val{5}}, differential_vars::Nothing) @owrenzen5pre5 return @inbounds (k[1] * b1Θdiff5 + k[3] * b3Θdiff5 + k[4] * b4Θdiff5 + k[5] * b5Θdiff5 + @@ -1474,7 +1474,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{OwrenZen5ConstantCache, OwrenZen5Cache}, - idxs, T::Type{Val{5}}) + idxs, T::Type{Val{5}}, differential_vars::Nothing) @owrenzen5pre5 (k[1][idxs] * b1Θdiff5 + k[3][idxs] * b3Θdiff5 + k[4][idxs] * b4Θdiff5 + k[5][idxs] * b5Θdiff5 + @@ -1483,7 +1483,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen5ConstantCache, OwrenZen5Cache}, - idxs::Nothing, T::Type{Val{5}}) + idxs::Nothing, T::Type{Val{5}}, differential_vars::Nothing) @owrenzen5pre5 @inbounds @.. broadcast=false out=(k[1] * b1Θdiff5 + k[3] * b3Θdiff5 + k[4] * b4Θdiff5 + k[5] * b5Θdiff5 + k[6] * b6Θdiff5 + k[7] * b7Θdiff5 + @@ -1497,7 +1497,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{OwrenZen5ConstantCache, OwrenZen5Cache}, - idxs, T::Type{Val{5}}) + idxs, T::Type{Val{5}}, differential_vars::Nothing) @owrenzen5pre5 @views @.. broadcast=false out=(k[1][idxs] * b1Θdiff5 + k[3][idxs] * b3Θdiff5 + k[4][idxs] * b4Θdiff5 + @@ -1538,7 +1538,7 @@ end end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::BS5ConstantCache, idxs::Nothing, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @bs5pre0 # return @.. broadcast=false y₀ + dt*Θ*k[1] + dt*(k[1]*b1Θ + k[3]*b3Θ + k[4]*b4Θ + k[5]*b5Θ + k[6]*b6Θ + k[7]*b7Θ + k[8]*b8Θ + k[9]*b9Θ + k[10]*b10Θ + k[11]*b11Θ) return @inbounds y₀ + dt * Θ * k[1] + @@ -1548,7 +1548,7 @@ end end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::BS5Cache, idxs::Nothing, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @bs5pre0 # return @.. broadcast=false y₀ + dt*Θ*k[1] + dt*(k[1]*b1Θ + k[3]*b3Θ + k[4]*b4Θ + k[5]*b5Θ + k[6]*b6Θ + k[7]*b7Θ + k[8]*b8Θ + k[9]*b9Θ + k[10]*b10Θ + k[11]*b11Θ) return @inbounds @.. broadcast=false y₀+dt*Θ*k[1]+ @@ -1560,7 +1560,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{BS5ConstantCache, BS5Cache}, idxs, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @bs5pre0 # return @.. broadcast=false y₀[idxs] + dt*Θ*k[1][idxs] + dt*(k[1][idxs]*b1Θ + k[3][idxs]*b3Θ + # k[4][idxs]*b4Θ + k[5][idxs]*b5Θ + k[6][idxs]*b6Θ + k[7][idxs]*b7Θ + @@ -1573,7 +1573,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{BS5ConstantCache, BS5Cache}, idxs::Nothing, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @bs5pre0 @inbounds @.. broadcast=false out=y₀ + dt * Θ * k[1] + dt * @@ -1588,7 +1588,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{BS5ConstantCache, BS5Cache}, idxs, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @bs5pre0 @views @.. broadcast=false out=y₀[idxs] + dt * Θ * k[1][idxs] + dt * @@ -1619,7 +1619,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{BS5ConstantCache, BS5Cache}, idxs::Nothing, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars::Nothing) @bs5pre1 # return @.. broadcast=false k[1] + k[1]*b1Θdiff + k[3]*b3Θdiff + k[4]*b4Θdiff + k[5]*b5Θdiff + k[6]*b6Θdiff + k[7]*b7Θdiff + k[8]*b8Θdiff + k[9]*b9Θdiff + k[10]*b10Θdiff + k[11]*b11Θdiff return @inbounds k[1] + k[1] * b1Θdiff + k[3] * b3Θdiff + k[4] * b4Θdiff + @@ -1630,7 +1630,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{BS5ConstantCache, BS5Cache}, idxs, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars::Nothing) @bs5pre1 # return @.. broadcast=false k[1][idxs] + k[1][idxs]*b1Θdiff + k[3][idxs]*b3Θdiff + # k[4][idxs]*b4Θdiff + k[5][idxs]*b5Θdiff + k[6][idxs]*b6Θdiff + @@ -1644,7 +1644,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{BS5ConstantCache, BS5Cache}, idxs::Nothing, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars::Nothing) @bs5pre1 @inbounds @.. broadcast=false out=k[1] + k[1] * b1Θdiff + k[3] * b3Θdiff + k[4] * b4Θdiff + k[5] * b5Θdiff + k[6] * b6Θdiff + @@ -1658,7 +1658,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{BS5ConstantCache, BS5Cache}, idxs, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars::Nothing) @bs5pre1 @views @.. broadcast=false out=k[1][idxs] + k[1][idxs] * b1Θdiff + k[3][idxs] * b3Θdiff + k[4][idxs] * b4Θdiff + @@ -1693,7 +1693,7 @@ end end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Vern6ConstantCache, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) @vern6pre0 #@.. broadcast=false y₀ + dt*(k[1]*b1Θ + k[4]*b4Θ + k[5]*b5Θ + k[6]*b6Θ + k[7]*b7Θ + k[8]*b8Θ + k[9]*b9Θ + k[10]*b10Θ + k[11]*b11Θ + k[12]*b12Θ) return @inbounds y₀ + @@ -1702,7 +1702,7 @@ end end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Vern6Cache, idxs::Nothing, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @vern6pre0 #@.. broadcast=false y₀ + dt*(k[1]*b1Θ + k[4]*b4Θ + k[5]*b5Θ + k[6]*b6Θ + k[7]*b7Θ + k[8]*b8Θ + k[9]*b9Θ + k[10]*b10Θ + k[11]*b11Θ + k[12]*b12Θ) return @inbounds @.. broadcast=false y₀+dt * (k[1] * b1Θ + k[4] * b4Θ + k[5] * b5Θ + @@ -1713,7 +1713,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Vern6ConstantCache, Vern6Cache}, idxs, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @vern6pre0 return y₀[idxs] + dt * (k[1][idxs] * b1Θ + k[4][idxs] * b4Θ + k[5][idxs] * b5Θ + @@ -1723,7 +1723,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Vern6ConstantCache, Vern6Cache}, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) @vern6pre0 @inbounds @.. broadcast=false out=y₀ + dt * @@ -1735,7 +1735,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Vern6ConstantCache, Vern6Cache}, idxs, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @vern6pre0 @views @.. broadcast=false out=y₀[idxs] + dt * @@ -1762,7 +1762,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Vern6ConstantCache, Vern6Cache}, - idxs::Nothing, T::Type{Val{1}}) + idxs::Nothing, T::Type{Val{1}}, differential_vars::Nothing) @vern6pre1 #@.. broadcast=false k[1]*b1Θdiff + k[4]*b4Θdiff + k[5]*b5Θdiff + k[6]*b6Θdiff + k[7]*b7Θdiff + k[8]*b8Θdiff + k[9]*b9Θdiff + k[10]*b10Θdiff + k[11]*b11Θdiff + k[12]*b12Θdiff return @inbounds k[1] * b1Θdiff + k[4] * b4Θdiff + k[5] * b5Θdiff + k[6] * b6Θdiff + @@ -1772,7 +1772,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Vern6ConstantCache, Vern6Cache}, idxs, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars::Nothing) @vern6pre1 return k[1][idxs] * b1Θdiff + k[4][idxs] * b4Θdiff + k[5][idxs] * b5Θdiff + k[6][idxs] * b6Θdiff + k[7][idxs] * b7Θdiff + k[8][idxs] * b8Θdiff + @@ -1782,7 +1782,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Vern6ConstantCache, Vern6Cache}, - idxs::Nothing, T::Type{Val{1}}) + idxs::Nothing, T::Type{Val{1}}, differential_vars::Nothing) @vern6pre1 @inbounds @.. broadcast=false out=k[1] * b1Θdiff + k[4] * b4Θdiff + k[5] * b5Θdiff + k[6] * b6Θdiff + k[7] * b7Θdiff + k[8] * b8Θdiff + @@ -1793,7 +1793,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Vern6ConstantCache, Vern6Cache}, idxs, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars::Nothing) @vern6pre1 @views @.. broadcast=false out=k[1][idxs] * b1Θdiff + k[4][idxs] * b4Θdiff + k[5][idxs] * b5Θdiff + k[6][idxs] * b6Θdiff + @@ -1834,7 +1834,7 @@ end end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Vern7ConstantCache, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) @vern7pre0 #@.. broadcast=false y₀ + dt*(k[1]*b1Θ + k[4]*b4Θ + k[5]*b5Θ + k[6]*b6Θ + k[7]*b7Θ + k[8]*b8Θ + k[9]*b9Θ + k[11]*b11Θ + k[12]*b12Θ + k[13]*b13Θ + k[14]*b14Θ + k[15]*b15Θ + k[16]*b16Θ) return @inbounds y₀ + @@ -1844,7 +1844,7 @@ end end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Vern7Cache, idxs::Nothing, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @vern7pre0 #@.. broadcast=false y₀ + dt*(k[1]*b1Θ + k[4]*b4Θ + k[5]*b5Θ + k[6]*b6Θ + k[7]*b7Θ + k[8]*b8Θ + k[9]*b9Θ + k[11]*b11Θ + k[12]*b12Θ + k[13]*b13Θ + k[14]*b14Θ + k[15]*b15Θ + k[16]*b16Θ) return @inbounds @.. broadcast=false y₀+dt * (k[1] * b1Θ + k[4] * b4Θ + k[5] * b5Θ + @@ -1856,7 +1856,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Vern7ConstantCache, Vern7Cache}, idxs, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @vern7pre0 return y₀[idxs] + dt * (k[1][idxs] * b1Θ + k[4][idxs] * b4Θ + k[5][idxs] * b5Θ + @@ -1867,7 +1867,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Vern7ConstantCache, Vern7Cache}, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) @vern7pre0 @inbounds @.. broadcast=false out=y₀ + dt * @@ -1880,7 +1880,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Vern7ConstantCache, Vern7Cache}, idxs, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @vern7pre0 @views @.. broadcast=false out=y₀[idxs] + dt * @@ -1912,7 +1912,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Vern7ConstantCache, Vern7Cache}, - idxs::Nothing, T::Type{Val{1}}) + idxs::Nothing, T::Type{Val{1}}, differential_vars::Nothing) @vern7pre1 #@.. broadcast=false k[1]*b1Θdiff + k[4]*b4Θdiff + k[5]*b5Θdiff + k[6]*b6Θdiff + k[7]*b7Θdiff + k[8]*b8Θdiff + k[9]*b9Θdiff + k[11]*b11Θdiff + k[12]*b12Θdiff + k[13]*b13Θdiff + k[14]*b14Θdiff + k[15]*b15Θdiff + k[16]*b16Θdiff return @inbounds k[1] * b1Θdiff + k[4] * b4Θdiff + k[5] * b5Θdiff + k[6] * b6Θdiff + @@ -1924,7 +1924,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Vern7ConstantCache, Vern7Cache}, idxs, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars::Nothing) @vern7pre1 return k[1][idxs] * b1Θdiff + k[4][idxs] * b4Θdiff + k[5][idxs] * b5Θdiff + k[6][idxs] * b6Θdiff + k[7][idxs] * b7Θdiff + k[8][idxs] * b8Θdiff + @@ -1935,7 +1935,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Vern7ConstantCache, Vern7Cache}, - idxs::Nothing, T::Type{Val{1}}) + idxs::Nothing, T::Type{Val{1}}, differential_vars::Nothing) @vern7pre1 @inbounds @.. broadcast=false out=k[1] * b1Θdiff + k[4] * b4Θdiff + k[5] * b5Θdiff + k[6] * b6Θdiff + k[7] * b7Θdiff + k[8] * b8Θdiff + @@ -1947,7 +1947,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Vern7ConstantCache, Vern7Cache}, idxs, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars::Nothing) @vern7pre1 @views @.. broadcast=false out=k[1][idxs] * b1Θdiff + k[4][idxs] * b4Θdiff + k[5][idxs] * b5Θdiff + k[6][idxs] * b6Θdiff + @@ -2001,7 +2001,7 @@ end end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Vern8ConstantCache, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) @vern8pre0 #@.. broadcast=false y₀ + dt*(k[1]*b1Θ + k[6]*b6Θ + k[7]*b7Θ + k[8]*b8Θ + k[9]*b9Θ + k[10]*b10Θ + k[11]*b11Θ + k[12]*b12Θ + k[14]*b14Θ + k[15]*b15Θ + k[16]*b16Θ + k[17]*b17Θ + k[18]*b18Θ + k[19]*b19Θ + k[20]*b20Θ + k[21]*b21Θ) return @inbounds y₀ + @@ -2014,7 +2014,7 @@ end end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Vern8Cache, idxs::Nothing, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @vern8pre0 #@.. broadcast=false y₀ + dt*(k[1]*b1Θ + k[6]*b6Θ + k[7]*b7Θ + k[8]*b8Θ + k[9]*b9Θ + k[10]*b10Θ + k[11]*b11Θ + k[12]*b12Θ + k[14]*b14Θ + k[15]*b15Θ + k[16]*b16Θ + k[17]*b17Θ + k[18]*b18Θ + k[19]*b19Θ + k[20]*b20Θ + k[21]*b21Θ) return @inbounds @.. broadcast=false y₀+dt * (k[1] * b1Θ + k[6] * b6Θ + k[7] * b7Θ + @@ -2028,7 +2028,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Vern8ConstantCache, Vern8Cache}, idxs, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @vern8pre0 return y₀[idxs] + dt * (k[1][idxs] * b1Θ + k[6][idxs] * b6Θ + k[7][idxs] * b7Θ + @@ -2041,7 +2041,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Vern8ConstantCache, Vern8Cache}, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) @vern8pre0 @inbounds @.. broadcast=false out=y₀ + dt * @@ -2055,7 +2055,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Vern8ConstantCache, Vern8Cache}, idxs, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @vern8pre0 @views @.. broadcast=false out=y₀[idxs] + dt * @@ -2107,7 +2107,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Vern8ConstantCache, Vern8Cache}, - idxs::Nothing, T::Type{Val{1}}) + idxs::Nothing, T::Type{Val{1}}, differential_vars::Nothing) @vern8pre1 #@.. broadcast=false k[1]*b1Θdiff + k[6]*b6Θdiff + k[7]*b7Θdiff + k[8]*b8Θdiff + k[9]*b9Θdiff + k[10]*b10Θdiff + k[11]*b11Θdiff + k[12]*b12Θdiff + k[14]*b14Θdiff + k[15]*b15Θdiff + k[16]*b16Θdiff + k[17]*b17Θdiff + k[18]*b18Θdiff + k[19]*b19Θdiff + k[20]*b20Θdiff + k[21]*b21Θdiff return @inbounds k[1] * b1Θdiff + k[6] * b6Θdiff + k[7] * b7Θdiff + k[8] * b8Θdiff + @@ -2121,7 +2121,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Vern8ConstantCache, Vern8Cache}, idxs, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars::Nothing) @vern8pre1 return k[1][idxs] * b1Θdiff + k[6][idxs] * b6Θdiff + k[7][idxs] * b7Θdiff + k[8][idxs] * b8Θdiff + k[9][idxs] * b9Θdiff + k[10][idxs] * b10Θdiff + @@ -2133,7 +2133,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Vern8ConstantCache, Vern8Cache}, - idxs::Nothing, T::Type{Val{1}}) + idxs::Nothing, T::Type{Val{1}}, differential_vars::Nothing) @vern8pre1 @inbounds @.. broadcast=false out=k[1] * b1Θdiff + k[6] * b6Θdiff + k[7] * b7Θdiff + k[8] * b8Θdiff + k[9] * b9Θdiff + k[10] * b10Θdiff + @@ -2147,7 +2147,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Vern8ConstantCache, Vern8Cache}, idxs, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars::Nothing) @vern8pre1 @views @.. broadcast=false out=k[1][idxs] * b1Θdiff + k[6][idxs] * b6Θdiff + k[7][idxs] * b7Θdiff + k[8][idxs] * b8Θdiff + @@ -2210,7 +2210,7 @@ end end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Vern9ConstantCache, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) @vern9pre0 #@.. broadcast=false y₀ + dt*(k[1]*b1Θ + k[2]*b8Θ + k[3]*b9Θ + k[4]*b10Θ + k[5]*b11Θ + k[6]*b12Θ + k[7]*b13Θ + k[8]*b14Θ + k[9]*b15Θ + k[11]*b17Θ + k[12]*b18Θ + k[13]*b19Θ + k[14]*b20Θ + k[15]*b21Θ + k[16]*b22Θ + k[17]*b23Θ + k[18]*b24Θ + k[19]*b25Θ + k[20]*b26Θ) return @inbounds y₀ + @@ -2223,7 +2223,7 @@ end end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Vern9Cache, idxs::Nothing, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @vern9pre0 #@.. broadcast=false y₀ + dt*(k[1]*b1Θ + k[2]*b8Θ + k[3]*b9Θ + k[4]*b10Θ + k[5]*b11Θ + k[6]*b12Θ + k[7]*b13Θ + k[8]*b14Θ + k[9]*b15Θ + k[11]*b17Θ + k[12]*b18Θ + k[13]*b19Θ + k[14]*b20Θ + k[15]*b21Θ + k[16]*b22Θ + k[17]*b23Θ + k[18]*b24Θ + k[19]*b25Θ + k[20]*b26Θ) return @inbounds @.. broadcast=false y₀+dt * (k[1] * b1Θ + k[2] * b8Θ + k[3] * b9Θ + @@ -2238,7 +2238,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Vern9ConstantCache, Vern9Cache}, idxs, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @vern9pre0 return y₀[idxs] + dt * (k[1][idxs] * b1Θ + k[2][idxs] * b8Θ + k[3][idxs] * b9Θ + @@ -2252,7 +2252,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Vern9ConstantCache, Vern9Cache}, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) @vern9pre0 @inbounds @.. broadcast=false out=y₀ + dt * @@ -2267,7 +2267,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Vern9ConstantCache, Vern9Cache}, idxs, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @vern9pre0 @views @.. broadcast=false out=y₀[idxs] + dt * @@ -2345,7 +2345,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Vern9ConstantCache, Vern9Cache}, - idxs::Nothing, T::Type{Val{1}}) + idxs::Nothing, T::Type{Val{1}}, differential_vars::Nothing) @vern9pre1 #@.. broadcast=false k[1]*b1Θdiff + k[2]*b8Θdiff + k[3]*b9Θdiff + k[4]*b10Θdiff + k[5]*b11Θdiff + k[6]*b12Θdiff + k[7]*b13Θdiff + k[8]*b14Θdiff + k[9]*b15Θdiff + k[11]*b17Θdiff + k[12]*b18Θdiff + k[13]*b19Θdiff + k[14]*b20Θdiff + k[15]*b21Θdiff + k[16]*b22Θdiff + k[17]*b23Θdiff + k[18]*b24Θdiff + k[19]*b25Θdiff + k[20]*b26Θdiff return @inbounds k[1] * b1Θdiff + k[2] * b8Θdiff + k[3] * b9Θdiff + k[4] * b10Θdiff + @@ -2359,7 +2359,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Vern9ConstantCache, Vern9Cache}, idxs, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars::Nothing) @vern9pre1 return k[1][idxs] * b1Θdiff + k[2][idxs] * b8Θdiff + k[3][idxs] * b9Θdiff + k[4][idxs] * b10Θdiff + k[5][idxs] * b11Θdiff + @@ -2374,7 +2374,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Vern9ConstantCache, Vern9Cache}, - idxs::Nothing, T::Type{Val{1}}) + idxs::Nothing, T::Type{Val{1}}, differential_vars::Nothing) @vern9pre1 @inbounds @.. broadcast=false out=k[1] * b1Θdiff + k[2] * b8Θdiff + k[3] * b9Θdiff + k[4] * b10Θdiff + k[5] * b11Θdiff + k[6] * b12Θdiff + @@ -2389,7 +2389,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Vern9ConstantCache, Vern9Cache}, idxs, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars::Nothing) @vern9pre1 @views @.. broadcast=false out=k[1][idxs] * b1Θdiff + k[2][idxs] * b8Θdiff + k[3][idxs] * b9Θdiff + k[4][idxs] * b10Θdiff + @@ -2449,7 +2449,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Vern9ConstantCache, Vern9Cache}, - idxs::Nothing, T::Type{Val{2}}) + idxs::Nothing, T::Type{Val{2}}, differential_vars::Nothing) @vern9pre2 #@.. broadcast=false k[1]*b1Θdiff + k[2]*b8Θdiff + k[3]*b9Θdiff + k[4]*b10Θdiff + k[5]*b11Θdiff + k[6]*b12Θdiff + k[7]*b13Θdiff + k[8]*b14Θdiff + k[9]*b15Θdiff + k[11]*b17Θdiff + k[12]*b18Θdiff + k[13]*b19Θdiff + k[14]*b20Θdiff + k[15]*b21Θdiff + k[16]*b22Θdiff + k[17]*b23Θdiff + k[18]*b24Θdiff + k[19]*b25Θdiff + k[20]*b26Θdiff return @inbounds (k[1] * b1Θdiff + k[2] * b8Θdiff + k[3] * b9Θdiff + k[4] * b10Θdiff + @@ -2464,7 +2464,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Vern9ConstantCache, Vern9Cache}, idxs, - T::Type{Val{2}}) + T::Type{Val{2}}, differential_vars::Nothing) @vern9pre2 return (k[1][idxs] * b1Θdiff + k[2][idxs] * b8Θdiff + k[3][idxs] * b9Θdiff + k[4][idxs] * b10Θdiff + k[5][idxs] * b11Θdiff + @@ -2479,7 +2479,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Vern9ConstantCache, Vern9Cache}, - idxs::Nothing, T::Type{Val{2}}) + idxs::Nothing, T::Type{Val{2}}, differential_vars::Nothing) @vern9pre2 @inbounds @.. broadcast=false out=(k[1] * b1Θdiff + k[2] * b8Θdiff + k[3] * b9Θdiff + k[4] * b10Θdiff + k[5] * b11Θdiff + k[6] * b12Θdiff + @@ -2494,7 +2494,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Vern9ConstantCache, Vern9Cache}, idxs, - T::Type{Val{2}}) + T::Type{Val{2}}, differential_vars::Nothing) @vern9pre2 @views @.. broadcast=false out=(k[1][idxs] * b1Θdiff + k[2][idxs] * b8Θdiff + k[3][idxs] * b9Θdiff + k[4][idxs] * b10Θdiff + @@ -2551,7 +2551,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Vern9ConstantCache, Vern9Cache}, - idxs::Nothing, T::Type{Val{3}}) + idxs::Nothing, T::Type{Val{3}}, differential_vars::Nothing) @vern9pre3 #@.. broadcast=false k[1]*b1Θdiff + k[2]*b8Θdiff + k[3]*b9Θdiff + k[4]*b10Θdiff + k[5]*b11Θdiff + k[6]*b12Θdiff + k[7]*b13Θdiff + k[8]*b14Θdiff + k[9]*b15Θdiff + k[11]*b17Θdiff + k[12]*b18Θdiff + k[13]*b19Θdiff + k[14]*b20Θdiff + k[15]*b21Θdiff + k[16]*b22Θdiff + k[17]*b23Θdiff + k[18]*b24Θdiff + k[19]*b25Θdiff + k[20]*b26Θdiff return @inbounds (k[1] * b1Θdiff + k[2] * b8Θdiff + k[3] * b9Θdiff + k[4] * b10Θdiff + @@ -2566,7 +2566,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Vern9ConstantCache, Vern9Cache}, idxs, - T::Type{Val{3}}) + T::Type{Val{3}}, differential_vars::Nothing) @vern9pre3 return (k[1][idxs] * b1Θdiff + k[2][idxs] * b8Θdiff + k[3][idxs] * b9Θdiff + k[4][idxs] * b10Θdiff + k[5][idxs] * b11Θdiff + @@ -2581,7 +2581,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Vern9ConstantCache, Vern9Cache}, - idxs::Nothing, T::Type{Val{3}}) + idxs::Nothing, T::Type{Val{3}}, differential_vars::Nothing) @vern9pre3 @inbounds @.. broadcast=false out=(k[1] * b1Θdiff + k[2] * b8Θdiff + k[3] * b9Θdiff + k[4] * b10Θdiff + k[5] * b11Θdiff + k[6] * b12Θdiff + @@ -2596,7 +2596,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Vern9ConstantCache, Vern9Cache}, idxs, - T::Type{Val{3}}) + T::Type{Val{3}}, differential_vars::Nothing) @vern9pre3 @views @.. broadcast=false out=(k[1][idxs] * b1Θdiff + k[2][idxs] * b8Θdiff + k[3][idxs] * b9Θdiff + k[4][idxs] * b10Θdiff + @@ -2637,7 +2637,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Vern9ConstantCache, Vern9Cache}, - idxs::Nothing, T::Type{Val{4}}) + idxs::Nothing, T::Type{Val{4}}, differential_vars::Nothing) @vern9pre4 #@.. broadcast=false k[1]*b1Θdiff + k[2]*b8Θdiff + k[3]*b9Θdiff + k[4]*b10Θdiff + k[5]*b11Θdiff + k[6]*b12Θdiff + k[7]*b13Θdiff + k[8]*b14Θdiff + k[9]*b15Θdiff + k[11]*b17Θdiff + k[12]*b18Θdiff + k[13]*b19Θdiff + k[14]*b20Θdiff + k[15]*b21Θdiff + k[16]*b22Θdiff + k[17]*b23Θdiff + k[18]*b24Θdiff + k[19]*b25Θdiff + k[20]*b26Θdiff return @inbounds (k[1] * b1Θdiff + k[2] * b8Θdiff + k[3] * b9Θdiff + k[4] * b10Θdiff + @@ -2652,7 +2652,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Vern9ConstantCache, Vern9Cache}, idxs, - T::Type{Val{4}}) + T::Type{Val{4}}, differential_vars::Nothing) @vern9pre4 return (k[1][idxs] * b1Θdiff + k[2][idxs] * b8Θdiff + k[3][idxs] * b9Θdiff + k[4][idxs] * b10Θdiff + k[5][idxs] * b11Θdiff + @@ -2667,7 +2667,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Vern9ConstantCache, Vern9Cache}, - idxs::Nothing, T::Type{Val{4}}) + idxs::Nothing, T::Type{Val{4}}, differential_vars::Nothing) @vern9pre4 @inbounds @.. broadcast=false out=(k[1] * b1Θdiff + k[2] * b8Θdiff + k[3] * b9Θdiff + k[4] * b10Θdiff + k[5] * b11Θdiff + k[6] * b12Θdiff + @@ -2682,7 +2682,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Vern9ConstantCache, Vern9Cache}, idxs, - T::Type{Val{4}}) + T::Type{Val{4}}, differential_vars::Nothing) @vern9pre4 @views @.. broadcast=false out=(k[1][idxs] * b1Θdiff + k[2][idxs] * b8Θdiff + k[3][idxs] * b9Θdiff + k[4][idxs] * b10Θdiff + @@ -2701,7 +2701,7 @@ end """ @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{DP8ConstantCache, DP8Cache}, idxs::Nothing, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) Θ1 = 1 - Θ # return @.. broadcast=false y₀ + dt*Θ*(k[1] + Θ1*(k[2] + Θ*(k[3]+Θ1*(k[4] + Θ*(k[5] + Θ1*(k[6]+Θ*k[7])))))) return @inbounds y₀ + @@ -2713,7 +2713,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{DP8ConstantCache, DP8Cache}, idxs, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) Θ1 = 1 - Θ # return @.. broadcast=false y₀[idxs] + dt*Θ*(k[1][idxs] + Θ1*(k[2][idxs] + Θ*(k[3][idxs]+Θ1*(k[4][idxs] + Θ*(k[5][idxs] + Θ1*(k[6][idxs]+Θ*k[7][idxs])))))) return y₀[idxs] + @@ -2726,7 +2726,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{DP8ConstantCache, DP8Cache}, idxs::Nothing, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) Θ1 = 1 - Θ @inbounds @.. broadcast=false out=y₀ + dt * Θ * @@ -2743,7 +2743,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{DP8ConstantCache, DP8Cache}, idxs, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) Θ1 = 1 - Θ @views @.. broadcast=false out=y₀[idxs] + dt * Θ * @@ -2761,7 +2761,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{DP8ConstantCache, DP8Cache}, idxs::Nothing, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars::Nothing) @inbounds b1diff = @.. broadcast=false k[1]+k[2] @inbounds b2diff = @.. broadcast=false -2*k[2]+2*k[3]+2*k[4] @inbounds b3diff = @.. broadcast=false -3 * k[3]-6 * k[4]+3*k[5]+3*k[6] @@ -2777,7 +2777,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{DP8ConstantCache, DP8Cache}, idxs, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars::Nothing) b1diff = @.. broadcast=false k[1][idxs]+k[2][idxs] b2diff = @.. broadcast=false -2*k[2][idxs]+2*k[3][idxs]+2*k[4][idxs] b3diff = @.. broadcast=false -3 * k[3][idxs]-6 * k[4][idxs]+3*k[5][idxs]+3*k[6][idxs] @@ -2794,7 +2794,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{DP8ConstantCache, DP8Cache}, idxs::Nothing, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars::Nothing) # b1diff = k[1] + k[2] # b2diff = -2*k[2] + 2*k[3] + 2*k[4] # b3diff = -3*k[3] - 6*k[4] + 3*k[5] + 3*k[6] @@ -2814,7 +2814,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{DP8ConstantCache, DP8Cache}, idxs, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars::Nothing) # b1diff = k[1][idxs] + k[2][idxs] # b2diff = -2*k[2][idxs] + 2*k[3][idxs] + 2*k[4][idxs] # b3diff = -3*k[3][idxs] - 6*k[4][idxs] + 3*k[5][idxs] + 3*k[6][idxs] @@ -2873,7 +2873,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) @dprkn6pre0 return ArrayPartition(duprev + dt * Θ * @@ -2888,7 +2888,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, idxs, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @dprkn6pre0 return ArrayPartition(duprev[idxs] + dt * Θ * @@ -2905,7 +2905,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) @dprkn6pre0 @inbounds @.. broadcast=false out.x[2]=uprev + dt * Θ * @@ -2930,7 +2930,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, idxs, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars::Nothing) @dprkn6pre0 @inbounds @.. broadcast=false out.x[2]=uprev[idxs] + dt * Θ * diff --git a/src/dense/rosenbrock_interpolants.jl b/src/dense/rosenbrock_interpolants.jl index 7247938945..06b5950710 100644 --- a/src/dense/rosenbrock_interpolants.jl +++ b/src/dense/rosenbrock_interpolants.jl @@ -18,14 +18,14 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Rosenbrock23ConstantCache, Rosenbrock32ConstantCache}, idxs::Nothing, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars) @rosenbrock2332pre0 @inbounds y₀ + dt * (c1 * k[1] + c2 * k[2]) end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Rosenbrock23Cache, Rosenbrock32Cache}, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars) @rosenbrock2332pre0 @inbounds @.. broadcast=false y₀+dt * (c1 * k[1] + c2 * k[2]) end @@ -33,7 +33,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Rosenbrock23ConstantCache, Rosenbrock23Cache, Rosenbrock32ConstantCache, Rosenbrock32Cache, - }, idxs, T::Type{Val{0}}) + }, idxs, T::Type{Val{0}}, differential_vars) @rosenbrock2332pre0 @.. broadcast=false y₀[idxs]+dt * (c1 * k[1][idxs] + c2 * k[2][idxs]) end @@ -42,14 +42,14 @@ end cache::Union{Rosenbrock23ConstantCache, Rosenbrock23Cache, Rosenbrock32ConstantCache, Rosenbrock32Cache, - }, idxs::Nothing, T::Type{Val{0}}) + }, idxs::Nothing, T::Type{Val{0}}, differential_vars) @rosenbrock2332pre0 @inbounds @.. broadcast=false out=y₀ + dt * (c1 * k[1] + c2 * k[2]) out end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Rosenbrock23Cache{<:Array}, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars) @rosenbrock2332pre0 @inbounds @simd ivdep for i in eachindex(out) out[i] = y₀[i] + dt * (c1 * k[1][i] + c2 * k[2][i]) @@ -61,7 +61,7 @@ end cache::Union{Rosenbrock23ConstantCache, Rosenbrock23Cache, Rosenbrock32ConstantCache, Rosenbrock32Cache, - }, idxs, T::Type{Val{0}}) + }, idxs, T::Type{Val{0}}, differential_vars) @rosenbrock2332pre0 @views @.. broadcast=false out=y₀[idxs] + dt * (c1 * k[1][idxs] + c2 * k[2][idxs]) out @@ -77,7 +77,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Rosenbrock23ConstantCache, Rosenbrock23Cache, Rosenbrock32ConstantCache, Rosenbrock32Cache, - }, idxs::Nothing, T::Type{Val{1}}) + }, idxs::Nothing, T::Type{Val{1}}, differential_vars) @rosenbrock2332pre1 @.. broadcast=false c1diff * k[1]+c2diff * k[2] end @@ -85,7 +85,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Rosenbrock23ConstantCache, Rosenbrock23Cache, Rosenbrock32ConstantCache, Rosenbrock32Cache, - }, idxs, T::Type{Val{1}}) + }, idxs, T::Type{Val{1}}, differential_vars) @rosenbrock2332pre1 @.. broadcast=false c1diff * k[1][idxs]+c2diff * k[2][idxs] end @@ -94,7 +94,7 @@ end cache::Union{Rosenbrock23ConstantCache, Rosenbrock23Cache, Rosenbrock32ConstantCache, Rosenbrock32Cache, - }, idxs::Nothing, T::Type{Val{1}}) + }, idxs::Nothing, T::Type{Val{1}}, differential_vars) @rosenbrock2332pre1 @.. broadcast=false out=c1diff * k[1] + c2diff * k[2] out @@ -104,7 +104,7 @@ end cache::Union{Rosenbrock23ConstantCache, Rosenbrock23Cache, Rosenbrock32ConstantCache, Rosenbrock32Cache, - }, idxs, T::Type{Val{1}}) + }, idxs, T::Type{Val{1}}, differential_vars) @rosenbrock2332pre1 @views @.. broadcast=false out=c1diff * k[1][idxs] + c2diff * k[2][idxs] out @@ -114,34 +114,34 @@ end From MATLAB ODE Suite by Shampine """ @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Rodas4ConstantCache, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars) Θ1 = 1 - Θ @inbounds Θ1 * y₀ + Θ * (y₁ + Θ1 * (k[1] + Θ * k[2])) end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Rodas4Cache, idxs::Nothing, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars) Θ1 = 1 - Θ @inbounds @.. broadcast=false Θ1 * y₀+Θ * (y₁ + Θ1 * (k[1] + Θ * k[2])) end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Rodas4ConstantCache, Rodas4Cache}, idxs, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars) Θ1 = 1 - Θ @.. broadcast=false Θ1 * y₀[idxs]+Θ * (y₁[idxs] + Θ1 * (k[1][idxs] + Θ * k[2][idxs])) end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Rodas4ConstantCache, Rodas4Cache}, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars) Θ1 = 1 - Θ @.. broadcast=false out=Θ1 * y₀ + Θ * (y₁ + Θ1 * (k[1] + Θ * k[2])) out end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Rodas4Cache{<:Array}, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars) Θ1 = 1 - Θ @inbounds @simd ivdep for i in eachindex(out) out[i] = Θ1 * y₀[i] + Θ * (y₁[i] + Θ1 * (k[1][i] + Θ * k[2][i])) @@ -151,7 +151,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Rodas4ConstantCache, Rodas4Cache}, idxs, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars) Θ1 = 1 - Θ @views @.. broadcast=false out=Θ1 * y₀[idxs] + Θ * (y₁[idxs] + Θ1 * (k[1][idxs] + Θ * k[2][idxs])) @@ -160,19 +160,19 @@ end # First Derivative @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Rodas4ConstantCache, - idxs::Nothing, T::Type{Val{1}}) + idxs::Nothing, T::Type{Val{1}}, differential_vars) @inbounds (k[1] + Θ * (-2 * k[1] + 2 * k[2] - 3 * k[2] * Θ) - y₀ + y₁) / dt end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Rodas4Cache, idxs::Nothing, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars) @inbounds @.. broadcast=false (k[1] + Θ * (-2 * k[1] + 2 * k[2] - 3 * k[2] * Θ) - y₀ + y₁)/dt end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Rodas4ConstantCache, Rodas4Cache}, idxs, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars) @.. broadcast=false (k[1][idxs] + Θ * (-2 * k[1][idxs] + 2 * k[2][idxs] - 3 * k[2][idxs] * Θ) - y₀[idxs] + y₁[idxs])/dt @@ -180,7 +180,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Rodas4ConstantCache, Rodas4Cache}, - idxs::Nothing, T::Type{Val{1}}) + idxs::Nothing, T::Type{Val{1}}, differential_vars) @.. broadcast=false out=(k[1] + Θ * (-2 * k[1] + 2 * k[2] - 3 * k[2] * Θ) - y₀ + y₁) / dt out @@ -188,7 +188,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Rodas4ConstantCache, Rodas4Cache}, idxs, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars) @views @.. broadcast=false out=(k[1][idxs] + Θ * (-2 * k[1][idxs] + 2 * k[2][idxs] - @@ -200,20 +200,20 @@ end #- @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Rosenbrock5ConstantCache, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars) Θ1 = 1 - Θ @inbounds Θ1 * y₀ + Θ * (y₁ + Θ1 * (k[1] + Θ * (k[2] + Θ * k[3]))) end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Rosenbrock5Cache, idxs::Nothing, - T::Type{Val{0}}) + T::Type{Val{0}}, differential_vars) Θ1 = 1 - Θ @inbounds @.. broadcast=false Θ1 * y₀+Θ * (y₁ + Θ1 * (k[1] + Θ * (k[2] + Θ * k[3]))) end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Rosenbrock5ConstantCache, Rosenbrock5Cache}, - idxs, T::Type{Val{0}}) + idxs, T::Type{Val{0}}, differential_vars) Θ1 = 1 - Θ @.. broadcast=false Θ1 * y₀[idxs]+Θ * (y₁[idxs] + @@ -222,14 +222,14 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Rosenbrock5ConstantCache, Rosenbrock5Cache}, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars) Θ1 = 1 - Θ @.. broadcast=false out=Θ1 * y₀ + Θ * (y₁ + Θ1 * (k[1] + Θ * (k[2] + Θ * k[3]))) out end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Rosenbrock5Cache{<:Array}, - idxs::Nothing, T::Type{Val{0}}) + idxs::Nothing, T::Type{Val{0}}, differential_vars) Θ1 = 1 - Θ @inbounds @simd ivdep for i in eachindex(out) out[i] = Θ1 * y₀[i] + Θ * (y₁[i] + Θ1 * (k[1][i] + Θ * (k[2][i] + Θ * k[3][i]))) @@ -239,7 +239,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Rosenbrock5ConstantCache, Rosenbrock5Cache}, - idxs, T::Type{Val{0}}) + idxs, T::Type{Val{0}}, differential_vars) Θ1 = 1 - Θ @views @.. broadcast=false out=Θ1 * y₀[idxs] + Θ * (y₁[idxs] + @@ -249,14 +249,14 @@ end # First Derivative @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Rosenbrock5ConstantCache, - idxs::Nothing, T::Type{Val{1}}) + idxs::Nothing, T::Type{Val{1}}, differential_vars) @inbounds (k[1] + Θ * (-2 * k[1] + 2 * k[2] + Θ * (-3 * k[2] + 3 * k[3] - 4 * Θ * k[3])) - y₀ + y₁) / dt end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Rosenbrock5Cache, idxs::Nothing, - T::Type{Val{1}}) + T::Type{Val{1}}, differential_vars) @inbounds @.. broadcast=false (k[1] + Θ * (-2 * k[1] + 2 * k[2] + Θ * (-3 * k[2] + 3 * k[3] - 4 * Θ * k[3])) - y₀ + y₁)/dt @@ -264,7 +264,7 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{Rosenbrock5ConstantCache, Rosenbrock5Cache}, - idxs, T::Type{Val{1}}) + idxs, T::Type{Val{1}}, differential_vars) @.. broadcast=false (k[1][idxs] + Θ * (-2 * k[1][idxs] + 2 * k[2][idxs] + Θ * (-3 * k[2][idxs] + 3 * k[3][idxs] - 4 * Θ * k[3][idxs])) - @@ -273,7 +273,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Rosenbrock5ConstantCache, Rosenbrock5Cache}, - idxs::Nothing, T::Type{Val{1}}) + idxs::Nothing, T::Type{Val{1}}, differential_vars) @.. broadcast=false out=(k[1] + Θ * (-2 * k[1] + 2 * k[2] + Θ * (-3 * k[2] + 3 * k[3] - 4 * Θ * k[3])) - y₀ + y₁) / dt @@ -282,7 +282,7 @@ end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{Rosenbrock5ConstantCache, Rosenbrock5Cache}, - idxs, T::Type{Val{1}}) + idxs, T::Type{Val{1}}, differential_vars) @views @.. broadcast=false out=(k[1][idxs] + Θ * (-2 * k[1][idxs] + 2 * k[2][idxs] + Θ * diff --git a/src/integrators/type.jl b/src/integrators/type.jl index 6f96de274a..cdf3b077c1 100644 --- a/src/integrators/type.jl +++ b/src/integrators/type.jl @@ -85,7 +85,7 @@ For more info see the linked documentation page. mutable struct ODEIntegrator{algType <: Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}, IIP, uType, duType, tType, pType, eigenType, EEstT, QT, tdirType, ksEltype, SolType, F, CacheType, O, FSALType, EventErrorType, - CallbackCacheType, IA} <: + CallbackCacheType, IA, DV} <: DiffEqBase.AbstractODEIntegrator{algType, IIP, uType, tType} sol::SolType u::uType @@ -133,13 +133,14 @@ mutable struct ODEIntegrator{algType <: Union{OrdinaryDiffEqAlgorithm, DAEAlgori opts::O stats::DiffEqBase.Stats initializealg::IA + differential_vars::DV fsalfirst::FSALType fsallast::FSALType function ODEIntegrator{algType, IIP, uType, duType, tType, pType, eigenType, EEstT, tTypeNoUnits, tdirType, ksEltype, SolType, F, CacheType, O, FSALType, EventErrorType, CallbackCacheType, - InitializeAlgType}(sol, u, du, k, t, dt, f, p, uprev, uprev2, + InitializeAlgType, DV}(sol, u, du, k, t, dt, f, p, uprev, uprev2, duprev, tprev, alg, dtcache, dtchangeable, dtpropose, tdir, eigen_est, EEst, qold, q11, erracc, dtacc, @@ -154,7 +155,7 @@ mutable struct ODEIntegrator{algType <: Union{OrdinaryDiffEqAlgorithm, DAEAlgori accept_step, isout, reeval_fsal, u_modified, reinitialize, isdae, opts, stats, - initializealg) where {algType, IIP, uType, + initializealg, differential_vars) where {algType, IIP, uType, duType, tType, pType, eigenType, EEstT, tTypeNoUnits, tdirType, @@ -163,10 +164,10 @@ mutable struct ODEIntegrator{algType <: Union{OrdinaryDiffEqAlgorithm, DAEAlgori FSALType, EventErrorType, CallbackCacheType, - InitializeAlgType} + InitializeAlgType, DV} new{algType, IIP, uType, duType, tType, pType, eigenType, EEstT, tTypeNoUnits, tdirType, ksEltype, SolType, - F, CacheType, O, FSALType, EventErrorType, CallbackCacheType, InitializeAlgType, + F, CacheType, O, FSALType, EventErrorType, CallbackCacheType, InitializeAlgType, DV }(sol, u, du, k, t, dt, f, p, uprev, uprev2, duprev, tprev, alg, dtcache, dtchangeable, dtpropose, tdir, eigen_est, EEst, qold, q11, erracc, dtacc, success_iter, @@ -175,7 +176,7 @@ mutable struct ODEIntegrator{algType <: Union{OrdinaryDiffEqAlgorithm, DAEAlgori do_error_check, event_last_time, vector_event_last_time, last_event_error, accept_step, isout, reeval_fsal, u_modified, reinitialize, isdae, - opts, stats, initializealg) # Leave off fsalfirst and last + opts, stats, initializealg, differential_vars) # Leave off fsalfirst and last end end diff --git a/src/interp_func.jl b/src/interp_func.jl index 0f43051a87..d6be19dddb 100644 --- a/src/interp_func.jl +++ b/src/interp_func.jl @@ -1,7 +1,7 @@ abstract type OrdinaryDiffEqInterpolation{cacheType} <: DiffEqBase.AbstractDiffEqInterpolation end -struct InterpolationData{F, uType, tType, kType, cacheType} <: +struct InterpolationData{F, uType, tType, kType, cacheType, DV} <: OrdinaryDiffEqInterpolation{cacheType} f::F timeseries::uType @@ -9,9 +9,10 @@ struct InterpolationData{F, uType, tType, kType, cacheType} <: ks::kType dense::Bool cache::cacheType + differential_vars::DV end -struct CompositeInterpolationData{F, uType, tType, kType, cacheType} <: +struct CompositeInterpolationData{F, uType, tType, kType, cacheType, DV} <: OrdinaryDiffEqInterpolation{cacheType} f::F timeseries::uType @@ -20,6 +21,7 @@ struct CompositeInterpolationData{F, uType, tType, kType, cacheType} <: alg_choice::Vector{Int} dense::Bool cache::cacheType + differential_vars::DV end function DiffEqBase.interp_summary(interp::OrdinaryDiffEqInterpolation{ diff --git a/src/misc_utils.jl b/src/misc_utils.jl index 02d3d7c626..0c47a735e7 100644 --- a/src/misc_utils.jl +++ b/src/misc_utils.jl @@ -171,3 +171,29 @@ macro fold(arg) esc(:(Base.@assume_effects :foldable $arg)) end end + +struct DifferentialVarsUndefined end + +""" + get_differential_vars(f, idxs, timeseries::uType) + +Returns an array of booleans for which values are the differential variables +vs algebraic variables. Returns `nothing` for the cases where all variables +are differential variables. Returns `DifferentialVarsUndefined` if it cannot +be determined (i.e. the mass matrix is not diagonal). +""" +function get_differential_vars(f, u) + differential_vars = nothing + if hasproperty(f, :mass_matrix) + mm = f.mass_matrix + mm = mm isa MatrixOperator ? mm.A : mm + + if mm isa UniformScaling + return nothing + elseif !(mm isa SciMLOperators.AbstractSciMLOperator) && isdiag(mm) + differential_vars = reshape(diag(mm) .!= 0, size(u)) + else + return DifferentialVarsUndefined() + end + end +end \ No newline at end of file diff --git a/src/perform_step/nordsieck_perform_step.jl b/src/perform_step/nordsieck_perform_step.jl index 56c720ba1f..fa3b54f8bf 100644 --- a/src/perform_step/nordsieck_perform_step.jl +++ b/src/perform_step/nordsieck_perform_step.jl @@ -14,7 +14,7 @@ function initialize!(integrator, cache::AN5ConstantCache) end @muladd function perform_step!(integrator, cache::AN5ConstantCache, repeat_step = false) - @unpack t, dt, uprev, u, f, p = integrator + @unpack t, dt, uprev, u, f, p, differential_vars = integrator @unpack z, l, m, c_LTE, dts, tsit5tab = cache # handle callbacks, rewind back to order one. if integrator.u_modified @@ -28,11 +28,11 @@ end z[1] = integrator.uprev z[2] = integrator.k[1] * dt z[3] = ode_interpolant(t, dt, nothing, nothing, integrator.k, tsit5tab, nothing, - Val{2}) * dt^2 / 2 + Val{2}, differential_vars) * dt^2 / 2 z[4] = ode_interpolant(t, dt, nothing, nothing, integrator.k, tsit5tab, nothing, - Val{3}) * dt^3 / 6 + Val{3}, differential_vars) * dt^3 / 6 z[5] = ode_interpolant(t, dt, nothing, nothing, integrator.k, tsit5tab, nothing, - Val{4}) * dt^4 / 24 + Val{4}, differential_vars) * dt^4 / 24 z[6] = zero(cache.z[6]) fill!(dts, dt) perform_predict!(cache) @@ -105,7 +105,7 @@ function initialize!(integrator, cache::AN5Cache) end @muladd function perform_step!(integrator, cache::AN5Cache, repeat_step = false) - @unpack t, dt, uprev, u, f, p, uprev2 = integrator + @unpack t, dt, uprev, u, f, p, uprev2, differential_vars = integrator @unpack z, l, m, c_LTE, dts, tmp, ratetmp, atmp, tsit5cache = cache # handle callbacks, rewind back to order one. if integrator.u_modified @@ -120,11 +120,11 @@ end @.. broadcast=false z[1]=integrator.uprev @.. broadcast=false z[2]=integrator.k[1] * dt ode_interpolant!(z[3], t, dt, nothing, nothing, integrator.k, tsit5cache, nothing, - Val{2}) + Val{2}, differential_vars) ode_interpolant!(z[4], t, dt, nothing, nothing, integrator.k, tsit5cache, nothing, - Val{3}) + Val{3}, differential_vars) ode_interpolant!(z[5], t, dt, nothing, nothing, integrator.k, tsit5cache, nothing, - Val{4}) + Val{4}, differential_vars) @.. broadcast=false z[3]=z[3] * dt^2 / 2 @.. broadcast=false z[4]=z[4] * dt^3 / 6 @.. broadcast=false z[5]=z[5] * dt^4 / 24 @@ -198,7 +198,7 @@ function initialize!(integrator, cache::JVODEConstantCache) end @muladd function perform_step!(integrator, cache::JVODEConstantCache, repeat_step = false) - @unpack t, dt, uprev, u, f, p = integrator + @unpack t, dt, uprev, u, f, p, differential_vars = integrator @unpack z, l, m, c_LTE, dts, tsit5tab = cache # handle callbacks, rewind back to order one. if integrator.u_modified || integrator.iter == 1 @@ -265,7 +265,7 @@ function initialize!(integrator, cache::JVODECache) end @muladd function perform_step!(integrator, cache::JVODECache, repeat_step = false) - @unpack t, dt, uprev, u, f, p, uprev2 = integrator + @unpack t, dt, uprev, u, f, p, uprev2, differential_vars = integrator @unpack z, l, m, c_LTE, dts, tmp, ratetmp, atmp, tsit5cache = cache # handle callbacks, rewind back to order one. if integrator.u_modified || integrator.iter == 1 diff --git a/src/solve.jl b/src/solve.jl index 5d73383239..ecf5672cd8 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -404,15 +404,16 @@ function DiffEqBase.__init(prob::Union{DiffEqBase.AbstractODEProblem, stop_at_next_tstop) stats = DiffEqBase.Stats(0) + differential_vars = prob isa DAEProblem ? prob.differential_vars : get_differential_vars(f, u) if _alg isa OrdinaryDiffEqCompositeAlgorithm - id = CompositeInterpolationData(f, timeseries, ts, ks, alg_choice, dense, cache) + id = CompositeInterpolationData(f, timeseries, ts, ks, alg_choice, dense, cache, differential_vars) sol = DiffEqBase.build_solution(prob, _alg, ts, timeseries, dense = dense, k = ks, interp = id, alg_choice = alg_choice, calculate_error = false, stats = stats) else - id = InterpolationData(f, timeseries, ts, ks, dense, cache) + id = InterpolationData(f, timeseries, ts, ks, dense, cache, differential_vars) sol = DiffEqBase.build_solution(prob, _alg, ts, timeseries, dense = dense, k = ks, interp = id, calculate_error = false, stats = stats) @@ -469,7 +470,7 @@ function DiffEqBase.__init(prob::Union{DiffEqBase.AbstractODEProblem, FType, cacheType, typeof(opts), fsal_typeof(_alg, rate_prototype), typeof(last_event_error), typeof(callback_cache), - typeof(initializealg)}(sol, u, du, k, t, tType(dt), f, p, + typeof(initializealg), typeof(differential_vars)}(sol, u, du, k, t, tType(dt), f, p, uprev, uprev2, duprev, tprev, _alg, dtcache, dtchangeable, dtpropose, tdir, eigen_est, EEst, @@ -485,7 +486,7 @@ function DiffEqBase.__init(prob::Union{DiffEqBase.AbstractODEProblem, last_event_error, accept_step, isout, reeval_fsal, u_modified, reinitiailize, isdae, - opts, stats, initializealg) + opts, stats, initializealg, differential_vars) if initialize_integrator if isdae diff --git a/test/interface/algebraic_interpolation.jl b/test/interface/algebraic_interpolation.jl new file mode 100644 index 0000000000..d9b5eba82f --- /dev/null +++ b/test/interface/algebraic_interpolation.jl @@ -0,0 +1,122 @@ +using OrdinaryDiffEq +using Test + +function rober_ip(du, u, p, t) + y₁, y₂, y₃ = u + k₁, k₂, k₃ = p + du[1] = -k₁ * y₁ + k₃ * y₂ * y₃ + du[2] = k₁ * y₁ - k₃ * y₂ * y₃ - k₂ * y₂^2 + du[3] = y₁ + y₂ + y₃ - 1 + nothing +end +function rober_op(u, p, t) + du = similar(u) + y₁, y₂, y₃ = u + k₁, k₂, k₃ = p + du[1] = -k₁ * y₁ + k₃ * y₂ * y₃ + du[2] = k₁ * y₁ - k₃ * y₂ * y₃ - k₂ * y₂^2 + du[3] = y₁ + y₂ + y₃ - 1 + du +end +M = [1.0 0 0 + 0 1.0 0 + 0 0 0] +fip = ODEFunction(rober_ip, mass_matrix = M) +prob_ip = ODEProblem(fip, [1.0, 0.0, 0.0], (0.0, 1e5), (0.04, 3e7, 1e4)) + +fop = ODEFunction(rober_op, mass_matrix = M) +prob_op = ODEProblem(fop, [1.0, 0.0, 0.0], (0.0, 1e5), (0.04, 3e7, 1e4)) + +ref_ip = solve(prob_ip, Rodas5P(), reltol = 1e-8, abstol = 1e-8) +ref_op = solve(prob_op, Rodas5P(), reltol = 1e-8, abstol = 1e-8) + +sol_ip = solve(prob_ip, FBDF(), reltol = 1e-8, abstol = 1e-8) +sol_op = solve(prob_op, FBDF(), reltol = 1e-8, abstol = 1e-8) + +# make sure interpolation changes don't accidentally break this test suite +# the intention is that ref uses a stiffness-aware interpolation, while sol uses hermite +@test DiffEqBase.interp_summary(sol_ip) == "3rd order Hermite" +@test DiffEqBase.interp_summary(sol_op) == "3rd order Hermite" +@test occursin("stiffness-aware", DiffEqBase.interp_summary(ref_ip)) +@test occursin("stiffness-aware", DiffEqBase.interp_summary(ref_op)) + +reltol = 1e-4 +abstol = 1e-4 +t = 1 +tv = [1, 10, 100] +idxs = 3 +idxsv = [2, 3] + +# primal, no index +@test isapprox(ref_ip(t), sol_ip(t), rtol = reltol, atol = abstol) # ip, t +@test isapprox(ref_ip(tv), sol_ip(tv), rtol = reltol, atol = abstol) # ip, tv +@test isapprox(ref_op(t), sol_op(t), rtol = reltol, atol = abstol) # op, t +@test isapprox(ref_op(tv), sol_op(tv), rtol = reltol, atol = abstol) # op, tv + +# primal, scalar index +@test isapprox(ref_ip(t, idxs=idxs), sol_ip(t, idxs=idxs), rtol = reltol, atol = abstol) # ip, t +@test isapprox(ref_ip(tv, idxs=idxs), sol_ip(tv, idxs=idxs), rtol = reltol, atol = abstol) # ip, tv +@test isapprox(ref_op(t, idxs=idxs), sol_op(t, idxs=idxs), rtol = reltol, atol = abstol) # op, t +@test isapprox(ref_op(tv, idxs=idxs), sol_op(tv, idxs=idxs), rtol = reltol, atol = abstol) # op, tv + +# primal, vector index +@test isapprox(ref_ip(t, idxs=idxsv), sol_ip(t, idxs=idxsv), rtol = reltol, atol = abstol) +@test isapprox(ref_ip(tv, idxs=idxsv), sol_ip(tv, idxs=idxsv), rtol = reltol, atol = abstol) +@test isapprox(ref_op(t, idxs=idxsv), sol_op(t, idxs=idxsv), rtol = reltol, atol = abstol) +@test isapprox(ref_op(tv, idxs=idxsv), sol_op(tv, idxs=idxsv), rtol = reltol, atol = abstol) + +abstol = 1e-3 +# derivative, no index +@test isapprox(ref_ip(t, Val{1}), sol_ip(t, Val{1}), rtol = reltol, atol = abstol) +@test isapprox(ref_ip(tv, Val{1}), sol_ip(tv, Val{1}), rtol = reltol, atol = abstol) +@test isapprox(ref_op(t, Val{1}), sol_op(t, Val{1}), rtol = reltol, atol = abstol) +@test isapprox(ref_op(tv, Val{1}), sol_op(tv, Val{1}), rtol = reltol, atol = abstol) + +# derivative, scalar index +@test isapprox(ref_ip(t, Val{1}, idxs=idxs), sol_ip(t, Val{1}, idxs=idxs), rtol = reltol, atol = abstol) +@test isapprox(ref_ip(tv, Val{1}, idxs=idxs), sol_ip(tv, Val{1}, idxs=idxs), rtol = reltol, atol = abstol) +@test isapprox(ref_op(t, Val{1}, idxs=idxs), sol_op(t, Val{1}, idxs=idxs), rtol = reltol, atol = abstol) +@test isapprox(ref_op(tv, Val{1}, idxs=idxs), sol_op(tv, Val{1}, idxs=idxs), rtol = reltol, atol = abstol) + +# derivative, vector index +@test isapprox(ref_ip(tv, Val{1}, idxs=idxsv), sol_ip(tv, Val{1}, idxs=idxsv), rtol = reltol, atol = abstol) +@test isapprox(ref_ip(t, Val{1}, idxs=idxsv), sol_ip(t, Val{1}, idxs=idxsv), rtol = reltol, atol = abstol) +@test isapprox(ref_op(t, Val{1}, idxs=idxsv), sol_op(t, Val{1}, idxs=idxsv), rtol = reltol, atol = abstol) +@test isapprox(ref_op(tv, Val{1}, idxs=idxsv), sol_op(tv, Val{1}, idxs=idxsv), rtol = reltol, atol = abstol) + +# higher derivatives should be zero +# second derivative, no index +@test (sol_ip(t, Val{2}) .== 0) == [false, false, true] +@test all(map(==([false, false, true]), sol_ip(tv, Val{2}) .== 0)) +@test (sol_op(t, Val{2}) .== 0) == [false, false, true] +@test all(map(==([false, false, true]), sol_op(tv, Val{2}) .== 0)) + +# second derivative, scalar index +@test sol_ip(t, Val{2}, idxs=idxs) == 0 +@test all(sol_ip(tv, Val{2}, idxs=idxs) .== 0) +@test sol_op(t, Val{2}, idxs=idxs) == 0 +@test all(sol_op(tv, Val{2}, idxs=idxs) .== 0) + +# second derivative, vector index +@test (sol_ip(t, Val{2}, idxs=idxsv) .== 0) == [false, true] +@test all(map(==([false, true]), sol_ip(tv, Val{2}, idxs=idxsv) .== 0)) +@test (sol_op(t, Val{2}, idxs=idxsv) .== 0) == [false, true] +@test all(map(==([false, true]), sol_op(tv, Val{2}, idxs=idxsv) .== 0)) + +# third derivative, no index +@test (sol_ip(t, Val{3}) .== 0) == [false, false, true] +@test all(map(==([false, false, true]), sol_ip(tv, Val{3}) .== 0)) +@test (sol_op(t, Val{3}) .== 0) == [false, false, true] +@test all(map(==([false, false, true]), sol_op(tv, Val{3}) .== 0)) + +# third derivative, scalar index +@test sol_ip(t, Val{3}, idxs=idxs) == 0 +@test all(sol_ip(tv, Val{3}, idxs=idxs) .== 0) +@test sol_op(t, Val{3}, idxs=idxs) == 0 +@test all(sol_op(tv, Val{3}, idxs=idxs) .== 0) + +# third derivative, vector index +@test (sol_ip(t, Val{3}, idxs=idxsv) .== 0) == [false, true] +@test all(map(==([false, true]), sol_ip(tv, Val{3}, idxs=idxsv) .== 0)) +@test (sol_op(t, Val{3}, idxs=idxsv) .== 0) == [false, true] +@test all(map(==([false, true]), sol_op(tv, Val{3}, idxs=idxsv) .== 0)) \ No newline at end of file diff --git a/test/interface/inplace_interpolation.jl b/test/interface/inplace_interpolation.jl index b5ea6459ca..88cdd084c3 100644 --- a/test/interface/inplace_interpolation.jl +++ b/test/interface/inplace_interpolation.jl @@ -35,6 +35,6 @@ out_VMF = vecarrzero(ntt, size(prob_ode_2Dlinear.u0)) # Vector{Matrix{Float64} @test sol_ODE_2D(out_VVF_1, tt; idxs = 3:3) isa Vector{Vector{Float64}} @test sol_ODE_2D(out_VVF_2, tt; idxs = 2:3) isa Vector{Vector{Float64}} @test sol_ODE_2D(out_VMF, tt) isa Vector{Matrix{Float64}} - @test sol_ODE_2D_interp.u == out_VMF + @test sol_ODE_2D_interp.u ≈ out_VMF end end diff --git a/test/interface/noindex_tests.jl b/test/interface/noindex_tests.jl index d9c4107cd0..ebb646292f 100644 --- a/test/interface/noindex_tests.jl +++ b/test/interface/noindex_tests.jl @@ -25,6 +25,8 @@ function Base.similar(bc::Base.Broadcast.Broadcasted{NoIndexStyle{N}}, end Base.Broadcast._broadcast_getindex(x::NoIndexArray, i) = x.x[i] Base.Broadcast.extrude(x::NoIndexArray) = x +using ArrayInterface +ArrayInterface.fast_scalar_indexing(::Type{<:NoIndexArray}) = false @inline function Base.copyto!(dest::NoIndexArray, bc::Base.Broadcast.Broadcasted{<:NoIndexStyle}) diff --git a/test/runtests.jl b/test/runtests.jl index 055edc4fdf..a893a6294b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -38,6 +38,7 @@ end @time @safetestset "Type Handling Tests" include("interface/type_handling.jl") @time @safetestset "Controller Tests" include("interface/controllers.jl") @time @safetestset "Inplace Interpolation Tests" include("interface/inplace_interpolation.jl") + @time @safetestset "Algebraic Interpolation Tests" include("interface/algebraic_interpolation.jl") end if !is_APPVEYOR && (GROUP == "All" || GROUP == "InterfaceII" || GROUP == "Interface")