Partial Integro-Differential Equations, PIDEs with integration over space possible?

[BernhardAhrens]

Is there a way to specify an integro-differential equation where one term integrates over the spatial domain (basically the sum of a term from xminto x)? I have a system where matter is deposited on the upper boundary of the system, whereby it basically shifts the reference frame upwards. The deposited matter can also decompose and thereby basically make the grid boxes smaller. I am representing the deposition and the decomposition with an extra advective flux that is the sum of the deposition and the decomposition over the grid.

I tried to boil it down to an MWE with NeuralPDE.jl and MethodOfLines.jl. In the MWE, I have not implemented advection completely but focused on the bit that would represent the additional advection velocity. This advection velocity is basically the cumulative sum of changes in u(t,x) from xmin to x. In the MWE, I solve this equation with NeuralPDE.jl and MethodOfLines.jl with an Integral Ix that represents the cumulative sum cumuSum(t,x) ~ Ix(u(t,x)). Without Ix, it works well with NeuralPDE.jl and MethodOfLines.jl. With Ix, it errors for MethodOfLines with ERROR: MethodError: Cannot `convert` an object of type SymbolicUtils.Add{Float64, Float64, Dict{Any, Number}, Nothing} to an object of type Float64

I wonder if my approach of replicating the NeuralPDE.jl way to specify integro-differential equations is the right starting point? I want to take the opportunity while writing this issue to applaud the developers of MethodOfLines for a great and easy-to-use environment and big steps forward in the last weeks 👏

using NeuralPDE, Flux, MethodOfLines, ModelingToolkit, OrdinaryDiffEq, Optimization, OptimizationOptimJL, DomainSets, Plots
import ModelingToolkit: Interval, infimum, supremum

@parameters t, x
@variables u(..) cumuSum(..)
Dt = Differential(t)
Dx = Differential(x)
xmin = 0.0
xmax = 2.0

Ix = Integral(x in DomainSets.ClosedInterval(xmin, x)) # basically cumulative sum from 0 to x

eq = [
    cumuSum(t,x) ~ Ix(u(t,x)), # with Integral from xmin to x: runs with NeuralPDE.jl, but not MethodOfLines.jl
    #cumuSum(t,x) ~ u(t,x), # without Integral works with NeuralPDE.jl and MethodOfLines.jl
    Dt(u(t,x)) + 2*u(t,x) + 5*Dx(cumuSum(t,x)) ~ 1
]
bcs = [u(0.,x) ~ 0.0, Dx(u(t,0.0)) ~ 0.0, Dx(u(t,2.0)) ~ 0]

domains = [t ∈ Interval(0.0,2.0), x ∈ Interval(xmin,xmax)]

@named pde_system = PDESystem(eq,bcs,domains,[t,x],[u(t,x), cumuSum(t,x)])

# solve with with NeuralPDE.jl
chain = Chain(Dense(2,15,Flux.tanh),Dense(15,15,Flux.tanh),Dense(15,1)) |> f64
strategy_ = GridTraining(0.1)
discretization = PhysicsInformedNN(chain,strategy_)

prob = NeuralPDE.discretize(pde_system,discretization)
callback = function (p,l)
    println("Current loss is: $l")
    return false
end
res = Optimization.solve(prob, BFGS(); callback = callback, maxiters=100)

ts,xs = [infimum(d.domain):0.1:supremum(d.domain) for d in domains]
phi = discretization.phi
u_predict  = [[first(phi([t,x],res.u)) for t in ts] for x in xs]

p1 = plot(ts, u_predict, legend = :inline)

# solve with MethodOfLines
order = 2
discretization = MOLFiniteDifference([x=>0.1], t, approx_order=order, grid_align=center_align)

prob = MethodOfLines.discretize(pde_system, discretization)
sol = solve(prob, QNDF(), saveat=0.1);

sol.u

p2 = plot(sol[u(t,x)], legend = :inline)

plot(p1,p2)

[xtalax]

With MOL, this method of defining an aux variable to represent the argument to the derivative being an integral is the right approach, though I'm currently working on a system that does this transformation automatically.

Internally, MOL represents the discretized form of each variable as a grid of symbolic variables, and the equations for each of these states are built up of expressions in terms of these by replacing derivatives by finite difference expressions, depending on certain replacement rules generated for each interior point.

It is easy to imagine a replacement rule that would replace integrals with their constant piecewise sum approximation, or some higher order quadrature scheme.

