Closed SKopecz closed 2 months ago
Todo:
d_prototype
and use it in all MPRK algorithms.d_prototype
.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 P
s and D
s in the algorithm caches have the same types as p_protoype
and d_prototype
. I suggest to add corresponding @assert
s 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
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
That's an interesting option! I'll use this.
See #104
Just like we use
p_prototype
to define the type of the production matrices used in the algorithms, there should be an analogd_prototype
to define the type of the destruction vector inPDSProblem
s.Otherwise the behavior of
PDSProblem
does not match the description in the docstring.