PDE parameters in MOLFiniteDifference

[michaelsachs]

Hi, many thanks for the recent fix with which MOLFiniteDifference now works with multiple dependent variables. I have noticed that I am encountering the error below when the PDESystem has parameters in it, here [Dn,Dp].

using OrdinaryDiffEq, ModelingToolkit, DiffEqOperators

# Parameters, variables, and derivatives
@parameters t x Dn Dp
@variables n(..) p(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

# 1D PDE and boundary conditions
eqs  = [Dt(n(t,x)) ~ Dn*Dxx(n(t,x)),
        Dt(p(t,x)) ~ Dp*Dxx(p(t,x))]

bcs = [n(0,x) ~ 2x,
        Dx(n(t,0)) ~ 0,
        Dx(n(t,10)) ~ 0,
        p(0,x) ~ 2x,
        Dx(p(t,0)) ~ 0,
        Dx(p(t,10)) ~ 0]

# Space and time domains
domains = [t ∈ IntervalDomain(0.0,100.0),
           x ∈ IntervalDomain(0.0,10.0)]

# PDE system
pdesys = PDESystem(eqs,bcs,domains,[t,x],[n(t,x),p(t,x)],[Dn,Dp])

# Method of lines discretization
Δx = 0.1
discretization = MOLFiniteDifference([x=>Δx],t)

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization) 

Errors with this stacktrace:

ERROR: LoadError: ArgumentError: Dict(kv): kv needs to be an iterator of tuples or pairs
Stacktrace:
 [1] Dict(kv::Vector{Any})
   @ Base .\dict.jl:132
 [2] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{Vector{Pair{Num, Float64}}, Num})
   @ DiffEqOperators ~\.julia\packages\DiffEqOperators\skiw0\src\MOLFiniteDifference\MOL_discretization.jl:153
 [3] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{Vector{Pair{Num, Float64}}, Num})
   @ DiffEqOperators ~\.julia\packages\DiffEqOperators\skiw0\src\MOLFiniteDifference\MOL_discretization.jl:158

Does this mean that defaults need to be supplied for the parameters in the PDESystem when discretizing using MOLFiniteDifference? I have tried to do that as follows but this gives me another error:

def = [Dn => 1.0, Dp => 1.0]
pdesys = PDESystem(eqs,bcs,domains,[t,x],[n(t,x),p(t,x)],[Dn,Dp],def)

# Method of lines discretization
Δx = 0.1
discretization = MOLFiniteDifference([x=>Δx],t)

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization) 

ERROR: LoadError: MethodError: Cannot `convert` an object of type 
  Vector{Pair{Num, Float64}} to an object of type
  Dict
Closest candidates are:
  convert(::Type{T}, ::Intervals.Interval{T, L, R} where {L<:Intervals.Bound, R<:Intervals.Bound}) where T at C:\Users\Michael\.julia\packages\Intervals\ua9cq\src\interval.jl:253
  convert(::Type{T}, ::Intervals.AnchoredInterval{P, T, L, R} where {L<:Intervals.Bounded, R<:Intervals.Bounded}) where {P, T} at C:\Users\Michael\.julia\packages\Intervals\ua9cq\src\anchoredinterval.jl:181
  convert(::Type{var"#s4"} where var"#s4"<:Dict, ::ChainRulesCore.Composite{var"#s3", var"#s2"} where {var"#s3"<:Dict, var"#s2"<:Dict}) at C:\Users\Michael\.julia\packages\ChainRulesCore\D0go7\src\differentials\composite.jl:69
  ...
Stacktrace:
 [1] PDESystem(eqs::Vector{Equation}, bcs::Vector{Equation}, domain::Vector{ModelingToolkit.VarDomainPairing}, indvars::Vector{Num}, depvars::Vector{Num}, ps::Vector{Num}, defaults::Vector{Pair{Num, Float64}})
   @ ModelingToolkit ~\.julia\packages\ModelingToolkit\SVti5\src\systems\pde\pdesystem.jl:56

[ChrisRackauckas]

Yes, it hasn't been implemented yet.

[michaelsachs]

Alright, thanks for the quick reply!

[valentinsulzer]

Do we want to explicitly pass in Dn and Dp or should they be parsed from the equations?

[ChrisRackauckas]

They should be handled similarly to the PINNs, https://neuralpde.sciml.ai/dev/pinn/parm_estim/ where defaults should be required on them. For MOL we would then instead generate an ODEProblem with a parameter vector p with those default values and using the parameter vector p inside the generated equations (the reason for this is to then be compatible with automatic differentiation and optimization in the generated code).

[valentinsulzer]

Makes sense, is there anything blocking this or does it just need to be done? I can give it a go

[ChrisRackauckas]

It just needs to get done.

[valentinsulzer]

This already works fine:

pdesys = PDESystem(eqs,bcs,domains,[t,x],[n(t,x),p(t,x)],[Dn=>1,Dp=>1])

I'll add an example

[ChrisRackauckas]

oh interesting haha. Yeah I put a few things in there for it but I didn't think I completed it.

[valentinsulzer]

there's nothing like an issue that resolves itself on its own eh