The subtlety enters this problem in making sure that the rules are general enough to handle integrals Integral(x in Interval(a, b)) from any a to any b, including both symbols and concrete values both on and off grid. It won't be possible to generate a specific rule for each combination in general.

We could limit integrals to start to be only from the lower/upper end of the domain to x to simplify to start.

I expect that the easiest way to solve this in general would be to find a way to use SymbolicUtils.@rule to extract the limits of integration, programmatically determine if they are sane, and then generate an appropriate discretization.

If this turns out not to be possible, then we will need to recursively parse the expression tree for each equation and find any integrals, and build targeted closures which will return the correct replacement output for each index, then use these to generate the rules on the interior. [NB. It may be faster just to do this for all terms, I will evaluate this].

In any case, I don't see a way for either of the limits to be a parameter outside of creating a lazy system object that defers part of the discretization, but this is a significant amount of work and somewhat defeats the point of having remakable parameters - open to comments on this.

To get an idea of what a simple replacement rule looks like see here.

If you have any questions, ideas or need anything clarified, ask away :). No stupid questions, this is a somewhat different programming paradigm so it can take a bit to get your head around it.

See also my JuliaCon talk on the package, part explaining some of the internals from 20:03.

[sdwfrost]

Is this a repeat of https://github.com/SciML/MethodOfLines.jl/issues/94 ?

[BernhardAhrens]

Yes, #94 sounds very similar. Are you also interested in integration over space?

[xtalax]

It is, but it adds more detail so I will leave it up

[Qfl3x]

I've started work on it, but it seems like it's not parsing properly:

function euler_integral(II, s, b, jx, u, ufunc) #where {T,N,Wind,DX<:Number}
    j, x = jx
    ndims(u, s) != 1 && return Num(0)
    # unit index in direction of the derivative
    I1 = unitindex(ndims(u, s), j)
    # dx for multiplication
    dx = s.dxs[x]

    Itap = [II - I1 * i for i = 0:(II[j]-2)] # Sum from current position till CartesianIndex 2

    weights = [dx for i = 0:(II[j]-2)]
    @info weights
    return dot(weights, ufunc(u, Itap, x))
end

@inline function generate_integration_rules(II::CartesianIndex, s::DiscreteSpace, depvars, pmap, indexmap, terms)
    ufunc(u, I, x) = s.discvars[u][I]
    @info "Generating"
    return reduce(vcat, [vec([@rule Integral(x in DomainSets.ClosedInterval(s.grid[x][1], x))(u) => euler_integral(Idx(II, s, u, indexmap), s, pmap.map[operation(u)][x], (x2i(s, u, x), x), u, ufunc)
                              for x in params(u, s)])
                                for u in depvars])

end

The resulting error message indicates that the function isn't being called:

A scheme has been incorrectly applied to the following equation: Differential(t)(u(x, t)) - Differential(x)(u(x, t)) - 3.0u(x, t) ~ 0.

The following rules were constructed at index CartesianIndex(32,):
7-element Vector{Any}:
 -Differential(x)(u(x, t)) => 32.0(u(t))[32] - 32.0(u(t))[33]
  Differential(x)(u(x, t)) => 32.0(u(t))[32] - 32.0(u(t))[31]
                            (Integral(x in DomainSets.ClosedInterval((s.grid[x])[1], x)))(u) => euler_integral(Idx(II, s, u, indexmap), s, (pmap.map[operation(u)])[x], (x2i(s, u, x), x), u, ufunc)
                           ⋮
                            Pair{SymbolicUtils.Symbolic{Real}, Num}(u(x, t), (u(t))[32])
                            Pair{SymbolicUtils.Symbolic{Real}, Num}(x, 0.96875)

I always get "generating", but never the weights. Unless I explicitly call the function myself. Calling either functions seems to work, but using generate_finite_difference_rules doesn't.

[xtalax]

Looks like you need to parse the rules in to pairs, see here how its done:

function generate_central_difference_rules(II::CartesianIndex, s::DiscreteSpace, terms::Vector{<:Term}, indexmap::Dict)
    rules = [[@rule Differential(x)(u) => second_order_central_difference(Idx(II, s, u, indexmap), s, u, x) for x in params(u, s)] for u in depvars]

    rules = reduce(vcat, rules)

    # Parse the rules in to pairs that can be used with `substitute`, this can be copy pasted.
    rule_pairs = []
    for t in terms
        for r in rules
            if r(t) !== nothing
                push!(rule_pairs, t => r(t))
            end
        end
    end
    return rule_pairs
