Using a custom Number type that throws an error when it is divided by zero (as opposed to returning Inf as Float64 does), an error is thrown during the initialization of the derivatives vector for non-fsal integrators (e.g. Vern integrators).
The derivatives vector is initialized with the correct size and data type.
Minimal Reproducible Example 👇
using OrdinaryDiffEq
using DACE
function kepler_ode!(du, u, p, t)
du[1:3] .= u[4:6]
du[4:6] .= -u[1:3] * p ./ sum(u[1:3].^2)^(3/2)
end
DACE.init(2,6)
xs = [1.0, 0.0, 0.0, 0.0, 1.0, 0.0] .+ [DA(i,1.0) for i in 1:6]
prob = ODEProblem(kepler_ode!, xs, [0.0, 2π], p=1.0)
sol = solve(prob, Vern7(), abstol=1e-12, reltol=1e-12)
Describe the bug 🐞
Using a custom
Number
type that throws an error when it is divided by zero (as opposed to returningInf
asFloat64
does), an error is thrown during the initialization of the derivatives vector for non-fsal integrators (e.g. Vern integrators).This is caused by line 51 in src/initdt.jl, i.e.:
and can be overcome by replacing the above with:
or
oneunit_tType
I guess.Expected behavior
The derivatives vector is initialized with the correct size and data type.
Minimal Reproducible Example 👇
Error & Stacktrace ⚠️
Environment (please complete the following information):
using Pkg; Pkg.status()
The DACE package is obtained from https://github.com/a-ev/DACE.jl
using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)
versioninfo()