Implement simple "controller" for adaptive time stepping

[SKopecz]

In most papers about MPRK schemes the Robertson problem is "adaptively" solved by setting $\Delta t_{i+1}=c\Delta t_i$ for some constant $c$. We should also have this option for comparisons.

[ranocha]

That's already possible right now with a callback:

julia> using OrdinaryDiffEq

julia> ode = ODEProblem((du, u, p, t) -> du .= -u, [1.0], (0.0, 1.0))
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 1.0)
u0: 1-element Vector{Float64}:
 1.0

julia> sol = solve(ode, ImplicitEuler(); 
                   adaptive = false, dt = 1 / 2^5, 
                   callback = DiscreteCallback(Returns(true), 
                       integrator -> set_proposed_dt!(integrator, 2 * get_proposed_dt(integrator));
                       save_positions = (false, false)));

julia> diff(sol.t)
6-element Vector{Float64}:
 0.03125
 0.0625
 0.125
 0.25
 0.5
 0.03125

But it would definitely be nice to document this explicitly.

[SKopecz]

This approach results in an OutOfMemoryError() when we try to solve prob_pds_robertson. Would be great to make this work. This should also be included in the tests.

julia> using PositiveIntegrators, OrdinaryDiffEq
julia> sol = solve(prob_pds_robertson, ImplicitEuler(); 
                          adaptive = false, dt = 1e-5, 
                          callback = DiscreteCallback(Returns(true), 
                          integrator -> set_proposed_dt!(integrator, 2 * get_proposed_dt(integrator));
                          save_positions = (false, false)));
ERROR: OutOfMemoryError()
Stacktrace:
  [1] sizehint!
    @ ./array.jl:1351 [inlined]
  [2] __init(prob::ODEProblem{…}, alg::ImplicitEuler{…}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{…}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::DiscreteCallback{…}, dense::Bool, calck::Bool, dt::Float64, dtmin::Float64, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{…}, abstol::Nothing, reltol::Nothing, qmin::Rational{…}, qmax::Int64, qsteady_min::Int64, qsteady_max::Rational{…}, beta1::Nothing, beta2::Nothing, qoldinit::Rational{…}, 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), progress_id::Symbol, userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::@Kwargs{})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/J5YGH/src/solve.jl:288
  [3] __init (repeats 5 times)
    @ ~/.julia/packages/OrdinaryDiffEq/J5YGH/src/solve.jl:10 [inlined]
  [4] __solve(::ODEProblem{…}, ::ImplicitEuler{…}; kwargs::@Kwargs{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/J5YGH/src/solve.jl:5
  [5] __solve
    @ ~/.julia/packages/OrdinaryDiffEq/J5YGH/src/solve.jl:1 [inlined]
  [6] solve_call(_prob::ODEProblem{…}, args::ImplicitEuler{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/eTCPy/src/solve.jl:608
  [7] solve_call
    @ DiffEqBase ~/.julia/packages/DiffEqBase/eTCPy/src/solve.jl:566 [inlined]
  [8] #solve_up#42
    @ DiffEqBase ~/.julia/packages/DiffEqBase/eTCPy/src/solve.jl:1057 [inlined]
  [9] solve_up
    @ DiffEqBase ~/.julia/packages/DiffEqBase/eTCPy/src/solve.jl:1043 [inlined]
 [10] #solve#40
    @ DiffEqBase ~/.julia/packages/DiffEqBase/eTCPy/src/solve.jl:980 [inlined]
 [11] top-level scope
    @ REPL[2]:1
Some type information was truncated. Use `show(err)` to see complete types.

[ranocha]

We need a kind of ugly hack here:

julia> sol = solve(prob_pds_robertson, ImplicitEuler(); 
                   adaptive = false, dt = Inf, 
                   callback = DiscreteCallback(Returns(true), 
                       integrator -> set_proposed_dt!(integrator, 2 * get_proposed_dt(integrator));
                       save_positions = (false, false),
                       initialize = (c, u, t, integrator) -> set_proposed_dt!(integrator, 1.0e-5)));

should work. The reason is that OrdinaryDiffEq.jl wants to allocate memory for the solution. With adaptive = false, it divides the time span by the given dt to get the number of steps. This is of course completely wrong in this case since the initial dt is tiny compared to tspan. The hack above avoids this.

[ranocha]

Documented in #110