end

EDIT: In fact since you aren't using rule syntax, you can remove the @rule tag and this should work. you only need this if you are using slot variables like ~a, see some of the other rulegen for examples

[Qfl3x]

An update, after more debugging, I found a bug in the Integral expression:

Integral(x in DomainSets.ClosedInterval(0.0, x))(u)

It just gives this error:

LoadError: promotion of types Float64 and Sym{Real, Base.ImmutableDict{DataType, Any}} failed to change any arguments

For this one, I still have no idea how to solve it.

I also modified the generate_finite_difference_rules code to accomodate for the II[j] == 1 scenario. This was the only way to avoid the function returning a Float.

[Qfl3x]

Update 2: Issue is specifically in the definition of the Domain. The implementation of x confuses DomainSets apparently.

[xtalax]

I bet you need to wrap it with Num(x)

[Qfl3x]

Tried wrapping with Num, same issue. Wrapping did Nothing as Sym appears to be a child type of Num

Nevermind, made a mistake when saving. Now it works. New problem that occured is the pairs code stopped working, with error code LoadError: MethodError: objects of type Pair{Term{Real, Nothing}, Num} are not callable.

The type of the Integral expression is Term{Real, Nothing}, while the type of the usual Differential expression is Term{SymbolicUtils.FnType{Tuple, Real}, Nothing}.

Anyways, if I remove the pairs loop (it seems it was removed from centered_difference), the program now generates all the rules! But then crashes at the solution:

LoadError: MethodError: Cannot `convert` an object of type Term{Float64, Nothing} to an object of type Float64

Full stacktrace: https://pastebin.com/TpRrLy40 as this stacktrace touches mostly on packages other than this one in the SciML ecosystem.

Line 31 in question:

sol = solve(prob, QNDF(), saveat=0.1);

[xtalax]

Yes, you only need to parse the rules if they are defined with @rule, pairs work without.

[xtalax]

I think I've seen this before, have a look at the system (you can find it in prob.f.sys.eqs) and check that there isn't anything that shouldn't be there i.e. no references to spatial dimensions and no integral terms - it may help to post this (low point count) so I can look over this. Getting caught out on types is one of the most common symbolics gotchas though.

Can you two @BernhardAhrens @Qfl3x please post PRs so there's no repeated work and we can collaborate easier?

[BernhardAhrens]

I have asked @Qfl3x if he could work on this issue a bit - I started initially but then did not have time to work further on it. From what I have seen on @Qfl3x's GitHub, his work is much more progressed, so it is probably best to go with his attempt at it.

[Qfl3x]

Odd thing happened now, it works (?!?); Code used to test it:

using MethodOfLines, ModelingToolkit, OrdinaryDiffEq, DomainSets, Plots
import ModelingToolkit: Interval, infimum, supremum

using SymbolicUtils

@parameters t, x
@variables u(..) cumuSum(..)
Dt = Differential(t)
Dx = Differential(x)
xmin = 0.0
xmax = 2.0

Ix = Integral(x in DomainSets.ClosedInterval(xmin, x)) # basically cumulative sum from 0 to x

eq = [
    cumuSum(t, x) ~ Ix(u(t, x)), # with Integral from xmin to x: runs with NeuralPDE.jl, but not MethodOfLines.jl
    #cumuSum(t, x) ~ u(t, x), # without Integral works with NeuralPDE.jl and MethodOfLines.jl
    Dt(u(t, x)) + 2 * u(t, x) + 5 * Dx(cumuSum(t, x)) ~ 1
]
bcs = [u(0.0, x) ~ 0.0, Dx(u(t, 0.0)) ~ 0.0, Dx(u(t, 2.0)) ~ 0]

domains = [t ∈ Interval(0.0, 2.0), x ∈ Interval(xmin, xmax)]

@named pde_system = PDESystem(eq, bcs, domains, [t, x], [u(t, x), cumuSum(t, x)])

order = 2
discretization = MOLFiniteDifference([x => 0.1], t, approx_order=order, grid_align=center_align)

prob = MethodOfLines.discretize(pde_system, discretization)
sol = solve(prob, QNDF(), saveat=0.1);

