`MPRK22(α)` for α≠1/2 applied to a `ConservativePDSProblem`?

[JoshuaLampert]

Thank you very much for this nice package! I'm interested in using PositiveIntegrators.jl to simulate Fokker-Planck equations, which can be written in conservative production-destruction form. However, I noticed that MPRK22(α) does not seem to be conservative for $\alpha\neq 1/2$ (while it is conservative for $\alpha = 1/2$). Solving the same system with classical RK methods like Heun() and also with MPE() yields discrete conservatism. You can find the code below. Unfortunately, I was not able to reproduce this for a simpler system.

Code

```julia using PositiveIntegrators using OrdinaryDiffEq using QuadGK using SimpleUnPack: @unpack # parameters I = (-1.0, 1.0) N = 80 w_interfaces = LinRange(first(I), last(I), N + 1) dw = w_interfaces[2] - w_interfaces[1] w_cells = w_interfaces[1:end - 1] .+ dw/2 σ2 = 0.2 # = σ^2 # functions in the model D = @. σ2/2 * (1 - w_interfaces^2)^2 # diffusion term Dprime = @. -2 * σ2 * w_interfaces * (1 - w_interfaces^2) # D', i.e. first derivative of diffusion term # aggregation dynamics operator (use simple quadrature rule to compute integral) # To avoid allocations, also pass `w_cells` and `dw` to avoid global variables function B(f_cells, w_cells, w, dw) res = 0.0 for i in 1:length(f_cells) res += (w - w_cells[i]) * f_cells[i] end res *= dw return res end λ = zeros(N) # on interfaces, first entry will not be used, don't need to allocate for last entry C = zeros(N) # on interfaces, first entry will not be used, don't need to allocate for last entry δ = zeros(N) # on interfaces, first entry will not be used, don't need to allocate for last entry BB = zeros(N) # on interfaces, first entry will not be used, don't need to allocate for last entry u_tilde = zeros(N) fluxes = zeros(N + 1) production = zeros(N, N) # production matrix p = (w_interfaces = w_interfaces, w_cells = w_cells, dw = dw, D = D, Dprime = Dprime, B = B, λ = λ, C = C, δ = δ, BB = BB, u_tilde = u_tilde, fluxes = fluxes, production = production) # initial condition c = 30.0 f0_1(w) = exp(-c*(w + 0.5)^2) + exp(-c*(w - 0.5)^2) oneoverβ, _ = quadgk(f0_1, first(I), last(I), rtol = 1e-12) β = 1/oneoverβ f0(w) = β * f0_1(w) u0 = f0.(w_cells) # stationary solution for reference f0w(w) = w * f0(w) # call this uu instead of u because u is used as unknown in the `ODEProblem` uu, _ = quadgk(f0w, first(I), last(I), rtol = 1e-12) f_stat_1(w) = 1/(1 - w^2)^2 * ((1 + w)/(1 - w))^(uu/(2 * σ2))*exp(-(1 - uu*w)/(σ2 * (1 - w^2))) oneoverK, _ = quadgk(f_stat_1, first(I), last(I), rtol = 1e-12) K = 1/oneoverK f_stat(w) = K * f_stat_1(w) u_stat = f_stat.(w_cells) tspan = (0.0, 10.0) # time span dt = dw^2/(2*σ2) function calculate_values!(BB, λ, C, δ, u_tilde, u, w_cells, w_interfaces, dw, D, Dprime, i) BB[i] = B(u, w_cells, w_interfaces[i], dw) λ[i] = dw * (BB[i] + Dprime[i])/D[i] C[i] = λ[i] * D[i]/dw δ[i] = isapprox(λ[i], 0.0) ? 0.5 : 1/λ[i] - 1/expm1(λ[i]) u_tilde[i] = (1 - δ[i]) * u[i] + δ[i] * u[i - 1] end # Out-of-place implementation of the P matrix for the Fokker-Planck model function P(u, p, t) @unpack w_interfaces, w_cells, dw, D, Dprime, B, λ, C, δ, BB, u_tilde, production = p N = length(w_interfaces) - 1 for i = 2:N calculate_values!(BB, λ, C, δ, u_tilde, u, w_cells, w_interfaces, dw, D, Dprime, i) end fill!(production, 0.0) for i = 2:(N - 1) production[i, i + 1] = (max(0, C[i + 1]) * u_tilde[i + 1] + D[i + 1] * u[i + 1]/dw)/dw production[i, i - 1] = (-min(0, C[i]) * u_tilde[i] + D[i] * u[i - 1]/dw)/dw end production[1, 2] = (max(0, C[2]) * u_tilde[2] + D[2] * u[2]/dw)/dw production[N, N - 1] = (-min(0, C[N]) * u_tilde[N] + D[N] * u[N - 1]/dw)/dw return production end # Create Fokker-Planck problem prob = ConservativePDSProblem(P, u0, tspan, p) saveat = range(tspan..., length = 100) sol_heun = solve(prob, Heun(), dt = dt, adaptive = false, save_everystep = false, saveat = saveat) sol_mpe = solve(prob, MPE(), dt = dt, adaptive = false, save_everystep = false, saveat = saveat) sol_mprk = solve(prob, MPRK22(1.0), dt = dt, adaptive = false, save_everystep = false, saveat = saveat) @show (sum.(sol_heun.u) .- sum(sol_heun.u[1])) @show (sum.(sol_mpe.u) .- sum(sol_mpe.u[1])) @show (sum.(sol_mprk.u) .- sum(sol_mprk.u[1])) ```

In addition, MPRK22 sometimes produces spurious oscillations for large timesteps. By using the code above for MPRK(0.5) and setting dt = 0.1, we can observe oscillations around 0, which should not exist (and do not exist for MPE()). Using MPRK(1.0) with dt = 1.0 produces even worse results.

[ranocha]

The problem is that you use production from the parameters (cache) and modify it when calling P again. The out-of-place implementations assume that you can call

P1 = P(u1, p, t)
# u2 = something_computed_from_P1
P2 = P(u2, p, t)
# P1 should still be the same as before!

but you change P1 when calling P(u2, p, t). If you replace the line

fill!(production, 0.0)

by

production = zeros(N, N)

you should get the correct result (but not with optimal performance, of course).

[JoshuaLampert]

Thanks a lot for finding that @ranocha! I was putting the production matrix in the cache because of #34.

[JoshuaLampert]

And thanks for fixing #34!