Questions About Saveat

When I run the following code with saveat I get an error:

using DifferentialEquations
function f(du,u,p,t)
  du[1] = 1.01*u[1]
  du[2] = 2.01*u[2]
end
ICs = [.1 3];
tspan = (0.0f0,5.0f0);
SaveTimes = [0.1:.33:4.0];
prob = ODEProblem(f,ICs,tspan);
function prob_func(prob,i,repeat)
  remake(prob,u0=ICs*rand())
end
ensemble_prob = EnsembleProblem(prob, prob_func = prob_func, safetycopy=false)
@time sim = solve(ensemble_prob,Tsit5(),EnsembleDistributed(),trajectories=20, saveat=SaveTimes)

I’m having interpreting the stack trace:

ERROR: LoadError: MethodError: no method matching isless(::Float32, ::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64})
Closest candidates are:
  isless(::T, ::T) where T<:Union{Float16, Float32, Float64} at float.jl:424
  isless(::AbstractVector, ::AbstractVector) at abstractarray.jl:2610
  isless(::AbstractFloat, ::AbstractFloat) at operators.jl:184
  ...
Stacktrace:
  [1] <(x::Float32, y::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64})
    @ Base ./operators.jl:356
  [2] initialize_saveat
    @ ~/.julia/packages/OrdinaryDiffEq/eqC5v/src/solve.jl:562 [inlined]
  [3] __init(prob::ODEProblem{Matrix{Float64}, Tuple{Float32, Float32}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Vector{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float32, dtmin::Nothing, dtmax::Float32, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{Int64}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/eqC5v/src/solve.jl:205
  [4] #__solve#562
    @ ~/.julia/packages/OrdinaryDiffEq/eqC5v/src/solve.jl:4 [inlined]
  [5] #solve_call#28
    @ ~/.julia/packages/DiffEqBase/KouNZ/src/solve.jl:437 [inlined]
  [6] #solve_up#34
    @ ~/.julia/packages/DiffEqBase/KouNZ/src/solve.jl:780 [inlined]
  [7] #solve#33
    @ ~/.julia/packages/DiffEqBase/KouNZ/src/solve.jl:760 [inlined]
  [8] batch_func(i::Int64, prob::EnsembleProblem{ODEProblem{Matrix{Float64}, Tuple{Float32, Float32}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, typeof(prob_func), typeof(SciMLBase.DEFAULT_OUTPUT_FUNC), typeof(SciMLBase.DEFAULT_REDUCTION), Nothing}, alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; kwargs::Base.Pairs{Symbol, Vector{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Vector{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}}})
    @ SciMLBase ~/.julia/packages/SciMLBase/BC0Kw/src/ensemble/basic_ensemble_solve.jl:92
  [9] #500
    @ ~/.julia/packages/SciMLBase/BC0Kw/src/ensemble/basic_ensemble_solve.jl:129 [inlined]
 [10] (::Base.var"#929#934"{SciMLBase.var"#500#501"{Base.Pairs{Symbol, Vector{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Vector{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}}}, EnsembleProblem{ODEProblem{Matrix{Float64}, Tuple{Float32, Float32}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, typeof(prob_func), typeof(SciMLBase.DEFAULT_OUTPUT_FUNC), typeof(SciMLBase.DEFAULT_REDUCTION), Nothing}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}})(r::Base.RefValue{Any}, args::Tuple{Int64})
    @ Base ./asyncmap.jl:100
 [11] macro expansion
    @ ./asyncmap.jl:234 [inlined]
 [12] (::Base.var"#945#946"{Base.var"#929#934"{SciMLBase.var"#500#501"{Base.Pairs{Symbol, Vector{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Vector{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}}}, EnsembleProblem{ODEProblem{Matrix{Float64}, Tuple{Float32, Float32}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, typeof(prob_func), typeof(SciMLBase.DEFAULT_OUTPUT_FUNC), typeof(SciMLBase.DEFAULT_REDUCTION), Nothing}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}}, Channel{Any}, Nothing})()
    @ Base ./task.jl:484
Stacktrace:
  [1] (::Base.var"#939#941")(x::Task)
    @ Base ./asyncmap.jl:177
  [2] foreach(f::Base.var"#939#941", itr::Vector{Any})
    @ Base ./abstractarray.jl:2772
  [3] maptwice(wrapped_f::Function, chnl::Channel{Any}, worker_tasks::Vector{Any}, c::UnitRange{Int64})
    @ Base ./asyncmap.jl:177
  [4] wrap_n_exec_twice
    @ ./asyncmap.jl:153 [inlined]
  [5] #async_usemap#924
    @ ./asyncmap.jl:103 [inlined]
  [6] #asyncmap#923
    @ ./asyncmap.jl:81 [inlined]
  [7] pmap(f::Function, p::Distributed.CachingPool, c::UnitRange{Int64}; distributed::Bool, batch_size::Int64, on_error::Nothing, retry_delays::Vector{Any}, retry_check::Nothing)
    @ Distributed /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Distributed/src/pmap.jl:126
  [8] solve_batch(prob::EnsembleProblem{ODEProblem{Matrix{Float64}, Tuple{Float32, Float32}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, typeof(prob_func), typeof(SciMLBase.DEFAULT_OUTPUT_FUNC), typeof(SciMLBase.DEFAULT_REDUCTION), Nothing}, alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, ensemblealg::EnsembleDistributed, II::UnitRange{Int64}, pmap_batch_size::Int64; kwargs::Base.Pairs{Symbol, Vector{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Vector{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}}})
    @ SciMLBase ~/.julia/packages/SciMLBase/BC0Kw/src/ensemble/basic_ensemble_solve.jl:128
  [9] macro expansion
    @ ./timing.jl:381 [inlined]
 [10] __solve(prob::EnsembleProblem{ODEProblem{Matrix{Float64}, Tuple{Float32, Float32}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, typeof(prob_func), typeof(SciMLBase.DEFAULT_OUTPUT_FUNC), typeof(SciMLBase.DEFAULT_REDUCTION), Nothing}, alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, ensemblealg::EnsembleDistributed; trajectories::Int64, batch_size::Int64, pmap_batch_size::Int64, kwargs::Base.Pairs{Symbol, Vector{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Vector{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}}})
    @ SciMLBase ~/.julia/packages/SciMLBase/BC0Kw/src/ensemble/basic_ensemble_solve.jl:56
 [11] #solve#35
    @ ~/.julia/packages/DiffEqBase/KouNZ/src/solve.jl:796 [inlined]
 [12] top-level scope
    @ ./timing.jl:263
 [13] include(fname::String)
    @ Base.MainInclude ./client.jl:476
 [14] top-level scope
    @ REPL[1]:1
in expression starting at /Users/yoshi/Dropbox/xtest2.jl:14

Also, how does the integrator work? Does it dynamically create the appropriate stepsize, and along the way, makes sure to save at the particular values specified? I’m a little confused about how to interpret the object that solve returns.

You don’t need to put the range in a vector. Try

SaveTimes = 0.1:.33:4.0;

The error messages says there is no way to compute if a Float32 is less than a StepRange during a search for the first save point.

Here is the documentation for the solution object, and here is a description of how step sizes are selected during integration.

2 Likes