fix 2N adaptive step size methods

[ranocha]

The first part fixes one of my typos from #1278 (f(k, u, p, t) instead of f(k, u, p, t+dt)).

As far as I understand, setting @.. tmp = dt*k after a step is only correct if the next step size dt will be the same - or do I misunderstand something here, @ChrisRackauckas, @kanav99 ?

While there are no 2N methods with error-based step size adaptation currently in DiffEq, users could vary the step size via callbacks, which should definitely be supported, e.g. for CFL-based step size control.

[ChrisRackauckas]

No that should be safe for dt changes since k is independent of dt.

[ranocha]

If we set @.. tmp = dt_n * k at the end of step n but want to use @.. tmp = dt_{n+1} * k at the beginning of step n+1, there seems to be something wrong? Okay, k is independent of dt, but @.. u = u + B1*tmp at the beginning of step n+1 should use dt_{n+1}, shouldn't it?

[ChrisRackauckas]

That's what the FSAL correction codes handle.

[ranocha]

But the low-storage methods are not marked as FSAL?

OrdinaryDiffEq.jl/src/alg_utils.jl

Line 17 in 821fc83

isfsal(alg::CarpenterKennedy2N54) = false

Anyway, if you say the current setup is correct, I need to believe you (or find a counterexample).

[ChrisRackauckas]

Oh then that would be an issue. It would need an if u_modified!... but really it should be marked as FSAL if it's doing this.

[ranocha]

Okay, so here is a counterexample:

julia> using Revise, Test, OrdinaryDiffEq, DiffEqCallbacks

julia> function f(du, u, p, t)
           du .= t
           nothing
       end

julia> u0 = zeros(2)

julia> function affect!(integrator)
         dt = rand((0.0625, 0.125, 0.25))
         set_proposed_dt!(integrator, dt)
         integrator.opts.dtmax = dt
         integrator.dtcache = dt
         u_modified!(integrator, false)
         nothing
       end

julia> condition_discrete(u, t, integrator) = true

julia> callback_discrete = DiscreteCallback(condition_discrete, affect!, save_positions=(false,false));

julia> tspan = (0.0, 2.0)

julia> prob = ODEProblem(f, u0, tspan)

julia> sol = solve(prob, Tsit5(), dt=0.125, adaptive=false, callback=callback_discrete)

julia> @test all(idx -> sol.u[idx] ≈ 0.5*sol.t[idx]^2*ones(length(sol.u[idx])), eachindex(sol.t))

julia> sol = solve(prob, CarpenterKennedy2N54(williamson_condition=false), dt=0.125, adaptive=false, callback=callback_discrete)

julia> @test all(idx -> sol.u[idx] ≈ 0.5*sol.t[idx]^2*ones(length(sol.u[idx])), eachindex(sol.t))

Master:

julia> sol = solve(prob, CarpenterKennedy2N54(williamson_condition=false), dt=0.125, adaptive=false, callback=callback_discrete)
retcode: Success
Interpolation: 3rd order Hermite
t: 33-element Array{Float64,1}:
 0.0
 ...
 2.0
u: 33-element Array{Array{Float64,1},1}:
 [0.0, 0.0]
 ...
 [2.0012237287245322, 2.0012237287245322]

julia> @test all(idx -> sol.u[idx] ≈ 0.5*sol.t[idx]^2*ones(length(sol.u[idx])), eachindex(sol.t))
Test Failed at REPL[12]:1
  Expression: all((idx->begin
            sol.u[idx] ≈ 0.5 * sol.t[idx] ^ 2 * ones(length(sol.u[idx]))
        end), eachindex(sol.t))
ERROR: There was an error during testing

This PR:

julia> sol = solve(prob, CarpenterKennedy2N54(williamson_condition=false), dt=0.125, adaptive=false, callback=callback_discrete)
retcode: Success
Interpolation: 3rd order Hermite
t: 21-element Array{Float64,1}:
 0.0
 ...
 2.0
u: 21-element Array{Array{Float64,1},1}:
 [0.0, 0.0]
 ...
 [2.0, 2.0]

julia> @test all(idx -> sol.u[idx] ≈ 0.5*sol.t[idx]^2*ones(length(sol.u[idx])), eachindex(sol.t))
Test Passed

This RHS has to be solved exactly for a fourth-order accurate method. Note that I'm not talking about the modification of u but the modification of dt. I've used the code from the StepsizeLimiter callback in the callback library.

[ChrisRackauckas]

It would be good to figure out why this is not using the standard FSAL tools though, since this should be all handled by that mechanism. That is why the mechanism is there.

[ranocha]

Could you open an issue for that or ping people here? At least the line setting isfsal(alg::CarpenterKennedy2N54) = false was changed in a huge commit ed420fb last year... When I added the first 2N methods in c5aa1b0, there was isfsal(alg::OrdinaryDiffEqAlgorithm) = true. The commit introducing the wrong way of setting tmp is 1923618 by @kanav99.

[ChrisRackauckas]

@kanav99 's stuff in #691 and #1185 I think was the cause of this. I guess we can use a workaround for now and open an issue to clean up the rest of this with the incrementing ODE form.

[ranocha]

Sounds good to me 👍

It looks like the 2N methods haven't been used for quite a long time...

[ChrisRackauckas]

@kanav99 was getting them fully minimal memory so they could be used in Clima stuff. That'll be the real users of it. It's always going to be PDE folks using that portion of the library of course, so hooking it up to real PDE applications is what will complete them.

[ranocha]

Well, that's the reason I'm fixing errors here: We're using them in our hyperbolic PDE applications in Trixi.jl.

[ChrisRackauckas]

Yeah. I'll merge if you get tests passing so you can get your application going, but @kanav99 should take a deeper look soon.

[ranocha]

Hmm, looks like I only needed to mark more tests as passing 🎉

[ChrisRackauckas]

That's always good.

[kanav99]

I will take a deep look today, but this looks good to me at a quick glance.

[ranocha]

The build error on Travis

curl: (22) The requested URL returned error: 502 Bad Gateway
Downloading artifact: CompilerSupportLibraries

is unrelated to this PR, isn't it?

[ChrisRackauckas]

That just means it lost internet while downloading. No big deal and unrelated tests.