Correct dimensionality reduction for whole domain Integrals

[xtalax]

fixes #221
fixes #222

[codecov]

Codecov Report

Merging #225 (fe7ac77) into master (3957101) will increase coverage by 1.68%.
The diff coverage is 75.81%.

@@            Coverage Diff             @@
##           master     #225      +/-   ##
==========================================
+ Coverage   77.29%   78.98%   +1.68%     
==========================================
  Files          38       39       +1     
  Lines        1947     1989      +42     
==========================================
+ Hits         1505     1571      +66     
+ Misses        442      418      -24     

Impacted Files	Coverage Δ
src/error_analysis.jl	26.66% <0.00%> (-2.97%)	⬇️
...ion/schemes/upwind_difference/upwind_difference.jl	88.40% <60.00%> (ø)
src/MOL_symbolic_utils.jl	65.15% <65.15%> (ø)
src/scalar_discretization.jl	77.27% <71.42%> (+1.66%)	⬆️
src/MOL_utils.jl	79.26% <91.66%> (+11.89%)	⬆️
src/discretization/generate_bc_eqs.jl	79.87% <92.85%> (+0.34%)	⬆️
...discretization/generate_finite_difference_rules.jl	95.00% <94.11%> (+2.69%)	⬆️
src/discretization/discretize_vars.jl	86.66% <100.00%> (ø)
src/discretization/schemes/WENO/WENO.jl	92.98% <100.00%> (ø)
...schemes/centered_difference/centered_difference.jl	100.00% <100.00%> (ø)
... and 10 more

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

[sdwfrost]

Thanks! My code now runs, although isn't giving the correct answer - the PDE and the ODE should be similar (leaking mass at the upper boundary aside) - again, this is likely to be my naivete in implementing the boundary conditions properly.

using DifferentialEquations
using ModelingToolkit
using MethodOfLines
using DomainSets
using OrdinaryDiffEq
using Plots

β = 0.0005
γ = 0.25
S₀ = 990.0
I₀ = 10.0
R₀ = 0.0
tmin = 0.0
tmax = 40.0
dt = 0.1
amin = 0.0
amax = 40.0
da = 0.1

@parameters t a
@variables S(..) I(..) R(..)
Dt = Differential(t)
Da = Differential(a)
Ia = Integral(a in DomainSets.ClosedInterval(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))]

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

domains = [t ∈ (tmin, tmax), a ∈ (amin, amax)]

@named pde_system = PDESystem(eqs, bcs, domains, [a, t], [S(t), I(a, t), R(t)])
discretization = MOLFiniteDifference([a=>da], t)

prob = MethodOfLines.discretize(pde_system, discretization)
sol = solve(prob, saveat = dt)

Smat = sol.u[S(t)]
Imat = sol.u[I(a, t)]
Rmat = sol.u[R(t)]

plot(sol.t, vec(Smat),label="S",xlabel="Time",ylabel="Number")
plot!(sol.t, sum(Imat,dims=1)',label="I")
plot!(sol.t, vec(Rmat),label="R")

## The above should be comparable to below
function sir_ode(u,p,t)
    (S, I, R) = u
    (β, γ) = p
    dS = -β*I*S
    dI = β*I*S - γ*I
    dR = γ*I
    [dS, dI, dR]
end;

prob_ode = ODEProblem(sir_ode, [S₀, I₀, R₀], (tmin, tmax), [β,γ])
sol_ode = solve(prob_ode, Tsit5(), saveat=dt)
plot(sol_ode)