My fork: https://github.com/Qfl3x/MethodOfLines.jl/ works for this. Still need more testing though, will work on that now.

[xtalax]

try write a couple of tests using analytic integrals, i.e. remove u(t,x) and integrate t*cos(x) or something, and an analytic PIDE if you know one - look at the tests for examples. But feel free to submit the PR and push changes to the branch so I can more easily track your changes and review.

[Qfl3x]

I'm trying the following code:

using MethodOfLines, ModelingToolkit, OrdinaryDiffEq, DomainSets, Plots
import ModelingToolkit: Interval, infimum, supremum

using SymbolicUtils

@parameters t, x
@variables u(..) cumuSum(..)
Dt = Differential(t)
Dx = Differential(x)
xmin = 0.0
xmax = 2.0

Ix = Integral(x in DomainSets.ClosedInterval(xmin, x)) # basically cumulative sum from 0 to x

eq = [
    cumuSum(t, x) ~ Ix(t * cos(x)), 
    Dt(u(t, x)) + 2 * u(t, x) + 5 * Dx(cumuSum(t, x)) ~ 1
]
bcs = [u(0.0, x) ~ 0.0, Dx(u(t, 0.0)) ~ 0.0, Dx(u(t, 2.0)) ~ 0]

domains = [t ∈ Interval(0.0, 2.0), x ∈ Interval(xmin, xmax)]

@named pde_system = PDESystem(eq, bcs, domains, [t, x], [u(t, x), cumuSum(t, x)])

order = 2
discretization = MOLFiniteDifference([x => 0.1], t, approx_order=order, grid_align=center_align)

prob = MethodOfLines.discretize(pde_system, discretization)
#sol = solve(prob, QNDF(), saveat=0.1);

If I replace the Integration Ix with a Differentiation Dx, prob.f.sys.eqs doesn't change. Yet with a differentiation it works, and with an integration it crashes with:

ERROR: MethodError: Cannot `convert` an object of type SymbolicUtils.Add{Float64, Float64, Dict{Any, Number}, Nothing} to an object of type Float6

[xtalax]

Yes, there isn't any system transformation for integrals yet, so that composite argument won't work unless you wrap it in an aux variable.

Try:

using MethodOfLines, ModelingToolkit, DomainSets, OrdinaryDiffEq

# test integrals
@parameters t, x
@variables integrand(..) cumuSum(..)
Dt = Differential(t)
Dx = Differential(x)
xmin = 0.0
xmax = 2.0 * pi

Ix = Integral(x in DomainSets.ClosedInterval(xmin, x)) # basically cumulative sum from 0 to x

eqs = [cumuSum(t, x) ~ Ix(integrand(t, x))
    integrand(t, x) ~ t * cos(x)]

bcs = [cumuSum(t, xmin) ~ 0.0,
    cumuSum(t, xmax) ~ 0.0,
    Dx(integrand(t, xmin)) ~ 0.0,
    Dx(integrand(t, xmax)) ~ 0.0]

domains = [t ∈ Interval(0.0, 1.0),
    x ∈ Interval(xmin, xmax)]

@named pde_system = PDESystem(eq, bcs, domains, [t, x], [integrand(t, x), cumuSum(t, x)])

asf(t, x) = t*sin(x)

disc = MOLFiniteDifference([x => 20], t)

prob = discretize(pde_system, disc)

sol = solve(prob, Tsit5(), saveat=0.1)

xdisc = sol[x]
tdisc = sol[t]

cumuSumsol = sol[cumuSum(t, x)]

exact = [asf(t_, x_) for t_ in tdisc, x_ in xdisc]

@test cumuSumsol ≈ exact atol=1e-3

Then try to remove the bcs, they shouldn't be needed

[xtalax]

If you can comprehend descend_to_incompatible in src/system_parsing/pde_system_transformation.jl, make an edit to check for incompatible integral arguments and do this wrapping automatically. If not, I can do it, just create the PR and I'll push this to it

[Qfl3x]

I'm afraid there is still an issue. I've found an Integro-Differential Equation that has an analytical solution:

2u(x) + Dx(u(x)) + Integral(x, 0..x)(u(x)) = 1

This equation has a solution for u(0)=0 of u(x) = exp(-x)*x

However, if I input it I still got the error from earlier:

ERROR: MethodError: Cannot `convert` an object of type SymbolicUtils.Add{Float64, Float64, Dict{Any, Number}, Nothing} to an object of type Float64