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

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!