Problem with upper absorbing boundary with MethodOfLines.jl

Specific Domains Modelling & Simulations
1

Hi everyone, the following code runs, but my system leaks mass, due to the upper boundary - N should be constant.

using ModelingToolkit
using DifferentialEquations
using MethodOfLines
using DomainSets
using Plots

@parameters t a β γ S₀ I₀ R₀
@variables S(..) I(..) R(..)
Dt = Differential(t) 
Da = Differential(a)
tmin = 0.0       # Minimum time
tmax = 40.0      # Maximum time
dt = 0.05        # Time step
amin = 0.0       # Minimum age
amax = 10.0      # Maximum age
da = 1.0         # Age step
Ia = Integral(a in DomainSets.ClosedInterval(amin,amax))
domains = [t ∈ (tmin, tmax), a ∈ (amin, amax)]

eqs = [Dt(S(t)) ~ -β * S(t) * Ia(I(a,t)), 
       Dt(I(a,t)) + Da(I(a,t)) ~ -γ*I(a,t),
       Dt(R(t)) ~ γ*Ia(I(a,t)) + I(amax,t)]

bcs = [S(0) ~ S₀,
       I(0,t) ~ β*S(t)*Ia(I(a,t)),
       I(a,0) ~ I₀/(amax-amin),
       I(amax,t) ~ 0,
       R(0) ~ R₀]

p = Dict([β=>0.5, γ=>0.25, S₀=>0.99, I₀=>0.01, R₀=>0.0])
@named pde_system = PDESystem(eqs,
                              bcs,
                              domains,
                              [a, t],
                              [S(t), I(a, t), R(t)],
                              [β, γ, S₀, I₀, R₀];
                              defaults=p);

# Discretize using the method of lines
discretization_mol = MOLFiniteDifference([a=>da],
                                          t,
                                          approx_order=2,
                                          advection_scheme=UpwindScheme(2))
@time prob_mol = discretize(pde_system, discretization_mol)
@time sol_mol = solve(prob_mol, Rodas5(), dt=dt, saveat=0.1)

t_points = sol_mol.t
S_sol_mol = sol_mol[S(t)]
I_sol_mol = (sum(sol_mol[I(a, t)], dims=1)[1,:]).*da  # Sum over age groups
R_sol_mol = sol_mol[R(t)];
N_sol_mol = S_sol_mol .+ I_sol_mol .+ R_sol_mol;

plot(t_points, S_sol_mol, label="S", xlabel="Time", ylabel="Number")
plot!(t_points, I_sol_mol, label="I")
plot!(t_points, R_sol_mol, label="R")
plot!(t_points, N_sol_mol, label="N")

I can sort of see what’s going on - the boundary condition is setting I(amax,t) to zero, so my equation for R doesn’t see the input necessary for the conservation of mass. Does anyone have any suggestions on how to get conservation of mass over the grid?

2

Open an issue for this.