Dirac delta in boundary conditions for MethodOfLines.jl

As a basis for more complex models, I’m trying to get a simplified version of the Von Foerster equation - Wikipedia discretized and running in MethodOfLines.jl.

Here’s my model:

using ModelingToolkit
using MethodOfLines
using OrdinaryDiffEq
using DomainSets

@parameters t a μ
@variables v(..) u(..)
Dt = Differential(t)
Da = Differential(a)
p = [μ => 0.04]
tmax = 100.0
amax = 100.0
da = 0.1
dt = 0.1

pde_eqns = [Dt(u(t,a)) + Da(u(t,a)) ~ -μ*u(t,a)]
domains = [t ∈ Interval(0,tmax), a ∈ Interval(0,amax)]
function a0(a)
    if a > 0
        return 0.0
    else
        return 1.0
    end
end
@register a0(a)
bcs = [u(t,0) ~ 0.0, u(0,a) ~ a0(a)]
@named pde_sys = PDESystem(pde_eqns, bcs, domains, [t,a], [u(t,a)], p)
discretization = MOLFiniteDifference([a=>da], t, approx_order = 2, grid_align = center_align)
pde_prob = discretize(pde_sys, discretization)
pde_sol = solve(pde_prob, Tsit5(), saveat = dt)

The boundary conditions for u(0,a) are meant to be 1 when a=0 and 0 when a>0, which is what I’m trying to do with the registered function. However, when I look at the discretized problem, the initial conditions are 0 for all values of a:

julia> pde_prob.u0
1000-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0
 ⋮
 0.0
 0.0
 0.0
 0.0
 0.0

I think I understand what’s going on here - while the above function a0 is correct in the continuous sense, it doesn’t work in the discretized system. If I change function:

function a0(a)
    if a > da
        return 0.0
    else
        return 1.0
    end
end

This now does something more sensible:

julia> pde_prob.u0
1000-element Vector{Float64}:
 1.0
 0.0
 0.0
 0.0
 0.0
 ⋮
 0.0
 0.0
 0.0
 0.0
 0.0

However, now the solution gives nonsensical results - the integral over a should be an exponential decay, but it shows exponential growth, which again points to another boundary issue problem. Any tips on how to fix this? Full code, along with the ‘correct’ exponential decay as an ODE is given below.

using ModelingToolkit
using MethodOfLines
using OrdinaryDiffEq
using DomainSets

@parameters t a μ
@variables v(..) u(..)
Dt = Differential(t)
Da = Differential(a)
p = [μ => 0.04]
tmax = 100.0
amax = 100.0
da = 0.1
dt = 0.1

# ODE version
ode_eqns = [Dt(v(t)) ~ -μ*v(t)]
v0 = [v(t) => 1.0]
@named ode_sys = ODESystem(ode_eqns)
ode_prob = ODEProblem(ode_sys, v0, (0,tmax), p)
ode_sol = solve(ode_prob, Tsit5(), saveat = dt)

# PDE version
pde_eqns = [Dt(u(t,a)) + Da(u(t,a)) ~ -μ*u(t,a)]
domains = [t ∈ Interval(0,tmax), a ∈ Interval(0,amax)]
function a0(a)
    if a > da
        return 0.0
    else
        return 1.0
    end
end
@register a0(a)
bcs = [u(t,0) ~ 0.0, u(0,a) ~ a0(a)]
@named pde_sys = PDESystem(pde_eqns, bcs, domains, [t,a], [u(t,a)], p)
discretization = MOLFiniteDifference([a=>da], t, approx_order = 2, grid_align = center_align)
pde_prob = discretize(pde_sys, discretization)
pde_sol = solve(pde_prob, Tsit5(), saveat = dt)

I don’t know what’s happening here, but a suggestion:
Try

using IfElse
a0(a) = IfElse.ifelse(a>0, 0.0, 1.0)

Let me know if that works, if not I can take a closer look.

bcs = [u(t,0) ~ 0.0, u(0,a) ~ a0(a)]

Those seem contradictory to me :thinking:

1 Like

Also, a centered difference scheme is probably unstable for an advection equation. You should rather be using e.g. an upwind scheme.

That said, MOLFiniteDifference takes an upwind_order parameter, but the manual doesn’t say what its effect is.
@xtalax It appears that upwind discretization should be performed automatically…? Can you briefly comment? I’m having a hard time understanding the source :smile:

Yeah the source relies heavily on Symbolics.jl and SymbolicUtils.jl, not so easy to understand if you don’t have experience with it. To see where terms are recognised and schemes applied see this file, also useful to understand this is the documentation on @rule. Short answer, yes advection terms are automatically intercepted and the correct upwind scheme applied :). Any more questions on functionality welcome!

upwind_order is the aproximation order of the upwind scheme, note that this appears to lead to instability when it is above 1, which is its default, which I suspect is due to offside stencils near the boundaries.

Current work is on implementing flux limiters to ensure stability and higher order accuracy for advection terms in a wider set of circumstances.

1 Like

Ah another thing @sdwfrost, you want to use @register_symbolic rather than @register, there was an error in the docs - Sorry!

Thank you for pointing me in the right direction!

advection terms are automatically intercepted and the correct upwind scheme applied

Cool! That’s a neat feature, yet as far as I can tell not mentioned explicitly in the docs. Neither here How it Works · MethodOfLines.jl nor here MOLFiniteDifference · MethodOfLines.jl.

In OPs example, the discretizer seems to pick the wrong direction though:

symbolic_discretize(pde_sys, discretization)
(ODESystem(Equation[10.0u[3](t) + Differential(t)(u[2](t)) - 10.0u[2](t) ~ -μ*u[2](t),
[...]

Explains the instability. Indeed, bringing the Da term to the rhs solves the issue, i.e. this integrates fine:

pde_eqns = [Dt(u(t,a)) ~ -μ*u(t,a) - Da(u(t,a))]
domains = [t ∈ Interval(0,tmax), a ∈ Interval(0,amax)]
function a0(a)
    if da < a <= 2da # don't start right at the boundary.
        return 1.0/da
    else
        return 0.0
    end
end
@register_symbolic a0(a)
bcs = [u(t,0) ~ 0.0, u(0,a) ~ a0(a)]
@named pde_sys = PDESystem(pde_eqns, bcs, domains, [t,a], [u(t,a)], p)
discretization = MOLFiniteDifference([a=>da], t, approx_order = 2, grid_align = center_align)
pde_prob = discretize(pde_sys, discretization)

pde_sol = solve(pde_prob, Tsit5(), saveat = dt)

Cool! That’s a neat feature, yet as far as I can tell not mentioned explicitly in the docs. Neither here How it Works · MethodOfLines.jl nor here MOLFiniteDifference · MethodOfLines.jl.

You’re right, I will add this.

Explains the instability. Indeed, bringing the Da term to the rhs solves the issue, i.e. this integrates fine:

Good catch! That is indeed a big problem, pushing a hotfix as soon as possible.

1 Like

This appears to be related to advection terms appearing on the same side as the time derivative, this causes the ‘wrong’ scheme to be dispatched from the perspective of the case when that term is on the other side of the equation - both positive and negative coefficients fail in this case.

I will update the docs to reflect this limitation until such a time as we can automatically rearrange equations in to a set form pre discretization.