Use d_prototype in PDSProblems

[SKopecz]

Just like we use p_prototype to define the type of the production matrices used in the algorithms, there should be an analog d_prototype to define the type of the destruction vector in PDSProblems.

Otherwise the behavior of PDSProblem does not match the description in the docstring.

[SKopecz]

Todo:

• Implement d_prototype and use it in all MPRK algorithms.
• State in doc strings that copies of these prototypes are in within algorithms
• Add tests for d_prototype.

[SKopecz]

The following code defines a PDSProblem with given sparse p_prototype and d_prototype. The code runs as expected.

But I don't think that this qualifies as a test for d_prototype (or p_prototype). The prototypes could be completely ignored internally.

julia> using PositiveIntegrators, SparseArrays

julia> N = 100; # number of subintervals

julia> dx = 1/N; # mesh width

julia> x = LinRange(dx, 1.0, N); # discretization points x_1,...,x_N = x_0

julia> u0 = 0.0 .+ (0.4 .≤ x .≤ 0.6) .* 1.0; # initial solution

julia> tspan = (0.0, 1.0); # time domain

julia> function lin_adv_P!(P, u, p, t)
           P .= 0.0
           N = length(u)
           dx = 1 / N
           P[1, N] = u[N] / dx
           for i in 2:N
               P[i, i - 1] = u[i - 1] / dx
           end
           return nothing
       end
lin_adv_P! (generic function with 1 method)

julia> function lin_adv_D!(D, u, p, t)
           D .= 0.0
           return nothing
       end
lin_adv_D! (generic function with 1 method)

julia> prob = PDSProblem(lin_adv_P!, lin_adv_D!, u0, tspan); # create the PDS

julia> p_prototype = spdiagm(-1 => ones(eltype(u0), N - 1), N - 1 => ones(eltype(u0), 1));

julia> d_prototype = spzeros(N);

julia> prob_sparse = PDSProblem(lin_adv_P!, lin_adv_D!, u0, tspan; p_prototype=p_prototype, d_prototype = d_prototype);

julia> sol = solve(prob, MPRK43I(1.0, 0.5); save_everystep = false);

julia> sol_sparse = solve(prob_sparse, MPRK43I(1.0, 0.5); save_everystep = false);

julia> @assert sol_sparse.t ≈ sol.t

julia> @assert sol_sparse.u ≈ sol.u

What would be proper tests?

As I see it, we want to ensure that the Ps and Ds in the algorithm caches have the same types as p_protoype and d_prototype. I suggest to add corresponding @asserts in the alg_cache function of each algorithm. In this case the above code would be a proper test. Otherwise we would need to create the algorithm caches manually in the tests and check the filed types.

@ranocha: What do you think?

Otherwise we would have to

[ranocha]

I think it's fine to use something like the test that you suggest above. Then, we can write the test like

julia> function lin_adv_D!(D, u, p, t)
           @test D isa SparseVector
           D .= 0.0
           return nothing
       end

[SKopecz]

That's an interesting option! I'll use this.

[ranocha]

See #104